Solve U=K^-1.F using stiffness matrix for vector fields

General Discussion
1

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 = whh*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 = Enu/((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+uucoef, y+vvcoef, z+ww*coef], label=ref2);

Thm = change(Thm, label=ref2);

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

image
image907×540 7.47 KB

2

Hello,
What are you trying to do exactly? Computing K^-1 explicitly should be avoided at all cost.

3

Hi prj,

My final goal is to use the Newmark-beta method for transient problems in elasticity.
So in order to do that, I need the Mass and Stiffness matrices for the body and operate on the same.
The equations to be solved consist of the inverse of these matrices at final.
To get the idea, I was first trying with this static deflection problem of beam.
So how should I do that in case K^-1 should not be used?
I am attaching the image of one of the equation.
[m]-mass matrix, [k] - stiffness matrix

newmark

4

In your code you have:

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

What do you do that, and not:

U=K^-1.F
5

In my code,
Kinv is an array real[int, int] which contains the elements of the inverse of stiffness matrix.
F is an array real[int] which contains the force vector for the whole body.
U is an array containing the displacements.
To carry out the matrix-vector multiplication, which is compatible in the way I use, I have used Kinv. So further, KU=F, which gives U = K^-1F = Kinv*F.

Also, can you let me know what problems can occur when I take an inverse of the matrix obtained from the varf.
The method I follow is:
I store the elements of obtained matrix in an array real[int,int] and then use lapack to get the inverse.

6

But why do you need the inverse explicitly, and why can’t you just compute U=A^-1*F, using matrix A assembled from your varf?

7

For the Newmark method that I mentioned earlier, I will require the inverse of addition of matrices explicitly, that is:
[ aM + bC + K ]^-1
where,
M - Mass matrix for the body
C - Damping Matrix
K - Stiffness matrix of the body
a,b,c - real-constants
For doing this, I am not able to create a form to be fed to varf.
This is the reason I was searching for an explicit method. It is a time-stepping method and I need to loop it.
Can you let me know any other way for this method?

8

What is your damping matrix? I don’t see the problem assembling through a varf the term aM + K.
You can then compute the sum of both matrix, and use set(sum, solver = sparsesolver).

9

At present, the damping matrix is zero, but I will be considering a model as follows for damping matrix:
C = nM + mK;
where, M,C,K represent the same matrices as earlier and m,n are real constants,
From your reply, as far as I understand, I should compute the K and M matrices using varf, then store their addition as required in a matrix sum = aM + b C + K, and then use
set(sum, solver = sparsesolver)
Am I correct?

10

You are indeed correct.

11

I was able to execute the whole code much faster using the above way.
Also, I rectified an error regarding the FE array size in my earlier code.
Now my code seems to work completely as expected.
Thanks a lot, @prj!

12

What are uucoef and vvcoef in the OP?

13

Those are the coefficients for amplification in visualization, they do not affect the results but help in visualization by amplifying the displacements.

1 Like