Evaluation of function with options affecting complexity

Question

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?

Answer 1

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 *)

Answer 2

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
]
...