How to calculate energy norm in linear elasticity promblem?

It is defined in terms of exact stress σ, computed stress σ^h, and coefficient matrix D written in Voigt notation as:

``````//Parameters
real Rho = 8000.;		//Density
real E = 210.e9;		//Young modulus
real Nu = 0.27;			//Poisson ratio

real Gravity = -9.81;	//Gravity

//Mesh
real nn = 10;			//Mesh quality
real L = 20.;			//Beam length
real H = 1.;			//Beam height
int Fixed = 1;			//Beam fixed label
int Free = 2;			//Beam free label
border b1(t=0., L){x=t; y=0.; label=Free;};
border b2(t=0., H){x=L; y=t; label=Fixed;};
border b3(t=0., L){x=L-t; y=H; label=Free;};
border b4(t=0., H){x=0.; y=H-t; label=Fixed;};

int nnL = max(2., nn*L);
int nnH = max(2., nn*H);
mesh Th = buildmesh(b1(nnL) + b2(nnH) + b3(nnL) + b4(nnH));

//Fespace
func Pk = P2;
fespace Uh(Th, [Pk, Pk]);
Uh [ux, uy];

//Macro
real sqrt2 = sqrt(2.);
macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //

//Problem
real Mu = E/(2.*(1.+Nu));
real Lambda = E*Nu/((1.+ Nu)*(1.-2.*Nu));

varf vElasticity ([ux,uy], [vx, vy])
= int2d(Th)(
Lambda * Divergence(vx, vy) * Divergence(ux, uy)
+ 2. * Mu * (
Epsilon(vx, vy)' * Epsilon(ux, uy)
)
)
+ int2d(Th)(
Rho * Gravity * vy
)
+ on(Fixed, ux=0, uy=0)
;

matrix<real> Elasticity = vElasticity(Uh, Uh, solver=sparsesolver);
real[int] ElasticityBoundary = vElasticity(0, Uh);
ux[] = Elasticity^-1 * ElasticityBoundary;

//Movemesh
Th = movemesh(Th, [x+ux, y+uy]);
[ux, uy] = [ux, uy];

//Plot
plot([ux, uy], value=true, cmm="u");
``````

In this 2d linear elasticity code, I want to calculate the energy norm as in the above expression.

ok, they is no difficulty

if the displacement is P2 [ux,uy] the discret stress is P1dc by construction

and the discrete stress is

``````fespace Wh3(Th,[P1dc, P1dc, P1dc]) ;
Wh3 [sigmaxx, sigmayy, sigmaxy] ;
func sigmah = [sigmaxx, sigmayy, sigmaxy] ;
func I = [1.,1.,0.];
sigmah = (Lambda * Divergence(ux, uy))*I + 2*Mu* Epsilon(ux,uy);
``````

now you have acess to `exact sigma` et `discret sigmah` so you can compute you error.

1 Like