Stress gradient: Linear elastic medium, 3D

Define P(F), (Subscript[∇, F])P

SetDirectory[NotebookDirectory[]]

Import 3d flux definitions

NotebookEvaluate[Directory[] <> "/../../FluxDefinitions3d.nb"];

Strain energy

Constitutive relation

c[i_,j_,k_,l_] := λ δ[[i,j]] δ[[k,l]] + μ(δ[[i,k]] δ[[j,l]]+δ[[i,l]] δ[[j,k]])

Engineering strain in terms of the deformation gradient

ϵ=Table[(Subscript[F, idx2[i,j]])-δ[[i,j]],{i,d},{j,d}]

Strain energy

ℰ =(1/2)Simplify[ Sum[c[i,j,k,l] ϵ[[i,j]] ϵ[[k,l]],{i,d},{j,d},{k,d},{l,d}] /. LinearizeFTensor]

Hessian of strain energy

Hℰ=Table[Simplify[D[D[ℰ,(Subscript[F, i])],(Subscript[F, j])]],{i,(d^2)},{j,(d^2)}];

TableForm[Hℰ]

Prules = Flatten[Table[(Subscript[P, i,j])->Hℰ[[i,j]],{i,(d^2)},{j,(d^2)}]]

Average Hessian at interface between left/right states

Define integration path between left and right states

FpathRules=Table[(Subscript[F, i])->(Subscript[F, L,i])+((Subscript[F, R,i])-(Subscript[F, L,i]))ξ,{i,(d^2)}]

Also define a parametrization for variation of material properties

MaterialPathRules={λ->lambdaL+(lambdaR-lambdaL)ξ,μ->muL+(muR-muL)ξ}

Hℰpath[ξ_] =Table[Expand[ Hℰ[[i,j]] /. FpathRules  /. MaterialPathRules],{i,(d^2)},{j,(d^2)}];

Manipulate[Hℰpath[ξ][[i,j]],{i,1,(d^2),1},{j,1,(d^2),1}]

If[FileExistsQ["Hav.mx"],Print["Reading Hav.mx ..."]; <<Hav.mx; Print["Done."],Hav=Monitor[ Table[Simplify[Integrate[Hℰpath[ξ][[i,j]],{ξ,0,1}] ],{i,(d^2)},{j,(d^2)}],{i,j}];DumpSave["Hav.mx",Hav]];

Manipulate[Hav[[i,j]],{i,1,(d^2),1},{j,1,(d^2),1}]

M=Table[H[[k]].U[[k]],{k,d}];
Map[TableForm[#]&,M]

Prules = Flatten[Table[(Subscript[P, i,j])->Hav[[i,j]],{i,(d^2)},{j,(d^2)}]]

MP = Monitor[Table[M[[k,i,j]] /. Prules,{k,1,d,1},{i,1,d,1},{j,1,d,1}],{k,i,j}];

Manipulate[MP[[k,i,j]] ,{k,1,d,1},{i,1,d,1},{j,1,d,1}]

TableForm[MP]

Generate Fortran code

Modify StressGradient Fortran template to take into account specific material properties

dPdF[]:=Module[{i,j,k,n,coef,instr,subr},subr="SUBROUTINE StressGradient(k,MatParams,F,M)\n";
subr=subr<>"  INTEGER,PARAMETER::MaxParams=200,L=1,R=2\n";
subr=subr<>"  INTEGER k\n";
subr=subr<>"  REAL (KIND=8) MatParams(MaxParams,2),F(9,2),M(3,3)\n";
subr=subr<>"  INTEGER i,j\n";
subr=subr<>"  REAL (KIND=8) lambdaL,lambdaR,muL,muR\n";
subr=subr<>"  lambdaL=MatParams(1,1); muL=MatParams(2,1)\n";
subr=subr<>"  lambdaR=MatParams(1,2); muR=MatParams(2,2)\n";
subr=subr<>"  SELECT CASE(k)\n";
For[k=1,k≤d,k++,
subr=subr<>"    CASE(" <> ToString[k] <>")\n";
For[i=1,i≤d,i++,For[j=i,j≤d,j++,coef=ToString[FortranForm[M[[k,i,j]]/.Prules]];
instr="      M("<>ToString[i]<>","<>ToString[j]<>") = "<>coef<>"\n";
subr=subr<>instr;];
For[j=1,j<i,j++,
subr=subr<> "      M("<>ToString[i]<>","<>ToString[j]<>") = M(" <>ToString[j]<>","<>ToString[i]<>")\n"];
];
];
subr=subr<>" END SELECT\n";
subr=subr<>"END SUBROUTINE StressGradient\n";
Export["StressGradient.f90",subr,"Text"];
Return[subr]]

dPdF[];