> restart; read(`NLElast.map`);

> 

######################################################################

# /home/mheil/maple/NLElast.map

#

# Derive the equations of nonlinear elasticity for selected undeformed 

# (reference) coordinate systems.

# 

# Notation: 

# ---------

# -- Symbolic vectors are expressed w.r.t. to a spatially fixed 

#    cartesian coordinate system. The cartesian index 

#    always comes first. 

# 

#       E.g.       rpos[i]    <-->   {\bf r}

#

#                  gvec[i,j]  <-->   {\bf g}_j

#

#

# -- Represent raised indices by neg. numbers

#

#       E.g.       x[-i]       <-->   x^i

#

#                  tau[-i,-j]  <-->   \tau^{ij}

#

# -- To make expressions more compact, the quantities in the

#    deformed configuration are also displayed in the following

#    abbreviated form

#

#          U1      =  u^1

#

#          U12     =  \frac{\partial u^1}{\partial x^2}

#

#          U123    =  \frac{\partial^2 u^1}{\partial x^2 \partial x^3}

#

#                etc.

#

#

######################################################################

> 

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> 

#-----------------------------------------------------

# Full output of even the most horrifying expressions?

#-----------------------------------------------------

> fulldisplay:=false:

> 

#--------------------------------------------------

# Set spatial dimension of problem and define

# the Lagrangian coordinates x^i = x[-i]

#--------------------------------------------------

> ndim:=2;

> x:=array(-3..-1):

> coords:=seq(x[-i],i=1..ndim);

> 

> 

#--------------------------------------------------

# Impose symmetry of stress tensor

#--------------------------------------------------

> tau[-2,-1]:=tau[-1,-2]:

> tau[-3,-1]:=tau[-1,-3]:

> tau[-3,-2]:=tau[-2,-3]:

> 

> 

#-----------------------------------------------------

# Choose coordinate system

#-----------------------------------------------------

> coordsyst:=1:

> 

> 

#-----------------------------------------------------

# Position vector in the undeformed coordinate system

#-----------------------------------------------------

> 

#-----------------------------------------------------

# Cartesian

#-----------------------------------------------------

> if (coordsyst=1) then 

>   x[-1]:=X;

>   x[-2]:=Y;

>   x[-3]:=Z;

>   assume(X,positive):

>   assume(Y,positive):

>   assume(Z,positive):

>   rpos[1](coords):=x[-1];

>   rpos[2](coords):=x[-2];

>   rpos[3](coords):=x[-3];

> fi;

#-----------------------------------------------------

# Cylindrical polars

#-----------------------------------------------------

> if (coordsyst=2) then 

>   x[-1]:=r;

>   x[-2]:=phi;

>   x[-3]:=z;

>   assume(r,positive):

>   rpos[1](coords):=x[-1]*cos(x[-2]);

>   rpos[2](coords):=x[-1]*sin(x[-2]);

>   rpos[3](coords):=x[-3];

> fi; 

#-----------------------------------------------------

# Spherical polars

#-----------------------------------------------------

> if (coordsyst=3) then 

>   x[-1]:=r;

>   x[-2]:=theta;

>   x[-3]:=phi;

>   assume(r,positive):

>   assume(theta,positive):

>   additionally(theta < Pi/2):

>   assume(phi,positive):

>   additionally(phi < Pi/2):

>   rpos[1](coords):=x[-1]*sin(x[-2])*cos(x[-3]);

>   rpos[2](coords):=x[-1]*sin(x[-2])*sin(x[-3]);

>   rpos[3](coords):=x[-1]*cos(x[-2]);

> fi; 

> 

> 

#----------------------------------------------------

# Covariant undeformed base vectors (tangent vectors)

# to coordinate lines)

#----------------------------------------------------

> gvec:=array(-ndim..ndim,-ndim..ndim):

> gvecaux:=array(1..ndim,1..ndim):

> 

> for i from 1 to ndim do 

> for j from 1 to ndim do

>   gvec[i,j]:=diff(rpos[i](coords),x[-j]);

>   gvecaux[i,j]:=gvec[i,j];

> od;

> od;

> 

> print(`Covariant undeformed base vectors (in columns of matrix)=`,gvecaux);

> 

#--------------------------------------------------

# Covariant undeformed metric tensor

#--------------------------------------------------

> gmet:=array(-ndim..ndim,-ndim..ndim):

> gmetaux:=array(1..ndim,1..ndim):

> gmetinv:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     gmet[i,j]:=simplify(add(gvec[k,i]*gvec[k,j],k=1..ndim));

>     gmetaux[i,j]:=gmet[i,j];

> od;

> od;

> 

> print(`Covariant undeformed metric tensor=`,gmetaux);

> 

> 

#--------------------------------------------------

# Contravariant metric tensor is the inverse

