Dear all,

I am trying to solve problem of 3-dimensional elasticity. I have extracted the stiffness matrix (K) and force vector (F) using the variational formulation. Now, I have solved the equation for displacement U, using U=K^-1.F.

Now, displacement is a vector field, so I need to again distribute U into its respective forms, U1, U2, U3. I did it but the plot is not what I should get for a tip-loaded beam.

Please find my code below. Thanks in advance.

load “mshmet”

load “msh3”

load “medit”

load “TetGen”

include “cube.idp”

load “lapack”

real l = 10; // lenght

real w = 1; // breadth

real h = 1; // height

real I = w*h*h*h/12; // area moment of inertia

int[int] Nxyz=[20,2,2];

real [int,int] Bxyz=[[0.,l],[0.,w],[0.,h]];

int [int,int] Lxyz=[[1,2],[3,4],[5,6]];

mesh3 Th=Cube(Nxyz,Bxyz,Lxyz);

real E = 210.0e9;

real nu = 0.3;

real rho = 7850;

real mu = E/(2*(1+nu));

real lambda = E*nu/((1+nu)*(1-2*nu));

// Fespace

fespace Vh(Th,[P1,P1,P1]);

Vh [u1,u2,u3], [v1,v2,v3];

real sqrt2=sqrt(2.);

macro epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM

macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

varf a([u1,u2,u3],[v1,v2,v3],tgv=1e50)=

```
int3d(Th)(
lambda*div(u1,u2,u3)*div(v1,v2,v3)
+2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) )
)
+ on(1,u1=0,u2=0,u3=0)
```

;

varf f([u1,u2,u3],[v1,v2,v3])=

```
int2d(Th,2)(60000*v3);
```

matrix A= a(Vh,Vh,solver=“SPARSESOLVER”); //stiffness matrix

real[int] F = f(0,Vh);

real[int,int] K(A.m,A.m);

real c = 0;

for (int i=0; i<A.m;i++){

```
for(int j=0; j<A.m; j++){
K(i,j)=A(i,j);
}
```

}

real[int,int] a1realinv(A.m,A.n);

a1realinv = K^-1;

matrix Kinv(A.m,A.n); /*inverse of matrix a1*/

for (int e=0; e<A.m; e++){

```
for (int ee=0;ee<A.n; ee++){
Kinv(e,ee) = a1realinv(e,ee);
}
```

}

real[int] U(A.n); /*Displacement*/

U=Kinv*F; //U=K^-1.F

real[int] U1(A.m/3),U2(A.m/3),U3(A.m/3);

for(int e=0; e<A.m/3; e++){

```
U1(e) = U(3*e);
U2(e) = U(3*e+1);
U3(e) = U(3*e+2);
```

}

real coef = 1000;

int[int] ref2 = [1, 0, 2, 0];

Vh [uu,vv,ww];

uu[] = U1;

vv[] = U2;

ww[] = U3;

mesh3 Thm = movemesh3(Th, transfo=[x+uu*coef, y+vv*coef, z+ww*coef], label=ref2);

Thm = change(Thm, label=ref2);

plot(Thm, wait=true, cmm="coef amplification = "+coef);