-- FreeFem++ v4.12 (Sat 06 May 2023 07:56:14 PM CEST - git v4.12-72-gb0015903)
   file : ./Elasticity-simple-support-BC.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : /*
    2 :   Exemple of homogenious Dirichket Boundary condition on some Edge in 3d ..
    3 :  in linear elasticity.
    4 :  The two edge are intersection of face with label 2,5 (resp.  4,5)
    5 :  
    6 : */
    7 : load "msh3" 
    8 : mesh3 Th=cube(20,5,5,[x*4,y,z]);
    9 : fespace Wh(Th,[P2,P2,P2]) ;
   10 : // appuie simple sur les bord interfece face 5 et4 ( arete x == 0 )   et 5 et 2  et (x == 1) ( arete x == 0 ) 
   11 : 
   12 : real E = 21.5e4;
   13 : real sigma = 0.29;
   14 : real mu = E/(2*(1+sigma));
   15 : real lambda = E*sigma/((1+sigma)*(1-2*sigma));
   16 : real gravity = -100;
   17 : 
   18 : real sqrt2=sqrt(2.);
   19 : 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
   20 : macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) )  )  // EOM
   21 : Wh [u1,u2,u3],[v1,v2,v3];
   22 : varf von5([u1,u2,u3],[v1,v2,v3]) = on(5,u1=1,u2=1,u3=1);
   23 : varf von24([u1,u2,u3],[v1,v2,v3]) = on(2,4,u1=1,u2=1,u3=1);
   24 : 
   25 : Wh [au1,au2,au3];
   26 : {
   27 : 	real[int] w5=von5(0,Wh, tgv=1); //  find dof on face 5 
   28 : 	real[int] w24=von24(0,Wh, tgv=1); //  find dof on face 2 and 4 
   29 : 	au1[] = w5 .* w24; // 1 do intersect  face 5 and 2 4 
   30 : 	plot(au3,wait=1); // see for debugging 
   31 : }
   32 : // so array au1[] is non zero on dof on edges  (2,5) and (4,5)
   33 : 
   34 : varf  Lame([u1,u2,u3],[v1,v2,v3])=
   35 :   int3d(Th)(  
   36 : 	     lambda*div(u1,u2,u3)     ( dx(u1)+dy(u2)+dz(u3) ) *div(v1,v2,v3)     ( dx(v1)+dy(v2)+dz(v3) ) 	
   37 : 	    +2.*mu*( 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] '*epsilon(v1,v2,v3)     [dx(v1),dy(v2),dz(v3),(dz(v2)+dy(v3))/sqrt2,(dz(v1)+dx(v3))/sqrt2,(dy(v1)+dx(v2))/sqrt2]  ) //')
   38 : 	      )
   39 :   + int3d(Th) (gravity*v3)
   40 :   
   41 :   ;
   42 : 
   43 :   matrix A = Lame(Wh,Wh,sym=1,positive=1,solver="CG");
   44 :   cout << " half "<< A.half << " nnz "<< A.nnz << endl; 
   45 :   real[int] b= Lame(0,Wh); 
   46 :   //  put tgv = -2 on matrix A 
   47 :   setBC(A,au1[],-2);   // 1 on diagonal and  0 on row  i et column i if  au1[][i] !=0  
   48 :   b =  au1[] ? 0 : b;
   49 :   set(A,solver= "CHOLMOD"); 
   50 :   cout << A.nnz << endl; 
   51 :    u1[] = A^-1*b;
   52 :  cout << " || u|| =" << u1[].linfty << endl;
   53 :  real cc = 0.5/u1[].linfty;
   54 :  mesh3 Thm = movemesh(Th,[x+u1*cc,y+u2*cc,z+u3*cc]) ;
   55 :  plot(Th,Thm,wait=1);
   56 :   
   57 : 
   58 :  sizestack + 1024 =3512  ( 2488 )

  -- Cube  nv=756 nt=3000 nbe=900 kind= 6
  -- Build Nodes/DF on mesh :   n.v. 756, n. elmt. 3000, n b. elmt. 900
     nb of Nodes 4961    nb of DoF   14883  DFon=3300
  -- FESpace: Nb of Nodes 4961 Nb of DoF 14883
  Plot::  Sorry no ps version for this type of plot 6
 half 1 nnz 572601
572601
CHOLMOD warning: matrix not positive definite. file: ../Supernodal/t_cholmod_super_numeric.c line: 911
cholmod error: file: ../Supernodal/t_cholmod_super_numeric.c line: 911 status: 1: matrix not positive definite

 || u|| =nan
 Bad orientation of tet  before cleanmesh 0 mes = 0
  current line = 54
Assertion fail : (l)
	line :422, in file ../femlib/GQuadTree.cpp
Assertion fail : (l)
	line :422, in file ../femlib/GQuadTree.cpp
 err code 6 ,  mpirank 0