# of the covariant metric tensor

#--------------------------------------------------

> gmetinv:=inverse(gmetaux):

> print(`Contravariant undeformed metric tensor=`,gmetinv):

> 

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     gmet[-i,-j]:=simplify(gmetinv[i,j]);

> od;

> od;

> 

> 

#--------------------------------------------------

# Contravariant undeformed basis vectors

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do 

>     gvec[i,-j]:=add(gvec[i,k]*gmet[-k,-j],k=1..ndim);

>     gvecaux[i,j]:=gvec[i,-j]

> od;

> od;

> 

> 

> print(`Contravariant undeformed base vectors (in columns of matrix)=`,gvecaux)

> ;

> 

> 

#--------------------------------------------------

# Test contravariant undeformed basis vectors

#--------------------------------------------------

> CheckKronecker:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do 

>     CheckKronecker[i,j]:=simplify(add(gvec[ii,i]*gvec[ii,-j],ii=1..ndim));

> od;

> od;

> 

> print(`Test orthogonality of Co- and Contravariant undeformed base vectors=`,C

> heckKronecker);

> 

> 

#--------------------------------------------------

# Undeformed Christoffel symbols

#--------------------------------------------------

> Gamma:=array(1..ndim,1..ndim,-ndim..-1):

> 

> for i from 1 to ndim do

> for j from 1 to ndim do

> for k from 1 to ndim do

> 

>   Gamma[i,j,-k]:=simplify(add(

>      diff(gvec[ii,i],x[-j])*gvec[ii,-k],ii=1..ndim));

> 

> od;

> od;

> od;

> 

> print(`Undeformed Christoffel symbols=`,Gamma);

> 

> 

> 

#--------------------------------------------------

# Trafo to physical components

#--------------------------------------------------

> PhysCompList:=[seq(

>                b[-i]=simplify(b."*"[-i]/sqrt(gmet[i,i]),

>                symbolic),i=1..ndim),

>                seq(

>                u[-i](coords)=simplify(u."*"[-i](coords)/sqrt(gmet[i,i]),

>                symbolic),i=1..ndim),

>                seq(seq(simplify(

>                tau[-i,-j]=T."*"[-i,-j]*sqrt(gmet[-j,-j])/sqrt(gmet[i,i]),

>                symbolic),j=1..ndim),i=1..ndim)];

> 

> 

#--------------------------------------------------

# Equations of Elasticity in undef config

#--------------------------------------------------

> for i from 1 to ndim do

>   LinElast[i]:=eval(

>      add(

>           Diff(tau[-i,-j],x[-j]) + 

>           add(

>                Gamma[j,m,-i]*tau[-m,-j]+Gamma[j,m,-j]*tau[-i,-m],m=1..ndim)

>      ,j=1..ndim)

>      +rho*b[-i]);

> od;

> 

> 

> 

> 

#--------------------------------------------------

# Equations of Elasticity in undef config

# in physical components

#--------------------------------------------------

> for i from 1 to ndim do

>    LinElastPhys[i]:=simplify(subs(PhysCompList,LinElast[i]));

> od;

> 

> 

> 

###############################################################

# Deformed Configuration	

###############################################################

> 

> 

> 

#--------------------------------------------------

# Deformed Position vector: displacement decomposed

# into the undeformed basis.

#--------------------------------------------------

> for i from 1 to ndim do

>   uvec[i]:=add(u[-j](coords)*gvec[i,j],j=1..ndim);

> od;

> 

> for i from 1 to ndim do

>   Rpos[i](coords):=rpos[i](coords) + uvec[i];

> od;

> 

> 

> 

#--------------------------------------------------

# Covariant deformed Base vectors (tangent vectors)

# to deformed coordinate lines)

#--------------------------------------------------

> Gvec:=array(-ndim..ndim,-ndim..ndim):

> Gvecaux:=array(1..ndim,1..ndim):

> 

> for i from 1 to ndim do 

> for j from 1 to ndim do

>   Gvec[i,j]:=simplify(diff(Rpos[i](coords),x[-j]));

>   Gvecaux[i,j]:=Gvec[i,j]

> od;

> od;

> 

> print(`Covariant deformed base vectors (in columns of matrix)=`,Gvecaux);

> 

#--------------------------------------------------

# Covariant deformed metric tensor

#--------------------------------------------------

> Gmet:=array(-ndim..ndim,-ndim..ndim):

> Gmetaux:=array(1..ndim,1..ndim):

> Gmetinv:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     Gmet[i,j]:=simplify(add(Gvec[k,i]*Gvec[k,j],k=1..ndim));

>     Gmetaux[i,j]:=Gmet[i,j];

> od;

> od;

> 

> print(`Covariant deformed metric tensor=`,Gmetaux);

> 

> 

> 

#--------------------------------------------------

