Print["This is a test"];
Unprotect[L,𝒟,δ,(Subscript[α, 0]),n,Δα, α, Ε,ν,Ω];
Clear["Global`*"]
ClearAll["Global`*"]
Clear["Global`*"]
Remove["Global`*"]
<<Notation`
Off[Symbolize::bsymbexs]
Symbolize[(NotationTemplateTag[𝒞])]; 
Symbolize[(NotationTemplateTag[(OverVector[X])])]; Symbolize[(NotationTemplateTag[(Subscript[X, 1])])]; Symbolize[(NotationTemplateTag[(Subscript[X, 2])])]; Symbolize[(NotationTemplateTag[(Subscript[X, 3])])] 
Symbolize[(NotationTemplateTag[(OverVector[U])])];Symbolize[(NotationTemplateTag[(Subscript[U, 1])])]; Symbolize[(NotationTemplateTag[(Subscript[U, 2])])]; Symbolize[(NotationTemplateTag[(Subscript[U, 3])])] 
Symbolize[(NotationTemplateTag[(Subscript[Un, 1])])]; Symbolize[(NotationTemplateTag[(Subscript[Un, 2])])]; Symbolize[(NotationTemplateTag[(Subscript[Un, 3])])] 
Symbolize[(NotationTemplateTag[(Subscript[Γ, (Subscript[g, _])])])];
Symbolize[(NotationTemplateTag[(Subscript[w, 0])])];
Symbolize[(NotationTemplateTag[(Subscript[a, 0])])];
Symbolize[(NotationTemplateTag[(Subscript[α, 0])])];
Symbolize[(NotationTemplateTag[δa])];
Symbolize[(NotationTemplateTag[(Subscript[n, 1])])]; Symbolize[(NotationTemplateTag[(Subscript[n, 2])])];Symbolize[(NotationTemplateTag[(Subscript[n, 3])])] 
Notation[(NotationTemplateTag[(Subscript[x_, i_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·k_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_,k_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·k_·l_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_,k_,l_]]])]


Notation[(NotationTemplateTag[(((x_))^(-1))  ]) ⟺ (NotationTemplateTag[ Inverse[x_]])]
Notation[(NotationTemplateTag[A_⊙B_  ]) ⟺ (NotationTemplateTag[ Sum[A_[[i,j]]B_[[i,j]],{i,1,3},{j,1,3}]])]
Notation[(NotationTemplateTag[A_\[Transpose]  ]) ⟺ (NotationTemplateTag[ Transpose[A_]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·k_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_,k_]]])]
Notation[(NotationTemplateTag[(Subscript[x_, i_ ·j_·k_·l_·])  ]) ⟺ (NotationTemplateTag[ x_[[i_,j_,k_,l_]]])]

(*<<HaneeshPackages`*)
PrntRls={(Subscript[U, 1])-> Subscript[U,1],(Subscript[U, 2])-> Subscript[U,2],(Subscript[U, 3])-> Subscript[U,3],(Subscript[X, 1])-> Subscript[X,1],(Subscript[X, 2])-> Subscript[X,2],(Subscript[X, 3])-> Subscript[X,3]};

Needs["NDSolve`FEM`"]
dd=DateValue[{"Hour12","Minute","AMPM"}];
Text[ToString@dd[[1]]<>":"<>ToString@(dd[[2]])<>" "<>ToString@dd[[3]]]

L=Quantity[160,"micrometers"];L::usage="Total length of the beam";
𝒟=Quantity[50,"micrometers"];𝒟::usage="Diameter of the spicule";
Ε=Quantity[25,"GPa"]; Ε::usage="Young's modulus";
(Subscript[a, 0])=
δ=𝒟/L;δ::usage="Non-dimensional spicule diameter";
(Subscript[w, 0])= 𝒟/100; (Subscript[w, 0])::usage="Indentation depth";
ν=0.3;ν::usage="Poisson's ratio";
ϵ=0.1δ;ϵ::usage="a small number";
(Subscript[α, 0])::usage="Notch length";
n=
Δα= (δ/2) (1/n);
(Subscript[α, 0])=δ/10;(*((Subscript[a, 0])/L)*)
Δα::usage="steps size of virtual crack increment for computing numerical derivative of energy"
Ω=Cylinder[{{0,0,0},{1/2,0,0}},δ/2];
Protect[L,𝒟,δ,(Subscript[α, 0]),Δα,  Ε,ν,Ω];

mesh=ToElementMesh[Ω,MaxCellMeasure ->{"Length"->.03}]
Basis={Graphics3D[Arrow[{{0,0,0},{2δ,0,0}}]],Graphics3D[Text["(Subscript[(UnderBar[(OverHat[e])]), 1])",{2δ+ϵ,0,0}]],
Graphics3D[Arrow[{{0,0,0},{0,2δ,0}}]],Graphics3D[Text["(Subscript[(UnderBar[(OverHat[e])]), 2])",{0,2δ+ϵ,0}]],
Graphics3D[Arrow[{{0,0,0},{0,0,2δ}}]],Graphics3D[Text["(Subscript[(UnderBar[(OverHat[e])]), 3])",{0,0,2δ+ϵ}]],
Graphics3D[{Red,Arrow[{{1/2,δ,0},{1/2,δ/2,0}}]}],Graphics3D[Text["(Subscript[(UnderBar[(OverHat[e])]), 3])",{0,0,2δ+ϵ}]],
Graphics3D[{Green,Cuboid[{0,δ,δ},{-δ,-δ,-δ}]}]
};

(*  *)
Show[mesh["Wireframe"],Basis, 
ViewPoint->{0,10,10}]

