Add all code

Author: singularitti

Elastic Formulae/cubic.wls

@@ -0,0 +1,25 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+cubic = {{c11, c12, c12, 0, 0, 0}, {0, c11, c12, 0, 0, 0}, {0, 0, c11,
+        0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c44, 0}, {0, 0, 0, 0, 0,
+        c44}} // Symmetrise
+
+f = ElasticEnergy[cubic]
+
+cj = {c11, c12, c44};
+
+(*Length[cj] == 3*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/hexagonal.wls

@@ -0,0 +1,25 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+hex = {{c11, c12, c13, 0, 0, 0}, {0, c11, c13, 0, 0, 0}, {0, 0, c33, 
+      0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c44, 0}, {0, 0, 0, 0, 0,
+       (c11 - c12) / 2}} // Symmetrise
+
+f = ElasticEnergy[hex]
+
+cj = {c11, c33, c12, c13, c44};
+
+(*Length[cj] == 5*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/monoclinic.wls

@@ -0,0 +1,26 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+mono = {{c11, c12, c13, 0, 0, c16}, {0, c22, c23, 0, 0, c26}, {0, 0, 
+    c33, 0, 0, c36}, {0, 0, 0, c44, c45, 0}, {0, 0, 0, 0, c55, 0}, {0, 0,
+     0, 0, 0, c66}} // Symmetrise
+
+f = ElasticEnergy[mono]
+
+cj = {c11, c22, c33, c12, c13, c23, c44, c55, c66, c16, c26, c36, c45
+    };
+
+(*Length[cj] == 13*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/monoclinic2.wls

@@ -0,0 +1,26 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+mono = {{c11, c12, c13, 0, c15, 0}, {0, c22, c23, 0, c25, 0}, {0, 0, 
+    c33, 0, c35, 0}, {0, 0, 0, c44, 0, c46}, {0, 0, 0, 0, c55, 0}, {0, 0,
+     0, 0, 0, c66}} // Symmetrise
+
+f = ElasticEnergy[mono]
+
+cj = {c11, c22, c33, c12, c13, c23, c44, c55, c66, c15, c25, c35, c46
+    };
+
+(*Length[cj] == 13*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/orthorhombic.wls

@@ -0,0 +1,25 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+orth = {{c11, c12, c13, 0, 0, 0}, {0, c22, c23, 0, 0, 0}, {0, 0, c33,
+     0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c55, 0}, {0, 0, 0, 0, 0,
+     c66}} // Symmetrise
+
+cj = {c11, c22, c33, c12, c13, c23, c44, c55, c66};
+
+f = ElasticEnergy[orth]
+
+(*Length[cj] == 9*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/sij.wl

@@ -0,0 +1,37 @@
+(* ::Package:: *)
+
+BeginPackage["Sij`"]
+
+Symmetrise::usage = "Construct a symmetric matrix from an upper triangular matrix."
+
+ElasticEnergy::usage = "Calculate the free energy formula."
+
+Sij::usage = "Calculate matrix element S[i, j]."
+
+S::usage = "Calculate matrix S."
+
+Examine::usage = "Examine whether the S matrix is correct."
+
+Begin["`Private`"]
+
+\[Epsilon] = {Global`\[Epsilon]1, Global`\[Epsilon]2, Global`\[Epsilon]3,
+     Global`\[Epsilon]4, Global`\[Epsilon]5, Global`\[Epsilon]6};
+
+Symmetrise[upperTriMatrix_] :=
+    UpperTriangularize[upperTriMatrix, 1] + Transpose[upperTriMatrix]
+
+ElasticEnergy[c_] :=
+    \[Epsilon] . c . \[Epsilon] / 2 // Expand
+
+Sij[f_, \[Epsilon]i_, cj_] :=
+    D[f, \[Epsilon]i, cj]
+
+S[f_, c_] :=
+    Table[Sij[f, \[Epsilon][[i]], c[[j]]], {i, 1, 6}, {j, 1, Length[c
+        ]}] // Simplify
+
+Examine[S_, c_, f_] := PossibleZeroQ[S . c . \[Epsilon] / 2 - f]
+
+End[]
+
+EndPackage[]

Elastic Formulae/tetragonal.wls

@@ -0,0 +1,25 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+tet = {{c11, c12, c13, 0, 0, 0}, {0, c11, c13, 0, 0, 0}, {0, 0, c33, 
+        0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c44, 0}, {0, 0, 0, 0, 0,
+         c66}} // Symmetrise
+
+cj = {c11, c33, c12, c13, c44, c66};
+
+f = ElasticEnergy[tet]
+
+(*Length[cj] == 6*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/tetragonal2.wls

@@ -0,0 +1,25 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+tet = {{c11, c12, c13, 0, 0, c16}, {0, c11, c13, 0, 0, -c16}, {0, 0, 
+        c33, 0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c44, 0}, {0, 0, 0, 
+        0, 0, c66}} // Symmetrise
+
+cj = {c11, c33, c12, c13, c16, c44, c66};
+
+(*Length[cj] == 7*)
+
+f = ElasticEnergy[tet]
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/triclinic.wls

@@ -0,0 +1,26 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+triclinic = {{c11, c12, c13, c14, c15, c16}, {0, c22, c23, c24, c25, 
+    c26}, {0, 0, c33, c34, c35, c36}, {0, 0, 0, c44, c45, c46}, {0, 0, 0,
+     0, c55, c56}, {0, 0, 0, 0, 0, c66}} // Symmetrise
+
+f = ElasticEnergy[triclinic]
+
+cj = {c11, c12, c13, c14, c15, c16, c22, c23, c24, c25, c26, c33, c34,
+     c35, c36, c44, c45, c46, c55, c56, c66};
+
+(*Length[cj] == 21*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/trigonal.wls

@@ -0,0 +1,27 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+tri = {{c11, c12, c13, c14, 0, 0}, {0, c11, c13, -c14, 0, 0}, {0, 0, 
+        c33, 0, 0, 0}, {0, 0, 0, c44, 0, 0}, {0, 0, 0, 0, c44, c14}, {0, 0, 0,
+         0, 0, (c11 - c12) / 2}} // Symmetrise
+
+(*tri // SymmetricMatrixQ*)
+
+cj = {c11, c33, c12, c13, c44, c14};
+
+f = ElasticEnergy[tri]
+
+(*Length[cj] == 6*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Elastic Formulae/trigonal2.wls

@@ -0,0 +1,27 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+SetDirectory[NotebookDirectory[]];
+
+<<Sij`
+
+tri = {{c11, c12, c13, c14, c15, 0}, {0, c11, c13, -c14, -c15, 0}, {0,
+         0, c33, 0, 0, 0}, {0, 0, 0, c44, 0, -c15}, {0, 0, 0, 0, c44, c14}, {
+        0, 0, 0, 0, 0, (c11 - c12) / 2}} // Symmetrise
+
+(*tri // SymmetricMatrixQ*)
+
+cj = {c11, c33, c12, c13, c44, c14, c15};
+
+f = ElasticEnergy[tri]
+
+(*Length[cj] == 7*)
+
+s = S[f, cj];
+
+s // MatrixForm
+
+
+Examine[s, cj, f]

Sound Velocity/Cubic.wls

@@ -0,0 +1,33 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+GreenChristoffelEquations[{ qx_, qy_, qz_ }] := {{
+        c11 qx ^ 2 + c44(qy ^ 2 + qz ^ 2), (c12 + c44) qx qy, (c12 + c44) qx qz
+    }, {
+        (c12 + c44) qx qy, c11 qy ^ 2 + c44(qz ^ 2 + qx ^ 2), (c12 + c44) qy qz
+    }, {
+        (c12 + c44) qx qz, (c12 + c44) qy qz, c11 qz ^ 2 + c44(qx ^ 2 + qy ^ 2)
+}}
+
+qVectorOnDirection[{ a_, b_, c_ }] := q / Norm[{ a, b, c }] { a, b, c }
+
+generateEquationsOnDirection[{ a_, b_, c_ }] := GreenChristoffelEquations[qVectorOnDirection[{ a, b, c }]]
+
+phaseVelocity[{ a_, b_, c_ }] := Simplify[Sqrt[Eigenvalues[generateEquationsOnDirection[{ a, b, c }]] / \[Rho]] / q, Assumptions-> q > 0]
+
+particleDisplacementDirections[{ a_, b_, c_ }] := Eigenvectors[generateEquationsOnDirection[{ a, b, c }]]
+
+takePositiveSolutions[{ a_, b_, c_ }] := Module[{ v, indices },
+    v = phaseVelocity[{ a, b, c }];
+    indices = Simplify[v // Sign, Assumptions -> c11 > c12 && c11 + 2 c12 > 0 && c44 > 0 && \[Rho] > 0 && q > 0] // Position[#, 1] & // Flatten;
+    Map[Part[v, #] &, indices]
+]
+
+
+takePositiveSolutions[{1,1,1}]
+particleDisplacementDirections[{1,1,1}]
+
+
+Simplify[((generateEquationsOnDirection[{1,1,1}]//Eigenvalues) / \[Rho] // Sqrt) /q,Assumptions->q>0]

Sound Velocity/Orthorhombic.wls

@@ -0,0 +1,39 @@
+#!/usr/bin/env wolframscript
+(* ::Package:: *)
+
+Clear["Global`*"]
+
+indices := {{1,1,1,1},{1,1,2,2},{1,1,3,3},{2,2,1,1},{2,2,2,2},{2,2,3,3},{3,3,1,1},{3,3,2,2},{3,3,3,3},
+{2,3,2,3},{2,3,3,2},{3,2,3,2},{3,2,2,3},{1,3,1,3},{1,3,3,1},{3,1,3,1},{3,1,1,3},
+{1,2,1,2},{1,2,2,1},{2,1,2,1},{2,1,1,2}};
+
+dict = <|{1,1}->1,{2,2}->2,{3,3}->3,{2,3}->4,{3,2}->4,{1,3}->5,{3,1}->5,{1,2}->6,{2,1}->6|>;
+
+c[i_,j_,k_,l_] := c[dict[[Key[{i,j}]]],dict[[Key[{k,l}]]]]
+
+c[i_,j_] := {{c11,c12,c13,0,0,0},{c12,c22,c23,0,0,0},{c13,c23,c33,0,0,0},{0,0,0,c44,0,0},{0,0,0,0,c55,0},{0,0,0,0,0,c66}}[[i,j]];
+
+qs = {qx,qy,qz};
+
+GreenChristoffelEquations[qs_] := Table[Sum[c[\[Alpha],\[Gamma],\[Beta],\[Delta]]qs[[\[Gamma]]]qs[[\[Delta]]],{\[Gamma],1,3},{\[Delta],1,3}],{\[Alpha],1,3},{\[Beta],1,3}]
+
+
+qVectorOnDirection[{ a_, b_, c_ }] := q / Norm[{ a, b, c }] { a, b, c }
+
+generateEquationsOnDirection[{ a_, b_, c_ }] := GreenChristoffelEquations[qVectorOnDirection[{ a, b, c }]]
+
+phaseVelocity[{ a_, b_, c_ }] := Simplify[Sqrt[Eigenvalues[generateEquationsOnDirection[{ a, b, c }]] / \[Rho]] / q, Assumptions-> q > 0]
+
+particleDisplacementDirections[{ a_, b_, c_ }] := Eigenvectors[generateEquationsOnDirection[{ a, b, c }]]
+
+takePositiveSolutions[{ a_, b_, c_ }] := Module[{ v, indices },
+    v = phaseVelocity[{ a, b, c }];
+    indices = Simplify[v // Sign, Assumptions -> c11 > 0 && c11 c22 > c12^2 && c11 c22 c33 + 2 c12 c13 c23 > c11 c23^2 + c22 c13^2 + c33 c12^2 && c44 > 0 && c55 > 0 && c66 > 0 && \[Rho] > 0 && q > 0] // Position[#, 1] & // Flatten;
+    Map[Part[v, #] &, indices]
+]
+
+
+phaseVelocity[{0,1,1}]
+
+
+Simplify[((generateEquationsOnDirection[{1,1,1}]//Eigenvalues) / \[Rho] // Sqrt) /q,Assumptions->q>0]