I am writing a function that produces a tensor with certain symmetries and I use a function option indicating the symmetry. My function definition has the form:
ClearAll@f Options[f] = {symmetry -> "simple"};
f[head_, OptionsPattern[]] :=
Module[{dummy, generic, symmetry1, symmetry2},
generic = Array[head, {4, 4, 4}];
Which[dummy = generic; OptionValue[symmetry] == "simple", dummy,
symmetry1 = A1; dummy = symmetry1[dummy];
OptionValue[symmetry] == "more symmetric", dummy, symmetry2 = A2;
dummy = symmetry2[dummy];
OptionValue[symmetry] == "much more symmetric", dummy,
True, $Failed]]
The symmetry1
and symmetry2
are relatively complicated tensor contractions so instead of generating them in the main body of Module
, I generate them after Which
tries to match the value of the option symmetry
. There are a few things that make me a bit uneasy about my function (like redefining the dummy
tensor in terms of itself or relying on Which
to search through the cases incrementally) but it seems to work faster than generating symmetry1
and symmetry2
irrespective of whether they are needed.
The actual code is a symbolic generator for an elasticity tensor
ClearAll[createElasticityTensor];
Options[createElasticityTensor] = {Symmetry -> "Generic"};
createElasticityTensor[head_: c, OptionsPattern[]] :=
Module[{generic, dummy, contraction, solution, reflection, rotation},
contraction =
TensorContract[#, {{2, 9}, {4, 10}, {6, 11}, {8, 12}}] &;
solution =
Quiet@Solve[#1 == #2, Cases[DeleteDuplicates@Flatten@#1, _head]] &;
reflection[v_] := ReflectionMatrix[v]\[TensorProduct]
ReflectionMatrix[v]\[TensorProduct]
ReflectionMatrix[v]\[TensorProduct]
ReflectionMatrix[v];
rotation[θ_, v_] :=
RotationMatrix[θ, v]\[TensorProduct]
RotationMatrix[θ, v]\[TensorProduct]
RotationMatrix[θ, v]\[TensorProduct]
RotationMatrix[θ, v];
generic =
Array[head, {3, 3, 3, 3}] /. (head[a_, b_, d_, e_] /; d > e :>
head[a, b, e, d]) /. (head[a_, b_, d_, e_] /; a > b :>
head[b, a, d, e]) /. (head[a_, b_, d_,
e_] /; ((a > d && b >= e) || (a >= d && b > e) || (a > d &&
b < e)) :> head[d, e, a, b]);
Which[
dummy = generic;
OptionValue[Symmetry] == "Generic", dummy,
dummy =
dummy /.
First@solution[dummy,
contraction[reflection[{0, 0, 1}]\[TensorProduct]dummy]];
OptionValue[Symmetry] == "Monoclinic", dummy,
dummy =
dummy /.
First@solution[dummy,
contraction[reflection[{1, 0, 0}]\[TensorProduct]dummy]];
OptionValue[Symmetry] == "Orthotropic", dummy,
dummy =
dummy /.
First@solution[dummy,
contraction[
rotation[π/4, {0, 0, 1}]\[TensorProduct]dummy]];
OptionValue[Symmetry] == "Transverse Isotropic", dummy,
dummy =
dummy /.
First@solution[dummy,
contraction[
rotation[π/4, {1, 0, 0}]\[TensorProduct]dummy]];
OptionValue[Symmetry] == "Isotropic", dummy,
True, $Failed
]
];
Each new option has new increased symmetry. It is evaluated and contracted with the previously constructed tensor for which symbolic components are solved for based on the symmetry. That process of solving for the symmetry is an expensive operation so I would like to avoid repeating it. However, as this has part educational purposes, I would like to retain as much mathematical transparency as I can (i.e., I don't want to replace a tensor contraction with a Table
or Map
line).
My question is: is this the best way to go through the symmetries from the options and, if not, what would be a better way to write my function so that the symmetries are only evaluated when needed?
Switch[OptionValue[symmetry] , "simple", generic , "more symmetric" , A1[generic] , "much more symmetric" , A2 @ A1 @ generic , _, $Failed ]
? $\endgroup$ – Kuba♦ Mar 13 '18 at 13:50A1
,A2
are tensor products (of reflection symmetries) which get somehow contracted with thegeneric
tensor. Writing each one down requires a few lines of code. ApplyingA1
strips the tensor of some components. Then applyingA2
to the result, strips it of more components and so on. I don't want to calculate either unless I have to and it would be really cumbersome (using your method) to write downA1
for the first case, thenA1
andA2
for the second etc. Does that make sense? $\endgroup$ – gpap Mar 13 '18 at 17:38A1
andA2
defined somewhere else and refer to them instead of putting them explicitly in the code (that is the problem, right?). They do not depend onhead
, do they? $\endgroup$ – Kuba♦ Mar 13 '18 at 19:42