# Strain tensor

#--------------------------------------------------

> Strain:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     Strain[i,j]:=1/2*(Gmet[i,j]-gmet[i,j]);

> od;

> od;

> 

> print(`Full nonlinear strain tensor=`,Strain);

> 

> 

#--------------------------------------------------

# Linear Strain tensor

#--------------------------------------------------

> LinStrain:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     LinStrain[i,j]:=simplify(1/2*(add(

>                              gvec[ii,i]*diff(uvec[ii],x[-j])+

>                              gvec[ii,j]*diff(uvec[ii],x[-i]),ii=1..ndim)));

> od;

> od;

> 

> print(`Linear strain tensor=`,LinStrain);

> 

#--------------------------------------------------

# Linear Strain tensor in physical components

#--------------------------------------------------

> LinStrainPhys:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     LinStrainPhys[i,j]:=simplify(subs(PhysCompList,LinStrain[i,j]

>                         /sqrt(gmet[i,i]*gmet[j,j])));

> od;

> od;

> 

> print(`Linear strain tensor in phys. comp.=`,LinStrainPhys);

> 

> 

> 

#===================================================

# Test nonlinear Strain tensor (rigid body 

# rotation + displacement in cartesians) 

#===================================================

> StrainTest:=array(1..ndim,1..ndim):

> LinStrainTest:=array(1..ndim,1..ndim):

> if (coordsyst=1) then

>   utest[-1]:=simplify(U1+sqrt(x[-1]^2+x[-2]^2) * 

>              cos(arctan(x[-2]/x[-1])+Phi)-x[-1]);

>   utest[-2]:=simplify(U2+sqrt(x[-1]^2+x[-2]^2) * 

>              sin(arctan(x[-2]/x[-1])+Phi)-x[-2]);

>   utest[-3]:=0;

> 

>   for i from 1 to ndim do 

>   for j from 1 to ndim do

>      StrainTest[i,j]:=simplify(expand(eval(

>                       subs(u[-1](coords)=utest[-1],

>                            u[-2](coords)=utest[-2],

>                            u[-3](coords)=utest[-3],Strain[i,j]))));

>   od;

>   od;

>   print(`Nonlinear strain tensor for rigid body motion:`,StrainTest);

> 

> 

>   for i from 1 to ndim do 

>   for j from 1 to ndim do

>      LinStrainTest[i,j]:=simplify(

>                            subs(u[-1](coords)=utest[-1],

>                            u[-2](coords)=utest[-2],

>                            u[-3](coords)=utest[-3],LinStrain[i,j])

>                            ,symbolic);

>   od;

>   od;

> 

>   print(`Linear strain tensor for rigid body motion:`,LinStrainTest);

> 

> 

> 

>   for i from 1 to ndim do 

>   for j from 1 to ndim do

>      LinStrainTest[i,j]:=simplify(taylor(

>                            subs(u[-1](coords)=utest[-1],

>                            u[-2](coords)=utest[-2],

>                            u[-3](coords)=utest[-3],LinStrain[i,j]),

>                            Phi,3),symbolic);

>   od;

>   od;

> 

>   print(`Taylor expanded linear strain tensor for rigid body motion:`,LinStrai

> nTest);

> 

> 

> fi;

> ;

> 

> 

> 

#--------------------------------------------------

# Contravariant metric tensor is the inverse

# of the covariant metric tensor

#--------------------------------------------------

> Gmetinv:=inverse(Gmetaux):

> 

> if fulldisplay then

>   print(`Contravariant deformed metric tensor=`,Gmetinv):

> fi;

> 

> 

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     Gmet[-i,-j]:=simplify(Gmetinv[i,j]);

> od;

> od;

> 

> 

> 

#--------------------------------------------------

# Contravariant deformed basis vectors

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do 

>     Gvec[i,-j]:=simplify(add(Gvec[i,k]*Gmet[-k,-j],k=1..ndim));

>     Gvecaux[i,j]:=Gvec[i,-j];

> od;

> od;

> 

> 

> if fulldisplay then

>   print(`Contravariant deformed base vectors (in columns of matrix)=`,Gvecaux)

> ;

> fi;

> 

> 

#--------------------------------------------------

# Test contravariant deformed basis vectors

#--------------------------------------------------

> CheckKroneckerDef:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do 

>     CheckKroneckerDef[i,j]:=simplify(add(Gvec[ii,i]*Gvec[ii,-j],ii=1..ndim));

> od;

> od;

> 

> print(`Test orthogonality of Co- and Contravariant deformed base vectors=`,Che

> ckKroneckerDef);

> 

> 

#--------------------------------------------------

# Deformed Christoffel symbols (Note: The common

# factors in the denominator all come from

# the contravariant base vectors -- they carry

# through all the way to the final equations)

