Deal all,
In my example, I want to calculate e (uu, vv) '* D * e (uu, vv), but an error occurs.
Here’s my code
// Parameters
real E = 21.5;
real sigma = 0.29;
real gravity = -0.05;
// Mesh
border a(t=2, 0){x=0; y=t; label=1;}
border b(t=0, 10){x=t; y=0; label=2;}
border c(t=0, 2){ x=10; y=t; label=1;}
border d(t=0, 10){ x=10-t; y=2; label=3;}
mesh th = buildmesh(b(20) + c(5) + d(20) + a(5));
// Fespace
fespace Vh(th, [P1, P1]);
Vh [uu, vv];
Vh [w, s];
// Macro
fespace Zh(th,[P1,P1,P1]);
Zh [Sxx,Syy,Sxy];
real sqrt2 = sqrt(2.);
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2] //
macro div(u,v) (dx(u) + dy(v)) //
macro e(x1,x2) [dx(x1),0.5*(dx(x2)+dy(x1)), dy(x2)]//
// Problem
real mu = E/(2*(1+sigma));
real lambda = Esigma/((1+sigma)(1-2sigma));
real a11 = 2mu + lambda,a12 = 0,a13 = lambda;
real a22 = mu,a23 = 0;
real a33 = 2mu + lambda;
func D = [[a11,a12,a13],[a12,a22,a23],[a13,a23,a33]];
solve Elasticity ([uu, vv], [w, s])
= int2d(th)(
lambdadiv(w,s)*div(uu,vv)
+ 2.mu( epsilon(w,s)'epsilon(uu,vv) )
)
+ int2d(th)(
- gravitys
)
+ on(1, uu=0, vv=0)
;
[Sxx,Syy,Sxy]=D*e(uu,vv);
real[int] auxVec = (Sxx.n);
auxVec = Sxx;
real com = e(uu,vv)‘*Sxx;
// real c = Sxx’*auxVec;
//real k = e(uu,vv)'*Sxx;
// Movemesh
cout << auxVec << endl;