.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/continuum_mechanics/tensor_basis.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_continuum_mechanics_tensor_basis.py: The basis of a tensor: frames, convected bases and the metric ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The components of a ``Tensor2`` / ``Tensor4`` only mean something together with the basis they are written in. By default this is the fixed orthonormal lab basis: ``t.basis`` is ``None`` and nothing is stored. A ``simcoon.Basis`` makes another choice explicit: * an **orthonormal** basis, built from a rotation (a material frame, a corotational frame). Its metric is the identity; * a **natural** basis, any three independent vectors, typically the convected basis ``g_i = F G_i`` of a deformation gradient. Its metric ``g = A^T A`` enters the invariants. Three operations, three meanings: * ``t.to_basis(b)``: the same tensor, components re-expressed in ``b``; * ``t.with_basis(b)``: the same components in another basis, a transport; * ``t.rotate(R, active=False)``: the same tensor in the frame turned by ``R``. The type tag gives the variance of the components: stress contravariant, strain covariant, stiffness and compliance likewise on four indices, concentration tensors mixed. .. GENERATED FROM PYTHON SOURCE LINES 26-35 .. code-block:: Python import copy import pickle import numpy as np import simcoon as sim np.set_printoptions(precision=4, suppress=True) .. GENERATED FROM PYTHON SOURCE LINES 36-39 1. The lab basis is the default --------------------------------- Nothing changes for a tensor that never leaves the lab: ``basis`` is ``None``. .. GENERATED FROM PYTHON SOURCE LINES 39-47 .. code-block:: Python sigma = sim.Tensor2.stress(np.array([100.0, 50.0, 0.0, 30.0, 0.0, 0.0])) eps = sim.Tensor2.strain(np.array([0.010, -0.003, -0.003, 0.005, 0.002, 0.001])) L_iso = sim.Tensor4.stiffness(sim.L_iso([70000.0, 0.3], "Enu")) print("sigma.basis:", sigma.basis) print("(L_iso @ eps).basis:", (L_iso @ eps).basis) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma.basis: None (L_iso @ eps).basis: None .. GENERATED FROM PYTHON SOURCE LINES 48-53 2. An orthonormal basis and its description ---------------------------------------------- A ``Basis`` holds the basis vectors as the columns of ``matrix`` (their lab components). The ``name`` is a label for printing; two bases are the same when their vectors are. .. GENERATED FROM PYTHON SOURCE LINES 53-66 .. code-block:: Python rot = sim.Rotation.from_euler("zxz", [30.0, 20.0, 10.0], degrees=True) material = sim.Basis(rotation=rot, name="material") print(material) print("orthonormal:", material.orthonormal, "| single:", material.single, "| name:", material.name) print("rotation is kept:", material.rotation.equals(rot)) print("matrix (columns = basis vectors):\n", material.matrix) print("metric is the identity:", np.allclose(material.metric, np.eye(3))) print("same vectors, other name -> equal:", material.equals(sim.Basis(rotation=rot, name="other"))) .. rst-class:: sphx-glr-script-out .. code-block:: none Basis(orthonormal, name='material') orthonormal: True | single: True | name: material rotation is kept: True matrix (columns = basis vectors): [[ 0.7713 -0.6337 0.0594] [ 0.6131 0.7146 -0.3368] [ 0.171 0.2962 0.9397]] metric is the identity: True same vectors, other name -> equal: True .. GENERATED FROM PYTHON SOURCE LINES 67-71 3. ``to_basis``: the same tensor in another basis ---------------------------------------------------- The components change, the tensor does not: its invariants are those of the lab, and ``to_basis(None)`` gives the lab components back. .. GENERATED FROM PYTHON SOURCE LINES 71-84 .. code-block:: Python sigma_m = sigma.to_basis(material) print(sigma_m) print("material components:", sigma_m.voigt) print("same Mises:", np.isclose(sigma_m.mises(), sigma.mises()), "| same trace:", np.isclose(sigma_m.trace(), sigma.trace())) print("back to the lab:", np.allclose(sigma_m.to_basis(None).voigt, sigma.voigt)) # Between two frames directly other = sim.Basis(rotation=sim.Rotation.from_euler("z", 45.0, degrees=True)) print("material -> other == lab -> other:", np.allclose(sigma_m.to_basis(other).voigt, sigma.to_basis(other).voigt)) .. rst-class:: sphx-glr-script-out .. code-block:: none Tensor2(type='stress', basis=Basis(orthonormal, name='material')) material components: [106.6534 38.5216 4.825 -22.0923 -12.4457 -8.1218] same Mises: True | same trace: True back to the lab: True material -> other == lab -> other: True .. GENERATED FROM PYTHON SOURCE LINES 85-89 4. Tensors in different bases cannot be mixed ------------------------------------------------ Adding material-frame components to lab components is meaningless; it is refused instead of returning numbers. .. GENERATED FROM PYTHON SOURCE LINES 89-101 .. code-block:: Python for label, operation in [("sigma + sigma_m", lambda: sigma + sigma_m), ("L_iso @ eps_m", lambda: L_iso @ eps.to_basis(material)), ("sigma % sigma_m", lambda: sigma % sigma_m)]: try: operation() except ValueError as err: print(f"{label}: {str(err).split(' (')[0]}") print("comparison across bases is False:", sigma == sigma_m) print("in a common basis it works:", (sigma_m + sigma.to_basis(material)).basis) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma + sigma_m: Mixed basis: lab vs Basis(orthonormal, name='material') L_iso @ eps_m: Mixed basis: lab vs Basis(orthonormal, name='material') sigma % sigma_m: Mixed basis: lab vs Basis(orthonormal, name='material') comparison across bases is False: False in a common basis it works: Basis(orthonormal, name='material') .. GENERATED FROM PYTHON SOURCE LINES 102-107 5. A passive rotation ships its frame ---------------------------------------- ``rotate(R, active=False)`` returns the same tensor seen from the frame turned by ``R``. The result now remembers that frame; a second passive rotation is read in the current frame, so the frames compose. .. GENERATED FROM PYTHON SOURCE LINES 107-117 .. code-block:: Python sigma_p = sigma.rotate(rot, active=False) print("basis after a passive rotation:", sigma_p.basis) print("same numbers as to_basis:", np.allclose(sigma_p.voigt, sigma_m.voigt)) rot2 = sim.Rotation.from_euler("x", 35.0, degrees=True) sigma_pp = sigma_p.rotate(rot2, active=False) print("composed frame = rot * rot2:", sigma_pp.basis.rotation.equals(rot * rot2)) print("still the same tensor:", np.allclose(sigma_pp.to_basis(None).voigt, sigma.voigt)) .. rst-class:: sphx-glr-script-out .. code-block:: none basis after a passive rotation: Basis(orthonormal) same numbers as to_basis: True composed frame = rot * rot2: True still the same tensor: True .. GENERATED FROM PYTHON SOURCE LINES 118-123 6. An active rotation is a transport --------------------------------------- ``rotate(R)`` gives another tensor, ``Q X Q^T``. A lab tensor stays written in the lab. A tensor with its own basis keeps its components and has its basis vectors turned: the same result, without touching the numbers. .. GENERATED FROM PYTHON SOURCE LINES 123-133 .. code-block:: Python sigma_rot = sigma.rotate(rot2) print("lab tensor, active rotation -> basis:", sigma_rot.basis) turned = sigma_m.rotate(rot2) print("framed tensor: components kept:", np.array_equal(turned.voigt, sigma_m.voigt)) print("its basis is rot2 * rot:", turned.basis.rotation.equals(rot2 * rot)) print("same tensor as the lab rotation:", np.allclose(turned.to_basis(None).voigt, sigma_rot.voigt)) .. rst-class:: sphx-glr-script-out .. code-block:: none lab tensor, active rotation -> basis: None framed tensor: components kept: True its basis is rot2 * rot: True same tensor as the lab rotation: True .. GENERATED FROM PYTHON SOURCE LINES 134-138 7. A constitutive law written in its material frame ------------------------------------------------------- A cubic stiffness is known in the crystal frame. Tag it with that basis, bring the strain to the same basis, contract, and take the result back. .. GENERATED FROM PYTHON SOURCE LINES 138-149 .. code-block:: Python L_crystal = sim.Tensor4.stiffness( sim.L_cubic([200000.0, 0.3, 80000.0], "EnuG")).with_basis(material) sigma_crystal = L_crystal @ eps.to_basis(material) print("stress is in the crystal frame:", sigma_crystal.basis) L_lab = L_crystal.to_basis(None) # the same tensor, lab components print("same stress as the lab computation:", np.allclose(sigma_crystal.to_basis(None).voigt, (L_lab @ eps).voigt)) print("compliance keeps the basis:", L_crystal.inverse().basis is material) .. rst-class:: sphx-glr-script-out .. code-block:: none stress is in the crystal frame: Basis(orthonormal, name='material') same stress as the lab computation: True compliance keeps the basis: True .. GENERATED FROM PYTHON SOURCE LINES 150-155 8. Batches: one basis for all, or one per tensor --------------------------------------------------- A single basis serves a whole batch: it is stored once and shared by reference by every tensor derived from it. A batch of N rotations gives one basis per tensor (one orientation per grain or per Gauss point). .. GENERATED FROM PYTHON SOURCE LINES 155-178 .. code-block:: Python N = 1000 rng = np.random.default_rng(42) sigma_batch = sim.Tensor2.stress(rng.standard_normal((N, 6)) * 100.0) shared = sigma_batch.to_basis(material) print("one object for", N, "tensors:", shared.basis is material, "| also after indexing:", shared[17].basis is material, "| and after arithmetic:", (2.0 * shared - shared).basis is material) grains = sim.Basis(rotation=sim.Rotation.random(N, random_state=1), name="grains") per_grain = sigma_batch.to_basis(grains) print(per_grain) print("basis of tensor 17:", per_grain[17].basis) print("round trip error:", np.max(np.abs(per_grain.to_basis(None).voigt - sigma_batch.voigt))) # Stacking: equal bases stay shared, different ones become a per-tensor basis a = sigma.to_basis(material) b = sigma.to_basis(other) print("from_list, same basis:", sim.Tensor2.from_list([a, a]).basis) print("from_list, two bases:", sim.Tensor2.from_list([a, b]).basis) print("concatenate:", sim.Tensor2.concatenate([shared[:3], shared[10:12]]).basis) .. rst-class:: sphx-glr-script-out .. code-block:: none one object for 1000 tensors: True | also after indexing: True | and after arithmetic: True Tensor2(N=1000, type='stress', basis=Basis(orthonormal, N=1000, name='grains')) basis of tensor 17: Basis(orthonormal, name='grains') round trip error: 4.547473508864641e-13 from_list, same basis: Basis(orthonormal, name='material') from_list, two bases: Basis(orthonormal, N=2) concatenate: Basis(orthonormal, name='material') .. GENERATED FROM PYTHON SOURCE LINES 179-184 9. A natural basis and its metric ------------------------------------ Simple shear ``F = I + gamma e1 (x) e2``. The convected basis ``g_i = F e_i`` is not orthonormal: ``g_2 = (gamma, 1, 0)``. Its metric is the right Cauchy-Green tensor ``C = F^T F``. ``from_F`` keeps a reference to ``F``. .. GENERATED FROM PYTHON SOURCE LINES 184-201 .. code-block:: Python gamma = 0.5 F = np.array([[1.0, gamma, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) convected = sim.Basis.from_F(F, name="convected") print(convected, "| orthonormal:", convected.orthonormal) print("metric g_ij = C:\n", convected.metric) print("inverse metric g^ij:\n", convected.inverse_metric) print("volume det(A) = J:", convected.det) print("reciprocal vectors g^i (columns of A^-T):\n", convected.reciprocal) print("right stretch U = sqrt(g):\n", convected.stretch) print("closest orthonormal basis, R = A U^-1:", convected.polar) print(convected.polar.matrix) print("any three vectors work too:", sim.Basis(vectors=np.diag([2.0, 1.0, 0.5])).metric.diagonal()) .. rst-class:: sphx-glr-script-out .. code-block:: none Basis(natural, name='convected') | orthonormal: False metric g_ij = C: [[1. 0.5 0. ] [0.5 1.25 0. ] [0. 0. 1. ]] inverse metric g^ij: [[ 1.25 -0.5 0. ] [-0.5 1. 0. ] [ 0. 0. 1. ]] volume det(A) = J: 1.0 reciprocal vectors g^i (columns of A^-T): [[ 1. 0. 0. ] [-0.5 1. 0. ] [ 0. 0. 1. ]] right stretch U = sqrt(g): [[0.9701 0.2425 0. ] [0.2425 1.0914 0. ] [0. 0. 1. ]] closest orthonormal basis, R = A U^-1: Basis(orthonormal, name='convected') [[ 0.9701 0.2425 0. ] [-0.2425 0.9701 -0. ] [-0. 0. 1. ]] any three vectors work too: [4. 1. 0.25] .. GENERATED FROM PYTHON SOURCE LINES 202-207 10. ``with_basis``: transport by keeping the components ----------------------------------------------------------- The Kirchhoff stress has, in the convected basis, the components of the second Piola-Kirchhoff stress: ``tau^ij = S^IJ``. Likewise the Almansi strain has the components of the Green-Lagrange strain: ``e_ij = E_IJ``. .. GENERATED FROM PYTHON SOURCE LINES 207-220 .. code-block:: Python S = sim.Tensor2.stress(np.array([0.0, 80.0, 0.0, 0.0, 0.0, 0.0])) # PK2 = 80 e2 (x) e2 tau = S.with_basis(convected) print("same numbers:", np.array_equal(tau.voigt, S.voigt), "| basis:", tau.basis) print("lab components of tau:\n", tau.to_basis(None).mat) print("== F S F^T:", np.allclose(tau.to_basis(None).mat, F @ S.mat @ F.T)) E = sim.Tensor2.strain(0.5 * (F.T @ F - np.eye(3))) # Green-Lagrange almansi = E.with_basis(convected).to_basis(None) print("Almansi from the convected components:\n", almansi.mat) print("== 1/2 (I - b^-1):", np.allclose(almansi.mat, 0.5 * (np.eye(3) - np.linalg.inv(F @ F.T)))) .. rst-class:: sphx-glr-script-out .. code-block:: none same numbers: True | basis: Basis(natural, name='convected') lab components of tau: [[20. 40. 0.] [40. 80. 0.] [ 0. 0. 0.]] == F S F^T: True Almansi from the convected components: [[ 0. 0.25 0. ] [ 0.25 -0.125 0. ] [ 0. 0. 0. ]] == 1/2 (I - b^-1): True .. GENERATED FROM PYTHON SOURCE LINES 221-225 11. Invariants with the metric --------------------------------- In a natural basis the trace is ``g_ij tau^ij``, not the sum of the diagonal components. Deviator, norm, von Mises, determinant and eigenvalues follow. .. GENERATED FROM PYTHON SOURCE LINES 225-249 .. code-block:: Python tau_lab = tau.to_basis(None) print("trace: ", tau.trace(), "| naive diagonal sum:", tau.voigt[:3].sum(), "| lab:", tau_lab.trace()) print("von Mises: ", tau.mises(), "| lab:", tau_lab.mises()) print("norm: ", tau.norm(), "| lab:", tau_lab.norm()) # Determinant and principal values, on a full stress in a basis that also changes volume: # det T = det(T^ij) det(g), and the eigenvalues solve det(T^ij - lambda g^ij) = 0. F_gen = np.array([[1.20, 0.15, -0.05], [0.10, 0.90, 0.20], [0.00, -0.10, 1.10]]) general = sim.Basis.from_F(F_gen) full = sim.Tensor2.stress(np.array([100.0, 50.0, 20.0, 30.0, -10.0, 5.0])) full_g = full.to_basis(general) print("det: ", round(full_g.det(), 3), "| lab:", round(full.det(), 3), "| naive det of the components:", round(np.linalg.det(full_g.mat), 3)) print("eigenvalues:", full_g.eigvals(), "| lab:", np.linalg.eigvalsh(full.mat)) s = tau.dev() # tau - tr/3 g^-1, same basis print("deviator is traceless:", np.isclose(s.trace(), 0.0), "| basis kept:", s.basis is convected) print("deviator, lab components match:", np.allclose(s.to_basis(None).voigt, tau_lab.dev().voigt)) .. rst-class:: sphx-glr-script-out .. code-block:: none trace: 100.0 | naive diagonal sum: 80.0 | lab: 100.0 von Mises: 100.0 | lab: 100.0 norm: 100.0 | lab: 100.0 det: 71500.0 | lab: 71500.0 | naive det of the components: 49985.459 eigenvalues: [ 15.7109 39.7239 114.5652] | lab: [ 15.7109 39.7239 114.5652] deviator is traceless: True | basis kept: True deviator, lab components match: True .. GENERATED FROM PYTHON SOURCE LINES 250-254 12. Contractions: dual variances need no metric -------------------------------------------------- ``sigma : eps = sigma^ij eps_ij`` acts on the components in any basis. Two tensors of the same variance need the metric on both indices. .. GENERATED FROM PYTHON SOURCE LINES 254-264 .. code-block:: Python e = eps.to_basis(convected) # covariant components t = sigma.to_basis(convected) # contravariant components print(f"sigma : eps, lab: {sigma % eps:.6f}") print(f"sigma : eps, convected: {t % e:.6f} (plain sum of products)") print(f"sigma : sigma, lab: {sigma % sigma:.3f}") print(f"sigma : sigma, convected: {t % t:.3f} (with the metric)") print(f"naive sum of squares: {np.sum(t.mat * t.mat):.3f} (not an invariant)") print("double_contract:", sim.double_contract(t, e)) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma : eps, lab: 1.000000 sigma : eps, convected: 1.000000 (plain sum of products) sigma : sigma, lab: 14300.000 sigma : sigma, convected: 14300.000 (with the metric) naive sum of squares: 9356.250 (not an invariant) double_contract: [1.] .. GENERATED FROM PYTHON SOURCE LINES 265-269 13. The identity tensor is the metric ---------------------------------------- Written with contravariant components the identity is ``g^ij``; with covariant components it is ``g_ij``. .. GENERATED FROM PYTHON SOURCE LINES 269-276 .. code-block:: Python print("identity, stress-typed == g^ij:", np.allclose(sim.Tensor2.identity("stress", basis=convected).mat, convected.inverse_metric)) print("identity, strain-typed == g_ij:", np.allclose(sim.Tensor2.identity("strain", basis=convected).mat, convected.metric)) print("in an orthonormal basis:", sim.Tensor2.identity("stress", basis=material).voigt) .. rst-class:: sphx-glr-script-out .. code-block:: none identity, stress-typed == g^ij: True identity, strain-typed == g_ij: True in an orthonormal basis: [1. 1. 1. 0. 0. 0.] .. GENERATED FROM PYTHON SOURCE LINES 277-282 14. Fourth-order tensors in a natural basis ---------------------------------------------- The convected components of the spatial stiffness are the lab components of its pull-back. Contraction and inverse act on the components; the identity and the projectors are built from the metric. .. GENERATED FROM PYTHON SOURCE LINES 282-315 .. code-block:: Python L_c = L_iso.to_basis(convected) print("L in the convected basis == L.pull_back(F):", np.allclose(L_c.mat, L_iso.pull_back(F, metric=False).mat)) print("L : eps there, back in the lab == lab result:", np.allclose((L_c @ e).to_basis(None).voigt, (L_iso @ eps).voigt)) print("inverse is the compliance in that basis:", np.allclose(L_c.inverse().to_basis(None).mat, L_iso.inverse().mat)) I4 = sim.Tensor4.identity("stiffness", basis=convected) ginv = convected.inverse_metric print("identity 1/2 (g^ik g^jl + g^il g^jk) raises both indices:", np.allclose((I4 @ e).mat, ginv @ e.mat @ ginv)) print("it is not eye(6) any more:", not np.allclose(I4.mandel, np.eye(6))) print("the mixed (concentration) identity still is:", np.allclose(sim.Tensor4.identity("strain_concentration", basis=convected).mat, np.eye(6))) K, mu = 70000.0 / (3 * (1 - 2 * 0.3)), 70000.0 / (2 * (1 + 0.3)) P_vol = sim.Tensor4.volumetric("stiffness", basis=convected) P_dev = sim.Tensor4.deviatoric("stiffness", basis=convected) print("L = 3K P_vol + 2 mu P_dev with the metric projectors:", np.allclose(L_c.mat, (3 * K * P_vol + 2 * mu * P_dev).mat)) # A concentration tensor (mixed variance) changes basis too A = sim.Tensor4.strain_concentration(sim.A_R(F_gen)) print("concentration tensor, same result through the convected basis:", np.allclose((A.to_basis(convected) @ e).to_basis(None).voigt, (A @ eps).voigt)) try: L_c @ t # stiffness with a stress: not dual except ValueError as err: print("refused:", err) .. rst-class:: sphx-glr-script-out .. code-block:: none L in the convected basis == L.pull_back(F): True L : eps there, back in the lab == lab result: True inverse is the compliance in that basis: True identity 1/2 (g^ik g^jl + g^il g^jk) raises both indices: True it is not eye(6) any more: True the mixed (concentration) identity still is: True L = 3K P_vol + 2 mu P_dev with the metric projectors: True concentration tensor, same result through the convected basis: True refused: In a natural basis a Tensor4 with input indices contravariant contracts a covariant Tensor2, got contravariant (use to_variance()) .. GENERATED FROM PYTHON SOURCE LINES 316-323 15. Push-forward and pull-back carry the basis along ------------------------------------------------------- A push-forward is a transport: another tensor, in the current configuration. In the convected basis it has the components the original had in the reference basis. On a lab tensor ``push_forward`` returns lab components, as always. On a tensor that has its own basis it keeps the components and convects the basis: nothing is computed until lab components are asked for. .. GENERATED FROM PYTHON SOURCE LINES 323-351 .. code-block:: Python reference = sim.Basis(rotation=sim.Rotation.identity(), name="reference") S_ref = S.with_basis(reference) tau_lazy = S_ref.push_forward(F, metric=False) print("components untouched:", np.array_equal(tau_lazy.voigt, S.voigt)) print("basis is now F:", np.allclose(tau_lazy.basis.matrix, F)) print("lab components == eager push-forward:", np.allclose(tau_lazy.to_basis(None).voigt, S.push_forward(F, metric=False).voigt)) # metric=True adds the Piola weight 1/J (Kirchhoff -> Cauchy); use a volume change to see it F_vol = F @ np.diag([1.1, 1.0, 1.0]) cauchy_lazy = S_ref.push_forward(F_vol) # metric=True is the default print("only the 1/J weight touches the numbers:", np.allclose(cauchy_lazy.voigt, S.voigt / np.linalg.det(F_vol))) print("Cauchy, lab == eager:", np.allclose(cauchy_lazy.to_basis(None).voigt, S.push_forward(F_vol).voigt)) back = cauchy_lazy.pull_back(F_vol) print("pull_back returns to the reference basis and components:", np.allclose(back.basis.matrix, np.eye(3)), np.allclose(back.voigt, S.voigt)) # The same holds for a stiffness: the lazy transport costs nothing C_ref = L_iso.with_basis(reference) c_lazy = C_ref.push_forward(F, metric=False) print("Tensor4, lab == eager:", np.allclose(c_lazy.to_basis(None).mat, L_iso.push_forward(F, metric=False).mat)) .. rst-class:: sphx-glr-script-out .. code-block:: none components untouched: True basis is now F: True lab components == eager push-forward: True only the 1/J weight touches the numbers: True Cauchy, lab == eager: True pull_back returns to the reference basis and components: True True Tensor4, lab == eager: True .. GENERATED FROM PYTHON SOURCE LINES 352-356 16. Objectivity ------------------ A superposed rigid rotation turns the convected basis and leaves the components, the metric and every invariant unchanged. .. GENERATED FROM PYTHON SOURCE LINES 356-366 .. code-block:: Python Q = sim.Rotation.from_euler("zxz", [70.0, 25.0, -40.0], degrees=True) tau_Q = tau.rotate(Q) print("components unchanged:", np.array_equal(tau_Q.voigt, tau.voigt)) print("basis is Q F:", np.allclose(tau_Q.basis.matrix, Q.as_matrix() @ F)) print("metric unchanged:", np.allclose(tau_Q.basis.metric, tau.basis.metric)) print("Mises unchanged:", np.isclose(tau_Q.mises(), tau.mises())) print("lab components are those of Q tau Q^T:", np.allclose(tau_Q.to_basis(None).voigt, tau_lab.rotate(Q).voigt)) .. rst-class:: sphx-glr-script-out .. code-block:: none components unchanged: True basis is Q F: True metric unchanged: True Mises unchanged: True lab components are those of Q tau Q^T: True .. GENERATED FROM PYTHON SOURCE LINES 367-371 17. A convected basis per Gauss point ---------------------------------------- With one deformation gradient per point the basis is a batch; it refers to the caller's ``F`` array without copying it. .. GENERATED FROM PYTHON SOURCE LINES 371-386 .. code-block:: Python F_batch = np.eye(3) + 0.15 * rng.standard_normal((N, 3, 3)) conv_batch = sim.Basis.from_F(F_batch) S_batch = sim.Tensor2.stress(rng.standard_normal((N, 6)) * 50.0) tau_batch = S_batch.with_basis(conv_batch) print(tau_batch) print("no arithmetic: components are those of S:", np.array_equal(tau_batch.voigt, S_batch.voigt)) print("Mises per point matches the lab:", np.allclose(tau_batch.mises(), tau_batch.to_basis(None).mises())) print("lab components == batch push-forward:", np.allclose(tau_batch.to_basis(None).voigt, S_batch.push_forward(F_batch, metric=False).voigt)) print("metric of point 3:\n", tau_batch[3].basis.metric) .. rst-class:: sphx-glr-script-out .. code-block:: none Tensor2(N=1000, type='stress', basis=Basis(natural, N=1000)) no arithmetic: components are those of S: True Mises per point matches the lab: True lab components == batch push-forward: True metric of point 3: [[ 0.8223 0.1778 -0.2343] [ 0.1778 1.5633 0.1852] [-0.2343 0.1852 0.9741]] .. GENERATED FROM PYTHON SOURCE LINES 387-392 18. Housekeeping ------------------- ``np.asarray`` returns the components in the tensor's own basis. A basis survives copying and pickling. Tensor types without a variance cannot be written in a natural basis, and a name alone is not a basis. .. GENERATED FROM PYTHON SOURCE LINES 392-421 .. code-block:: Python print("np.asarray(tau):", np.asarray(tau)) clone = pickle.loads(pickle.dumps(sigma_m)) print("pickled:", clone, "| equal:", clone == sigma_m, "| deepcopy equal:", copy.deepcopy(sigma_m) == sigma_m) symmetric = sim.Tensor2.from_voigt(sigma.voigt, "symmetric") print("a symmetric tensor in an orthonormal basis is fine:", symmetric.to_basis(material).basis) for label, operation in [("symmetric in a natural basis", lambda: symmetric.with_basis(convected)), ("a name instead of a basis", lambda: sigma.with_basis("material"))]: try: operation() except (ValueError, TypeError) as err: print(f"{label}: {type(err).__name__}") # A two-point or non-symmetric tensor (F, R, PK1) takes the type "none" -- no # Voigt convention, 9 components: lab-lab components, rotation of both legs, # invariants, but no Voigt vector, no variance and no transport. F_t = sim.Tensor2.from_mat(F_gen, "none") print(F_t, "| stored exactly:", np.array_equal(F_t.mat, F_gen), "| det J =", round(F_t.det(), 6)) print("rotated with both legs:", np.allclose(F_t.rotate(Q).mat, Q.as_matrix() @ F_gen @ Q.as_matrix().T)) for label, operation in [("a symmetric type refuses F", lambda: sim.Tensor2.stress(F_gen)), ("F has no convected basis", lambda: F_t.with_basis(convected)), ("F is not transported", lambda: F_t.push_forward(F_gen))]: try: operation() except ValueError as err: print(f"{label}: ValueError") .. rst-class:: sphx-glr-script-out .. code-block:: none np.asarray(tau): [ 0. 80. 0. 0. 0. 0.] pickled: Tensor2(type='stress', basis=Basis(orthonormal, name='material')) | equal: True | deepcopy equal: True a symmetric tensor in an orthonormal basis is fine: Basis(orthonormal, name='material') symmetric in a natural basis: ValueError a name instead of a basis: TypeError Tensor2(type='none') | stored exactly: True | det J = 1.196 rotated with both legs: True a symmetric type refuses F: ValueError F has no convected basis: ValueError F is not transported: ValueError .. GENERATED FROM PYTHON SOURCE LINES 422-430 19. Bases that come from elsewhere: ply frames, shells, successive transports ------------------------------------------------------------------------------ A ``Basis`` is just the matrix of its vectors, so every case goes through the same three operations. A composite ply knows its stress in the fibre frame; pushing it forward convects that frame (the first vector becomes the stretched fibre ``F a0``). ``F`` and the rotations of ``rotate`` are always given in lab components; a deformation gradient known in ply axes is brought to the lab first, ``F = R F_hat R^T``. .. GENERATED FROM PYTHON SOURCE LINES 430-469 .. code-block:: Python ply = sim.Basis(rotation=sim.Rotation.from_euler("z", 30.0, degrees=True), name="ply") S_ply = sim.Tensor2.stress(np.array([500.0, 20.0, 0.0, 10.0, 0.0, 0.0])).with_basis(ply) # PK2, ply axes F_ply = np.array([[1.1, 0.4, 0.0], [0.0, 0.95, 0.0], [0.0, 0.0, 1.0]]) tau_ply = S_ply.push_forward(F_ply, metric=False) print("convected ply basis = F R:", np.allclose(tau_ply.basis.matrix, F_ply @ ply.matrix)) fibre = tau_ply.basis.matrix[:, 0] print("convected fibre F a0:", fibre, "| stretch:", round(np.linalg.norm(fibre), 4)) print("fibre component kept: tau^11 = S^11 =", tau_ply.voigt[0]) print("lab components:", tau_ply.to_basis(None).voigt) # The orthonormal frame closest to the convected one: its polar rotation (the ply # frame turned by the rotation of F) print("in the polar (rotated orthonormal) ply frame:", tau_ply.to_basis(tau_ply.basis.polar).voigt) # Successive transports compose, and the total pull-back restores the ply frame F1 = np.array([[1.05, 0.1, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) F2 = np.array([[1.0, 0.0, 0.2], [0.0, 1.1, 0.0], [0.0, 0.0, 0.9]]) two_steps = S_ply.push_forward(F1, metric=False).push_forward(F2, metric=False) print("F2 after F1 == F2 F1:", np.allclose(two_steps.basis.matrix, F2 @ F1 @ ply.matrix), "| components kept:", np.array_equal(two_steps.voigt, S_ply.voigt)) print("pull-back by F2 F1 restores the ply frame:", np.allclose(two_steps.pull_back(F2 @ F1, metric=False).basis.matrix, ply.matrix)) # A shell: the covariant surface basis (a_1, a_2, n) is a natural basis; a # membrane stress has contravariant components N^{ab} in it. a1 = np.array([1.0, 0.2, 0.0]) a2 = np.array([-0.1, 0.8, 0.3]) normal = np.cross(a1, a2) normal /= np.linalg.norm(normal) shell = sim.Basis(vectors=np.column_stack([a1, a2, normal]), name="shell") N_shell = sim.Tensor2.stress(np.array([120.0, 40.0, 0.0, 15.0, 0.0, 0.0])).with_basis(shell) print("surface metric a_ab:\n", shell.metric) print("membrane stress: trace", round(N_shell.trace(), 4), "| Mises", round(N_shell.mises(), 4), "| lab Mises", round(N_shell.to_basis(None).mises(), 4)) lamina = sim.Basis(rotation=sim.Rotation.from_matrix(np.linalg.qr(shell.matrix)[0]), name="lamina") print("the same stress in the orthonormal lamina frame:", N_shell.to_basis(lamina).voigt) .. rst-class:: sphx-glr-script-out .. code-block:: none convected ply basis = F R: True convected fibre F a0: [1.1526 0.475 0. ] | stretch: 1.2467 fibre component kept: tau^11 = S^11 = 500.0 lab components: [660.4113 134.1659 0. 278.9151 0. 0. ] in the polar (rotated orthonormal) ply frame: [776.2664 18.3108 0. 58.3352 0. 0. ] F2 after F1 == F2 F1: True | components kept: True pull-back by F2 F1 restores the ply frame: True surface metric a_ab: [[ 1.04 0.06 -0. ] [ 0.06 0.74 0. ] [-0. 0. 1. ]] membrane stress: trace 156.2 | Mises 117.835 | lab Mises 117.835 the same stress in the orthonormal lamina frame: [126.7385 29.4615 0. 15.1479 0. 0. ] .. GENERATED FROM PYTHON SOURCE LINES 470-476 20. The variance tag: raising and lowering indices ----------------------------------------------------- The variance of the components is a tag defaulted from the type (stress contravariant, strain covariant, stiffness (contra, contra), compliance (co, co), concentration tensors mixed). ``to_variance`` gives the same tensor with the other variance: nothing changes in the lab, the metric acts in a natural basis. .. GENERATED FROM PYTHON SOURCE LINES 476-514 .. code-block:: Python print("defaults:", sigma.variance, "|", eps.variance, "|", L_iso.variance, "|", sim.Tensor4.strain_concentration(np.eye(6)).variance) tau_low = tau.to_variance("covariant") # tau_ij = g_ik tau^kl g_lj print(tau_low) print("lowered with the metric:", np.allclose(tau_low.mat, convected.metric @ tau.mat @ convected.metric)) print("same tensor: same Mises", np.isclose(tau_low.mises(), tau.mises()), "| same lab components", np.allclose(tau_low.to_basis(None).voigt, tau.to_basis(None).voigt)) print("raised back:", np.allclose(tau_low.to_variance("contravariant").mat, tau.mat)) # Contractions read the tags: dual variances contract as they are, equal ones # through the metric; the value is the lab one whatever the tags s_conv, e_conv = sigma.to_basis(convected), eps.to_basis(convected) print("sigma : eps, lab:", round(sigma % eps, 6), "| sigma^sharp : eps_flat:", round(s_conv % e_conv, 6), "| sigma_flat : eps^sharp:", round(s_conv.to_variance("covariant") % e_conv.to_variance("contravariant"), 6), "| sigma_flat : eps_flat (metric):", round(s_conv.to_variance("covariant") % e_conv, 6)) # A stiffness with lowered indices contracts a contravariant strain L_low = L_iso.to_basis(convected).to_variance(("covariant", "covariant")) print(L_low) print("L_low : e^sharp back in the lab == L : eps:", np.allclose((L_low @ e_conv.to_variance("contravariant")).to_basis(None).voigt, (L_iso @ eps).voigt)) try: L_low @ e_conv except ValueError as err: print("not dual:", str(err).split(" (")[0]) # The variance is part of what a transport means: lowering the indices of a # lab tensor and pushing it forward is not the push-forward of the original. pushed = sigma.push_forward(F_gen, metric=False) pushed_low = sigma.to_variance("covariant").push_forward(F_gen, metric=False) print("F sigma F^T vs F^-T sigma F^-1 differ:", not np.allclose(pushed.mat, pushed_low.mat)) # A "symmetric" tensor has no default variance: declare it first sym = sim.Tensor2.from_voigt(sigma.voigt, "symmetric") print("symmetric, undeclared:", sym.variance, "| declared:", sym.to_variance("contravariant").variance) .. rst-class:: sphx-glr-script-out .. code-block:: none defaults: contravariant | covariant | ('contravariant', 'contravariant') | ('covariant', 'contravariant') Tensor2(type='stress', basis=Basis(natural, name='convected'), variance='covariant') lowered with the metric: True same tensor: same Mises True | same lab components True raised back: True sigma : eps, lab: 1.0 | sigma^sharp : eps_flat: 1.0 | sigma_flat : eps^sharp: 1.0 | sigma_flat : eps_flat (metric): 1.0 Tensor4(type='stiffness', basis=Basis(natural, name='convected'), variance=('covariant', 'covariant')) L_low : e^sharp back in the lab == L : eps: True not dual: In a natural basis a Tensor4 with input indices covariant contracts a contravariant Tensor2, got covariant F sigma F^T vs F^-T sigma F^-1 differ: True symmetric, undeclared: None | declared: contravariant .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.025 seconds) .. _sphx_glr_download_examples_continuum_mechanics_tensor_basis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tensor_basis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tensor_basis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tensor_basis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_