vd=NDSolve`VariableData[{"DependentVariables","Space"}->{{(Subscript[U, 1]),(Subscript[U, 2]),(Subscript[U, 3])},{(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])}}];
sd=NDSolve`SolutionData["Space"->ToNumericalRegion[mesh]];
methodData=InitializePDEMethodData[vd,sd];

c11 ={{λ+2μ,0,0},{0,μ,0},{0,0,μ}};
c12={{0,λ,0},{μ,0,0},{0,0,0}};
c13={{0,0,λ},{0,0,0},{μ,0,0}};

c22 ={{μ,0,0},{0,λ+2μ,0},{0,0,μ}};
c21={{0,μ,0},{λ,0,0},{0,0,0}};
c23={{0,0,0},{0,0,λ},{0,μ,0}};


c33 ={{μ,0,0},{0,μ,0},{0,0,λ+2μ}};
c31={{0,0,μ},{0,0,0},{λ,0,0}};
c32={{0,0,0},{0,0,μ},{0,λ,0}};

Block[{λ=(ν/((1+ν)(1-2ν))),μ=(1/(2(1+ν)))},pdeCoefficients=InitializePDECoefficients[vd,sd,"DiffusionCoefficients"->{{c11,c12,c13},{c21,c22,c23},{c31,c32,c33}}]];

discretePDE=DiscretizePDE[pdeCoefficients,methodData,sd];

(Subscript[Γ, (Subscript[g, 1])])=DirichletCondition[(Subscript[U, 1])[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]==0,(Subscript[X, 1])==0];
(Subscript[Γ, (Subscript[g, 2])])=DirichletCondition[(Subscript[U, 2])[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]==0,(Subscript[X, 1])==0];
(Subscript[Γ, (Subscript[g, 3])])=DirichletCondition[(Subscript[U, 3])[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]==0,(Subscript[X, 1])==0];
(Subscript[Γ, (Subscript[g, 5])])=DirichletCondition[(Subscript[U, 2])[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]==-((Subscript[w, 0])/L),  ((Subscript[X, 2])^2)+((Subscript[X, 3])^2)==(((δ/2))^2) &&    (Subscript[X, 1])== (1/2) && (Subscript[X, 2])>0.9 (δ/2)];

data={};
n=30
Protect[n];
For[α=(Subscript[α, 0]),α≤ (Subscript[α, 0])+ n Δα,α=α+Δα,
(Subscript[Γ, (Subscript[g, 4])])=DirichletCondition[(Subscript[U, 1])[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]==0,(Subscript[X, 1])==1/2 && (Subscript[X, 2])≥ -(δ/2)+α];
initBCs=InitializeBoundaryConditions[vd,sd,{{(Subscript[Γ, (Subscript[g, 1])]),(Subscript[Γ, (Subscript[g, 4])])},{(Subscript[Γ, (Subscript[g, 2])]),(Subscript[Γ, (Subscript[g, 5])])},{(Subscript[Γ, (Subscript[g, 3])])}}];
discreteBCs=DiscretizeBoundaryConditions[initBCs,methodData,sd];
{load,stiffness,damping,mass}=discretePDE["SystemMatrices"];
DeployBoundaryConditions[{load,stiffness},discreteBCs];
solution=LinearSolve[stiffness,load];

split=Span@@@Transpose[{Most[#+1],Rest[#]}&[methodData["IncidentOffsets"]]];
(Subscript[Un, 1])=ElementMeshInterpolation[{mesh}, solution[[split[[1]]]]];
(Subscript[Un, 2])=ElementMeshInterpolation[{mesh},  solution[[split[[2]]]]];
(Subscript[Un, 3])=ElementMeshInterpolation[{mesh},  solution[[split[[3]]]]];

Block[{λ=(ν/((1+ν)(1-2ν))),μ=(1/(2(1+ν)))},
(*( Ε (L^3) )*)
TotalEnergy=
 2(1/2)  NIntegrate[(μ (((Derivative[0,0,1][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+Derivative[0,1,0][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]))^2)+Derivative[0,0,1][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])] ((λ+2 μ) Derivative[0,0,1][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+λ Derivative[0,1,0][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+λ Derivative[1,0,0][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])])+Derivative[0,1,0][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])] (λ Derivative[0,0,1][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+(λ+2 μ) Derivative[0,1,0][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+λ Derivative[1,0,0][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])])+Derivative[1,0,0][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])] (λ Derivative[0,0,1][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+λ Derivative[0,1,0][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+(λ+2 μ) Derivative[1,0,0][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])])+μ (((Derivative[0,1,0][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+Derivative[1,0,0][(Subscript[Un, 2])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]))^2)+μ (((Derivative[0,0,1][(Subscript[Un, 1])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]+Derivative[1,0,0][(Subscript[Un, 3])][(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])]))^2)),{(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])}∈mesh];
];
AppendTo[data,{α-(Subscript[α, 0]) ,TotalEnergy}];]

data;
Π[δα_]=Fit[data, {1,δα,δα^2,δα^3,δα^4,δα^5,δα^6,δα^7,δα^8,δα^9,δα^10,δα^11}, δα];
G[δα_]=-D[Π[δα],δα];
Echo[UnitConvert[G[0](Ε (L^2)),"micro Newtons"],"Energy Release Rate :"];

Coeffs={−1/60, 3/20,−3/4, 0 ,3/4,−3/20, 1/60};
Echo[-UnitConvert[(Sum[data[[i,2]]Coeffs[[i]],{i,1,Length@Coeffs}]/Δα )(Ε (L^2)),"micro Newtons"],"Energy Release Rate :"];

Show[ListPlot[data],Plot[Π[δα],{δα,-3Δα,3Δα}]]
Show[Plot[{-(Sum[data[[i,2]]Coeffs[[i]],{i,1,Length@Coeffs}]/Δα ),G[δα]},{δα,-3Δα,3Δα}]]