MATH920 Lesson03: DNN-ELS

Least squares estimation of elastic moduli

Data input

Read data from solution of elasticity equations: get data file names, extract grid info, find number of time slices

SetDirectory["/home/student/courses/MATH920/models/DNN-ELS/interface/planestrain/data"];
qFiles=FileNames["q*.hdf"];  tFiles=FileNames["t*.hdf"];
tInfo = Map[Import[#,{"HDF","Data"}][[2]]&,tFiles];
nt = Length[tFiles]

Extract time values at each slice

t = tInfo[[;;,3]]

Extract grid step size (assume dx=dy=h, use grid spacing along x)

h =(tInfo[[1,14]]-tInfo[[1,10]])/tInfo[[1,6]]

Read data file with material properties used in BEARCLAW simulation.

setprob=ImportString[Import["../setprob.data","Text"],"Table"]

Extract numerical values of interest.

{cp1,cs1,ρ1,cp2,cs2,ρ2,x0,y0,r0,sig0} = Map[#[[1]]&,setprob]

Compute elastic moduli used for numerical simulation.

μ1=ρ1 (cs1^2); μ2=ρ2 (cs2^2); λ1=ρ1 (cp1^2)-2 μ1; ; λ2=ρ2 (cp2^2)-2 μ2;

Import the data files and extract number of time steps, components, grid sizes

q=Transpose[Table[Import[qFiles[[k]],{"HDF","Data"}],{k,nt}],{2,1,4,3}];
{nq,nt,m1,m2}=Dimensions[q]

Extract components of q into standard physics notation

{(Subscript[σ, 11]),(Subscript[σ, 22]),(Subscript[σ, 33]),(Subscript[σ, 12]),(Subscript[v, 1]),(Subscript[v, 2]),p}=q;

Define a subset of displacement values for visualization

xyuv=Table[Flatten[Table[{i1, i2,(Subscript[v, 1])[[k,i1,i2]], (Subscript[v, 2])[[k,i1,i2]]} ,{i1,1,m1,5},{i2,1,m2,5}],1],{k,nt}];

Visualize by plotting p=Tr(σ)/3 and displacement vectors

pdPlot[k_]:=Show[{ArrayPlot[Transpose[ p[[k]]],DataReversed->True ,PlotLegends->Automatic,ColorFunction->"Rainbow"],ListVectorPlot[xyuv[[k]]]}]

Investigate stress, displacement fields. Initial condition corresponds to a circular region in the lower half plane subjected to extension forces that are released at initial time. Hence material displaces initially towards circular region center.

pdPlot[4]

Process input data

(Subscript[X, 1])=Table[i1 h,{i1,m1},{i2,m2}];
(Subscript[X, 2])=Table[i2 h,{i1,m1},{i2,m2}];

Find displacements from velocity

zero=Table[0.,{i1,m1},{i2,m2}]; (Subscript[u, 1])={zero}; (Subscript[u, 2])={zero};
Monitor[For[k=2,k≤nt,k++,
         dt=(t[[k]]-t[[k-1]]);
         (Subscript[du, 1])= ((Subscript[v, 1])[[k-1]]+(Subscript[v, 1])[[k]]) dt/2; 
         (Subscript[du, 2])= ((Subscript[v, 2])[[k-1]]+(Subscript[v, 2])[[k]]) dt/2;
         AppendTo[(Subscript[u, 1]),(Subscript[u, 1])[[k-1]]+(Subscript[du, 1])];
         AppendTo[(Subscript[u, 2]),(Subscript[u, 2])[[k-1]]+(Subscript[du, 2])];
],k]

Define deformed coordinates, scaled for visualization

(Subscript[x, 1])={(Subscript[X, 1])}; (Subscript[x, 2])={(Subscript[X, 2])};
sc=2h/(Max[(Subscript[u, 1])[[nt]]]-Min[(Subscript[u, 1])[[nt]]]);
Monitor[
For[k=2,k≤nt,k++,
         AppendTo[(Subscript[x, 1]),(Subscript[X, 1])+sc*(Subscript[u, 1])[[k]]];
         AppendTo[(Subscript[x, 2]),(Subscript[X, 2])+sc*(Subscript[u, 2])[[k]]]; 
],k]

Define a function to plot grid lines

GridPlot[x_,y_]:= Module[{pts,xlines,ylines,mx,my},
{mx,my}=Dimensions[x];
pts =Table[{x[[i,j]],y[[i,j]]},{i,mx},{j,my}]; 
xlines =Map[Line[#]&,pts];
ylines= Map[Line[#]&,Transpose[pts]];
Graphics[{xlines,ylines}]
]

Animate[GridPlot[(Subscript[x, 1])[[k]],(Subscript[x, 2])[[k]]],{k,1,nt,1}]

Least squares problem

Compute pre-stress. Material starts with an initial isotropic stress (Subscript[σ, 11])=(Subscript[σ, 22])=(Subscript[σ, 33])=(Subscript[σ, 0]) in a circular region of radius r=0.125

(Subscript[σ, 0])=Table[If[ Sqrt[(((i1 h-1-x0))^2)+(((i2 h-1-y0))^2)]<r0,sig0,0],{i1,m1},{i2,m2}];

Compute strain and stress tensors at a specific time slice

k=4; d=2;
ϵ=Table[0,{j,d},{k,d},{i1,m1},{i2,m2}];
For[i1=2,i1<m1,i1++,
For[i2=2,i2<m2, i2++,
ϵ[[1,1,i1,i2]]=((Subscript[u, 1])[[k,i1+1,i2]]-(Subscript[u, 1])[[k,i1-1,i2]])/(2h);
ϵ[[1,2,i1,i2]]=((Subscript[u, 1])[[k,i1,i2+1]]-(Subscript[u, 1])[[k,i1,i2-1]]+
                                  (Subscript[u, 2])[[k,i1+1,i2]]-(Subscript[u, 2])[[k,i1-1,i2]])/(4h);
ϵ[[2,1,i1,i2]]= ϵ[[1,2,i1,i2]];
ϵ[[2,2,i1,i2]]=((Subscript[u, 2])[[k,i1,i2+1]]-(Subscript[u, 2])[[k,i1,i2-1]])/(2h);
]];
sig={{(Subscript[σ, 11])[[k]]-(Subscript[σ, 0]),(Subscript[σ, 12])[[k]]},{(Subscript[σ, 12])[[k]],(Subscript[σ, 22])[[k]]-(Subscript[σ, 0])}};

Form basis set over a region R((Subscript[i, 11]) to (Subscript[i, 12]),(Subscript[i, 21]) to (Subscript[i, 22])) of the simulation domain

i11=32; i12=48; i21=12; i22=28; δ=IdentityMatrix[2];
tr[a_]:=a[[1,1]]+a[[2,2]];
BR=Table[{tr[ϵ[[;;,;;,i1,i2]]] δ, 2 ϵ[[;;,;;,i1,i2]]},{i1,i11,i12},{i2,i21,i22}];
b1=Flatten[BR[[;;,;;,1,;;,;;]]];
b2=Flatten[BR[[;;,;;,2,;;,;;]]];
σR=Table[sig[[;;,;;,i1,i2]],{i1,i11,i12},{i2,i21,i22}];
s=Flatten[σR];
B=Transpose[{b1,b2}];
{Dimensions[B],Dimensions[s]}

Compute the QR decomposition and find the least squares solution

{Q,R}=QRDecomposition[B];
{λ,μ}=LinearSolve[R,Q.s]

Compare elastic moduli estimated by least squares with those used in simulation

TableForm[{{λ,λ1,λ2},{μ,μ1,μ2}},TableHeadings->{{"λ","μ"},{"Lsq","Mat 1","Mat 2"}}]

Define a function to investigate accuracy of least squares fit over various subregions of the domain

lsq[i10_,i20_,r_,k_]:= Module[{ϵ,i1,i2,i11,i12,i21,i22,d=2,sig},
ϵ=Table[0,{j1,d},{j2,d},{i1,m1},{i2,m2}];
For[i1=2,i1<m1,i1++,
For[i2=2,i2<m2, i2++,
ϵ[[1,1,i1,i2]]=((Subscript[u, 1])[[k,i1+1,i2]]-(Subscript[u, 1])[[k,i1-1,i2]])/(2h);
ϵ[[1,2,i1,i2]]=((Subscript[u, 1])[[k,i1,i2+1]]-(Subscript[u, 1])[[k,i1,i2-1]]+
                                  (Subscript[u, 2])[[k,i1+1,i2]]-(Subscript[u, 2])[[k,i1-1,i2]])/(4h);
ϵ[[2,1,i1,i2]]= ϵ[[1,2,i1,i2]];
ϵ[[2,2,i1,i2]]=((Subscript[u, 2])[[k,i1,i2+1]]-(Subscript[u, 2])[[k,i1,i2-1]])/(2h);
]];
sig={{(Subscript[σ, 11])[[k]]-(Subscript[σ, 0]),(Subscript[σ, 12])[[k]]},{(Subscript[σ, 12])[[k]],(Subscript[σ, 22])[[k]]-(Subscript[σ, 0])}};
i11=i10-r; i12=i10+r; i21=i20-r; i22=i20+r; δ=IdentityMatrix[d];
BR=Table[{tr[ϵ[[;;,;;,i1,i2]]] δ, 2 ϵ[[;;,;;,i1,i2]]},{i1,i11,i12},{i2,i21,i22}];
b1=Flatten[BR[[;;,;;,1,;;,;;]]];
b2=Flatten[BR[[;;,;;,2,;;,;;]]];
σR=Table[sig[[;;,;;,i1,i2]],{i1,i11,i12},{i2,i21,i22}];
s=Flatten[σR];
B=Transpose[{b1,b2}];
{Q,R}=QRDecomposition[B];
{λ,μ}=LinearSolve[R,Q.s];
TableForm[{{λ,λ1,λ2},{μ,μ1,μ2}},TableHeadings->{{"λ","μ"},{"Lsq","Mat 1","Mat 2"}}]
]

lsq[40,20,8,8]

lsqerr[i10_,i20_,r_,k_]:= Module[{ϵ,i1,i2,i11,i12,i21,i22,d=2,sig,yi20,err1,err2,err},
ϵ=Table[0,{j1,d},{j2,d},{i1,m1},{i2,m2}];
For[i1=2,i1<m1,i1++,
For[i2=2,i2<m2, i2++,
ϵ[[1,1,i1,i2]]=((Subscript[u, 1])[[k,i1+1,i2]]-(Subscript[u, 1])[[k,i1-1,i2]])/(2h);
ϵ[[1,2,i1,i2]]=((Subscript[u, 1])[[k,i1,i2+1]]-(Subscript[u, 1])[[k,i1,i2-1]]+
                                  (Subscript[u, 2])[[k,i1+1,i2]]-(Subscript[u, 2])[[k,i1-1,i2]])/(4h);
ϵ[[2,1,i1,i2]]= ϵ[[1,2,i1,i2]];
ϵ[[2,2,i1,i2]]=((Subscript[u, 2])[[k,i1,i2+1]]-(Subscript[u, 2])[[k,i1,i2-1]])/(2h);
]];
sig={{(Subscript[σ, 11])[[k]]-(Subscript[σ, 0]),(Subscript[σ, 12])[[k]]},{(Subscript[σ, 12])[[k]],(Subscript[σ, 22])[[k]]-(Subscript[σ, 0])}};
i11=i10-r; i12=i10+r; i21=i20-r; i22=i20+r; δ=IdentityMatrix[d];
BR=Table[{tr[ϵ[[;;,;;,i1,i2]]] δ, 2 ϵ[[;;,;;,i1,i2]]},{i1,i11,i12},{i2,i21,i22}];
b1=Flatten[BR[[;;,;;,1,;;,;;]]];
b2=Flatten[BR[[;;,;;,2,;;,;;]]];
σR=Table[sig[[;;,;;,i1,i2]],{i1,i11,i12},{i2,i21,i22}];
s=Flatten[σR];
B=Transpose[{b1,b2}];
{Q,R}=QRDecomposition[B];
{λ,μ}=LinearSolve[R,Q.s];
yi20=-1+(i20-0.5)h;
err1 =Sqrt[((((λ-λ1))^2)+(((μ-μ1))^2))/(((λ1)^2)+((μ1)^2))];
err2 =Sqrt[((((λ-λ2))^2)+(((μ-μ2))^2))/(((λ2)^2)+((μ2)^2))];
err=If[yi20<0,err2,err1];
Return[Log10[err]]
]

lsqerr[40,24,4,4]

Construct a plot of the relative error in elastic moduli identification by least squares

Off[LinearSolve::nosol]; Off[Set::shape];
fiterrors = Monitor[Table[{-1+(i1-0.5)h,-1+(i2-0.5)h,lsqerr[i1,i2,4,4]},{i1,8,m2,8},{i2,8,m2,8}],{i1,i2}];
ListContourPlot[Flatten[fiterrors,1],Contours->Table[l,{l,-3,0,0.5}],ContourLabels->All,PlotLegends->Automatic,PlotRange->All]

This is the relative error contour plot obtained for a 200 by 200 grid

This is the relative error contour plot obtained for an 80 by 80 grid