3
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ What about Switch[OptionValue[symmetry] , "simple", generic , "more symmetric" , A1[generic] , "much more symmetric" , A2 @ A1 @ generic , _, $Failed ] ? $\endgroup$ – Kuba Mar 13 '18 at 13:50
  • $\begingroup$ Hi, thanks for the reply. Maybe I stripped too much out of my actual function for this to make sense. A1, A2 are tensor products (of reflection symmetries) which get somehow contracted with the generic tensor. Writing each one down requires a few lines of code. Applying A1 strips the tensor of some components. Then applying A2 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 down A1 for the first case, then A1 and A2 for the second etc. Does that make sense? $\endgroup$ – gpap Mar 13 '18 at 17:38
  • 1
    $\begingroup$ I'm not sure I understand the problem, can't you just have A1 and A2 defined somewhere else and refer to them instead of putting them explicitly in the code (that is the problem, right?). They do not depend on head, do they? $\endgroup$ – Kuba Mar 13 '18 at 19:42
  • $\begingroup$ Thanks for the interest and sorry I took long to edit - it's been a busy few days. I have edited with the exact code I am interested in improving, hope it doesn't confuse things. $\endgroup$ – gpap Mar 16 '18 at 11:45
3
$\begingroup$

You could precalculate all tensors with a dummy head, then, when your function is called, simply replace this dummy head with appropriate one.

BeginPackage@"ElasticityTensor`"; Unprotect@"`*"; ClearAll@"`*"

CreateElasticityTensor::usage = "CreateElasticityTensor[head] ...";
CreateElasticityTensor::unknownSymmetry = "`1` is not one of known symmetry specifications: `2`.";

Begin@"`Private`"; ClearAll@"`*"

reflection@v_ := With[{r = ReflectionMatrix@v}, r\[TensorProduct]r\[TensorProduct]r\[TensorProduct]r]
rotation[θ_, v_] := With[{r = RotationMatrix[θ, v]}, r\[TensorProduct]r\[TensorProduct]r\[TensorProduct]r]

symmetrize[tens_, symm_] := tens /. First@Quiet[
  Solve[
    tens == TensorContract[symm\[TensorProduct]tens, {{2, 9}, {4, 10}, {6, 11}, {8, 12}}],
    DeleteDuplicates@Cases[Flatten@tens, _dummy]
  ],
  Solve::svars
]

symmetric@"Generic" = Array[dummy, {3, 3, 3, 3}] /.
  dummy[a_, b_, d_, e_] /; d > e :> dummy[a, b, e, d] /.
  dummy[a_, b_, d_, e_] /; a > b :> dummy[b, a, d, e] /.
  dummy[a_, b_, d_, e_] /; ((a > d && b >= e) || (a >= d && b > e) || (a > d && b < e)) :> dummy[d, e, a, b];

symmetric@"Monoclinic" =
  symmetrize[symmetric@"Generic", reflection@{0, 0, 1}];

symmetric@"Orthotropic" =
  symmetrize[symmetric@"Monoclinic", reflection@{1, 0, 0}];

symmetric@"Transverse Isotropic" =
  symmetrize[symmetric@"Orthotropic", rotation[\[Pi]/4, {0, 0, 1}]];

symmetric@"Isotropic" =
  symmetrize[symmetric@"Transverse Isotropic", rotation[\[Pi]/4, {1, 0, 0}]];

symmetric@sym_ := (
  Message[
    CreateElasticityTensor::unknownSymmetry,
    sym,
    Cases[DownValues@symmetric, (_@_@s_String :> _) :> s]
  ];
  $Failed
)

CreateElasticityTensor // Options = {"Symmetry" -> "Generic"};
CreateElasticityTensor[head_, OptionsPattern[]] :=
  symmetric@OptionValue@"Symmetry" /. dummy -> head

End[]; Protect@"`*"; EndPackage[];

Results are the same, but function is much faster.

res1 = createElasticityTensor[h]; // MaxMemoryUsed // RepeatedTiming
res2 = CreateElasticityTensor[h]; // MaxMemoryUsed // RepeatedTiming
res1 === res2
(* {0.00070, 23592} *)
(* {0.000093, 9152} *)
(* True *)

res1 = createElasticityTensor[h, "Symmetry" -> "Isotropic"]; // MaxMemoryUsed // AbsoluteTiming
res2 = CreateElasticityTensor[h, "Symmetry" -> "Isotropic"]; // MaxMemoryUsed // AbsoluteTiming
res1 === res2
(* {1.40616, 64547248} *)
(* {0.000121, 7344} *)
(* True *)

res1 = createElasticityTensor[h, "Symmetry" -> "foo"]; // MaxMemoryUsed // AbsoluteTiming
res2 = CreateElasticityTensor[h, "Symmetry" -> "foo"]; // MaxMemoryUsed // AbsoluteTiming
res1 === res2
(* {1.41578, 64548288} *)
(* CreateElasticityTensor::unknownSymmetry: foo is not one of known symmetry specifications: {Generic, Isotropic, Monoclinic, Orthotropic, Transverse Isotropic}. *)
(* {0.001234, 23416} *)
(* True *)
$\endgroup$
  • $\begingroup$ Ach, this makes much more sense - thanks so much! $\endgroup$ – gpap Mar 16 '18 at 20:30
2
$\begingroup$

Too long for a comment, what about:

...
add[action_]:= generic = First @ solution[
  generic,contraction[action\[TensorProduct]generic]
];

Switch[OptionValue[Symmetry],
  Except @ Alternatives["Generic","Monoclinic","Orthotropic","Transverse Isotropic","Isotropic"], 
  $Failed,
  "Generic", 
  generic,
  add[reflection[{0,0,1}]];  "Monoclinic", 
  generic,
  add[reflection[{1,0,0}]];  "Orthotropic", 
  generic,
  add[rotation[π/4,{0,0,1}]];"Transverse Isotropic", 
  generic,
  add[rotation[π/4,{1,0,0}]];"Isotropic",
  generic
]
...
$\endgroup$
  • $\begingroup$ This is not really that much different to what I have, other than subbing Which for Switch though :\ It is neater so +1 $\endgroup$ – gpap Mar 16 '18 at 20:32
  • $\begingroup$ @gpap yes it is rather a take on readablity and the flow than on math etc. p.s. Except at the beggining makes a difference. $\endgroup$ – Kuba Mar 16 '18 at 20:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.