#--------------------------------------------------

> GammaDef:=array(1..ndim,1..ndim,-ndim..-1):

> 

> for i from 1 to ndim do

> for j from 1 to ndim do

> for k from 1 to ndim do

> 

>   GammaDef[i,j,-k]:=simplify(add(

>      diff(Gvec[ii,i],x[-j])*Gvec[ii,-k],ii=1..ndim));

> 

> od;

> od;

> od;

> 

> 

> if fulldisplay then

>   print(`Deformed Christoffel symbols=`,GammaDef);

> fi;

> 

> 

> 

#--------------------------------------------------

# Equations of Elasticity in def config

# (Note: The common factors in the denominator all 

# come from the Christoffel symbols and hence from

# the contravariant base vectors)

#--------------------------------------------------

> for ihide from 1 to 1 do

> for i from 1 to ndim do

>   Elast[i]:=eval(

>      add(

>           Diff(tau[-i,-j],x[-j]) + 

>           add(

>                GammaDef[j,m,-i]*tau[-m,-j]+GammaDef[j,m,-j]*tau[-i,-m],m=1..nd

> im)

>      ,j=1..ndim)

>      +rho*b[-i]):

> od;

> od;

> 

> if fulldisplay then

>   print(`Equations of Elasticity=`,Elast);

> fi;

> 

> 

##########################################################

# REWRITE EXPRESSIONS IN COMPACT FORM

##########################################################

> 

> oldquiet:=interface(quiet):

> oldecho:=interface(echo):

> 

> interface(echo=0,quiet=true);

> 

> 

#--------------------------------------------------

# Deformed Position vector: 

#--------------------------------------------------

> for i from 1 to ndim do

>   CompactRpos[i](coords):=rename(rpos[i](coords) + uvec[i]):

> od;

> 

> 

#--------------------------------------------------

# Deformed Base vectors

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do

>   Gvecaux[i,j]:=rename(Gvec[i,j]);

> od;

> od;

> 

> print(`Covariant deformed base vectors (in columns of matrix)=`,Gvecaux);

> 

#--------------------------------------------------

# Covariant metric tensor

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     Gmetaux[i,j]:=rename(Gmet[i,j]);

> od;

> od;

> 

> print(`Covariant deformed metric tensor=`,Gmetaux);

> 

> 

#--------------------------------------------------

# Strain tensor

#--------------------------------------------------

> CompactStrain:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     CompactStrain[i,j]:=rename(Strain[i,j]):

> od;

> od;

> 

> print(`Full nonlinear strain tensor=`,CompactStrain);

> 

> 

#--------------------------------------------------

# Linear Strain tensor

#--------------------------------------------------

> CompactLinStrain:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     CompactLinStrain[i,j]:=rename(LinStrain[i,j]);

> od;

> od;

> 

> print(`Linear strain tensor=`,CompactLinStrain);

> 

#--------------------------------------------------

# Linear Strain tensor in physical components

#--------------------------------------------------

> CompactLinStrainPhys:=array(1..ndim,1..ndim):

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     CompactLinStrainPhys[i,j]:=rename(LinStrainPhys[i,j]);

> od;

> od;

> 

> print(`Linear strain tensor in phys. comp.=`,CompactLinStrainPhys);

> 

> 

#--------------------------------------------------

# Contravariant metric tensor is the inverse

# of the covariant metric tensor

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do

>     Gmetaux[i,j]:=rename(Gmetinv[i,j]);

> od;

> od;

> 

> if fulldisplay then

>   print(`Contravariant deformed metric tensor=`,Gmetaux):

> fi;

> 

> 

#--------------------------------------------------

# Contravariant basis vectors

#--------------------------------------------------

> for i from 1 to ndim do 

> for j from 1 to ndim do 

>     Gvecaux[i,j]:=rename(Gvec[i,-j]);

> od;

> od;

> 

> 

> print(`Contravariant deformed base vectors (in columns of matrix)=`,Gvecaux);

> 

> 

> 

#--------------------------------------------------

# Deformed Christoffel symbols

#--------------------------------------------------

> CompactGammaDef:=array(1..ndim,1..ndim,-ndim..-1):

> 

> for i from 1 to ndim do

> for j from 1 to ndim do

> for k from 1 to ndim do

> 

>   CompactGammaDef[i,j,-k]:=rename(GammaDef[i,j,-k]);

> 

> od;

> od;

> od;

> 

> print(`Deformed Christoffel symbols=`,CompactGammaDef);

> 

> 

#--------------------------------------------------

# Equations of Elasticity in def config

#--------------------------------------------------

> for i from 1 to ndim do

>   CompactElast[i]:=collect(rename(Elast[i]),[tau[-1,-1],tau[-2,-2],tau[-1,-2]]

> );

> od;

> 

> 

>