Welcome to the forums. Please post in English or French.

You are not logged in.

#1 2013-03-26 23:31:32

bashseb
Member
From: Autriche
Registered: 2012-03-28
Posts: 33
Website

How can I extract the right hand side of the Lagrange multiplier rows?

Dear all,

I'm playing around with extracting the stiffness matrix and the load vector just to understand FEA better smile . ... I can extract the stiffness matrix and the load vector using the command below.

My problem is, while the stiffness matrix contains the rows and columns of the (double) Lagrange multipliers, the load vector only contains the load for the physical d.o.f.

This is not needed for homogeneous Dirichlet b.c. because the rhs is zero. But  I'm interested also in the rhs of the Lagrange multipliers e.g. for the thermal problem or for LIAISON_DLL.

DEBUT()
import numpy
DEFI_MATERIAU..
MA=LIRE_MAILLAGE..
MO=AFFE_MODELE ..
CHMAT=AFFE_MATERIAU
CH0=AFFE_CHAR_MECA

# manual assembly of stiffness matrix
  ASSEMBLAGE(MODELE=MO,
                  CHAM_MATER=CHMAT,
                  # CARA_ELEM=CARELEM,
                  CHARGE=CH0,
                  NUME_DDL=CO('NUMEDDL'),
                  MATR_ASSE=(_F(MATRICE=CO('RIGIDITE'),
                                OPTION='RIGI_MECA',),
                             ),
                  VECT_ASSE=(_F(VECTEUR=CO("VECTASS"),
                      OPTION='CHAR_MECA')
                             ),
                  );
  # stiffness matrix
  rig=numpy.array(RIGIDITE.EXTR_MATR())
  numpy.set_printoptions(precision=6,suppress=True, linewidth=300,threshold=100000)

  print "Stiffness matrix size " + str(rig.shape) # this is including Lagrange multipliers
  rhs = VECTASS.EXTR_COMP() # the size of this vector is only the number of physical d.o.f. :(

thanks for any hints on this.

Last edited by bashseb (2013-04-08 14:38:35)

Offline

#2 2013-04-08 15:26:05

bashseb
Member
From: Autriche
Registered: 2012-03-28
Posts: 33
Website

Re: How can I extract the right hand side of the Lagrange multiplier rows?

Perhaps I should clarify what I mean...

Of the linear system `A.x=f` solved by MECA_STATIQUE, I want to extract `A` and `f` and store it in numpy variables for very small problems (1-10 elements). Then I call `numpy.linalg.solve(rig,load)` and compare the result to MECA_STATIQUE.

with the command above I can extract the stiffness matrix A and a vector (I think it is `f`). But the size of f is not the same as A. It is only defined at the physical nodes, whereas A, of course, contains rows and columns of the dualized Lagrange multipliers.

Here is the example a stiffness matrix of a plane, constrained single triangle (6 d.o.f. physical, 6 Lagrange multipliers)

Stiffness matrix size (12, 12)
 ['l0' 'l1' 'x1' 'y1' 'l2' 'l3' 'l4' 'x2' 'y2' 'l5' 'x3' 'y3']
[[-5.   0.   0.   5.   5.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.  -5.   5.   0.   0.   5.   0.   0.   0.   0.   0.   0. ]
 [ 0.   5.   7.5  2.5  0.   5.   0.  -5.  -2.5  0.  -2.5  0. ]
 [ 5.   0.   2.5  7.5  5.   0.   0.   0.  -2.5  0.  -2.5 -5. ]
 [ 5.   0.   0.   5.  -5.   0.   0.   0.   0.   0.   0.   0. ]
 [ 0.   5.   5.   0.   0.  -5.   0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   0.  -5.   0.   5.   5.   0.   0. ]
 [ 0.   0.  -5.   0.   0.   0.   0.   5.   0.   0.   0.   0. ]
 [ 0.   0.  -2.5 -2.5  0.   0.   5.   0.   2.5  5.   2.5  0. ]
 [ 0.   0.   0.   0.   0.   0.   5.   0.   5.  -5.   0.   0. ]
 [ 0.   0.  -2.5 -2.5  0.   0.   0.   0.   2.5  0.   2.5  0. ]
 [ 0.   0.   0.  -5.   0.   0.   0.   0.   0.   0.   0.   5. ]]
 ['l0' 'l1' 'x1' 'y1' 'l2' 'l3' 'l4' 'x2' 'y2' 'l5' 'x3' 'y3']

but `f` looks like this (three dof have been assigned with FORCE_NODALE):

[ 0.   0.   0.4  0.  -0.2  1.   0.   0. ]

So I use a trick. I noticed that the Lagrange multipliers result in a negative number on the diagonal.

# stiffness matrix is stored in `rig`
nodeValsInRig = rig.diagonal() >=0
# set up load vector
nddl=rig.shape[0]
load = numpy.zeros(nddl)
rhs = VECTASS.EXTR_COMP()
rhsval= rhs.valeurs
load[nodeValsInRig] = rhsval # write rhs only in physical positions of the load vector (leave lagrange nodes 0)
# it looks like this when I apply FORCE_NODALE on three physical d.o.f.:
# [ 0.   0.   0.   0.   0.   0.   0.   0.4  0.   0.   0.  -0.2  1.   0.   0.   0. ] 

# then, the following command yields the same result as MECA_STATIQUE
x = numpy.linalg.solve(rig,load)

This works when I use FORCE_NODALE in AFFE_CHAR_MECA and even for the LIAISON commands.

However, `f` is not always zero in the Lagrange multiplier rows as assumed above, right? For example, when I use LIAISON_DDL, the above works only with COEF_IMPO=0. For COEF_IMPO=1. the numpy result and the meca_statique result differ. The reason must be my incorrect `load` vector `f`.

So I'd like to know how I can extract the true `f` of the linear system without resorting to tricks that don't work in special cases. The complete comm file is attached. Many thanks in advance!

Last edited by bashseb (2013-04-09 08:22:24)


Attachments:
simple.comm, Size: 7.85 KiB, Downloads: 229

Offline

#3 2020-08-14 03:58:55

Nidish
Member
Registered: 2020-06-04
Posts: 18

Re: How can I extract the right hand side of the Lagrange multiplier rows?

Hello,

I'm very much interested in exactly this. Did you manage to find an answer to this yet?

Thank you,
Nidish

Offline

#4 2020-08-15 21:37:36

bashseb
Member
From: Autriche
Registered: 2012-03-28
Posts: 33
Website

Re: How can I extract the right hand side of the Lagrange multiplier rows?

Hello Nidish,

I'm sorry, but I couldn't figure this out -- I actually remember the topic, but I cannot understand the code_aster commands any more and hence don't know exactly what I even did seven years ago ... this was the first code_aster email I received in years.

Hope it helps nonetheless and another forum member can answer our question.

All the best!

Offline