Added examples elasticity

benchmarks/elasticity/conge2.geo

+
+Mesh.CharacteristicLengthFromCurvature = 0;
+Mesh.CharacteristicLengthFromPoints = 1;
+Mesh.CharacteristicLengthExtendFromBoundary = 1;
+
+unit = 1.0e-02 ;
+
+e1 =  4.5 * unit ;
+e2 =  6.0 * unit / 2.0 ;
+e3 =  5.0 * unit / 2.0 ;
+h1 =  5.0 * unit ;
+h2 = 10.0 * unit ;
+h3 =  5.0 * unit ;
+//h4 =  1.0 * unit ;
+h4 =  2.0 * unit ;
+h5 =  4.5 * unit ;
+R1 =  1.0 * unit ;
+R2 =  1.5 * unit ;
+//R2 =  1.5 * unit ;
+r  =  2 * unit ;
+ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
+ssin = ( 1.0 - ccos*ccos )^0.5 ;
+
+eps = 0.01 * unit;
+
+Lc1 = 0.001 ;
+Lc2 = 0.001 ;
+
+Point(1) = { -e1-e2, 0.0  , 0.0 , Lc1};
+Point(2) = { -e1-e2, h1   , 0.0 , Lc1};
+Point(3) = { -e3-r , h1   , 0.0 , Lc2};
+Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
+Point(5) = { -e3   , h1+r , 0.0 , Lc2};
+Point(6) = { -e3   , h1+h2-eps, 0.0 , Lc1};
+Point(7) = {  e3   , h1+h2, 0.0 , Lc1};
+Point(8) = {  e3   , h1+r , 0.0 , Lc2};
+Point(9) = {  e3+r , h1+r , 0.0 , Lc2};
+Point(10)= {  e3+r , h1   , 0.0 , Lc2};
+Point(11)= {  e1+e2, h1   , 0.0 , Lc1};
+Point(12)= {  e1+e2, 0.0  , 0.0 , Lc1};
+Point(13)= {  e2   , 0.0  , 0.0 , Lc1};
+
+Point(14)= {  R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
+Point(15)= {  0.0       , h5        , 0.0 , Lc2};
+Point(16)= { -R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
+Point(17)= { -e2        , 0.0       , 0.0 , Lc1};
+
+Point(18)= { -R2  , h1+h3   , 0.0 , Lc2};
+Point(19)= { -R2  , h1+h3+h4, 0.0 , Lc2};
+Point(20)= {  0.0 , h1+h3+h4, 0.0 , Lc2};
+Point(21)= {  R2  , h1+h3+h4, 0.0 , Lc2};
+Point(22)= {  R2  , h1+h3   , 0.0 , Lc2};
+Point(23)= {  0.0 , h1+h3   , 0.0 , Lc2};
+
+Point(24)= {  0 , h1+h3+h4+R2, 0.0 , Lc2};
+Point(25)= {  0 , h1+h3-R2,    0.0 , Lc2};
+
+Line(1)  = {1 ,17};    /* ux=uy=0 */
+Line(2)  = {17,16};
+Circle(3) = {16,15,14};
+Line(4)  = {14,13};
+Line(5)  = {13,12};    /* ux=uy=0 */
+Line(6)  = {12,11};
+Line(7)  = {11,10};
+Circle(8) = { 8, 9,10};
+Line(9)  = { 8, 7};
+Line(10) = { 7, 6};    /* T=10000 N */
+Line(11) = { 6, 5};
+Circle(12) = { 5, 4, 3};
+Line(13) = { 3, 2};
+Line(14) = { 2, 1};
+
+Line(15) = {18,19};
+Circle(16) = {24,20,21};
+Circle(17) = {19,20,24};
+Circle(18) = {18,23,25};
+Circle(19) = {22,23,25};
+Line(20) = {21,22};
+
+Line Loop(21) = {17,15,-18,19,20,16};
+Line Loop(23) = {11,12,13,14,1,2,3,4,5,6,7,-8,9,10};
+Plane Surface(24) = {23,21};
+
+//Physical Line(25) = {9,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,19,20,10};
+//Physical Surface(26) = {22,24};
+Physical Line(9) = {11};
+Physical Line(10) = {10};
+Physical Line(8) = {1, 5};
+Physical Surface(7) = {24};
+Physical Point(20) = {2};
+Physical Point(21) = {11};
+// Physical Line(25) = {11, 12, 13, 14};
+

benchmarks/elasticity/conge2.lua

+function neumannCondition (x,y,z)
+  return {-100e6,0,0}
+end
+
+
+function errorOnView (view, dim)
+  basis = polynomialBasis.find(2)
+  gaussPoints = fullMatrix(1,1)
+  gaussWeights = fullMatrix(1,1)
+  gaussIntegration.getTriangle(3,gaussPoints,gaussWeights)
+
+  nnodes = 3 -- only triangles
+  viewData = view:getData()
+  values = fullMatrix(nnodes,dim)
+  nodes = fullMatrix(nnodes,3)
+  gaussXYZ = fullMatrix(gaussPoints:size1(),3)
+  gaussValues = fullMatrix(gaussPoints:size1(),dim)
+  nEntities = viewData:getNumEntities(-1)
+  psi = fullMatrix(gaussPoints:size1(),3)
+  basis:f(gaussPoints,psi)
+  for i=0,nEntities-1 do 
+    nElements = viewData:getNumElements(-1,i)
+    for j=0,nElements-1 do
+      if(viewData:getNumNodes(0,i,j)==nnodes) then
+        viewData:getAllValuesForElement(0,i,j,values)
+        viewData:getAllNodesForElement(0,i,j,nodes)
+      end
+      gaussXYZ:gemm(psi,nodes,1,0)
+      gaussValues:gemm(psi,values,1,0)
+    end
+  end
+end
+
+
+
+m = GModel()
+m:load("conge.geo")
+m:load("conge.msh")
+
+print('---> creating elastic solver')
+e = elasticitySolver(m,1)
+e:addElasticDomain (7, 200e9, 0.3)
+e:addDirichletBC (1,8,0,0)
+e:addDirichletBC (1,8,1,0)
+e:addDirichletBC (1,8,2,0)
+e:addNeumannBCLua (1,9,"neumannCondition")
+sys = linearSystemPETSc()
+--sys = linearSystemCSRTaucs()
+
+print('----------> assembling')
+e:assemble(sys)
+
+print('----------> solving')
+--sys:systemSolve()
+e:solve()
+
+print('----------> writing results')
+format = 5 --gmsh format
+--view = e:buildVonMisesView("vonMises")
+--view = e:buildDisplacementView("displacement")
+view = e:buildStressesView("stresses")
+view:write("post.msh",format,false)
+
+errorOnView(view,3)

benchmarks/elasticity/cylinder.geo

+Mesh.CharacteristicLengthFactor=0.5;
+
+R = 0.3;
+H = 1.5;
+
+Point (1) = {0, 0, 0};
+Point (2) = {0, R, 0};
+Point (3) = {0, -R, 0};
+Point (4) = {0, 0, R};
+Point (5) = {0, 0, -R};
+
+Circle (1) = {2, 1, 4};
+Circle (2) = {4, 1, 3};
+Circle (3) = {3, 1, 5};
+Circle (4) = {5, 1, 2};
+
+Line Loop(5) = {4,1,2,3};
+Plane Surface(6) = {5};
+
+//here we create a volume by extruding
+Extrude {H,0,0} {
+  Surface{6};
+}
+
+//Compound Surface(100)={15,19,23,27};
+
+//Surface Loop(200)={100,6,28};
+//Volume(201)={200};
+
+Physical Surface(100)={6};
+Physical Surface(200)={28};
+Physical Surface(300)={15,19,23,27};
+Physical Volume(1)={1};

benchmarks/elasticity/cylinder.lua

+function neumann( XYZ, NXYZ, FCT )
+  for i=0,FCT:size1()-1 do
+    nx = NXYZ:get(i,0)
+    ny = NXYZ:get(i,1)
+    nz = NXYZ:get(i,2)
+    print('nx= ', nx, 'ny= ', ny, 'nz= ',nz);
+    FCT:set(i,0, nx)
+    FCT:set(i,1, ny)
+    FCT:set(i,2, nz) 
+  end
+end
+
+E = 200e5
+nu = 0.3
+xyz = functionCoordinates.get()
+nxyz = functionNormals.get()
+
+NEUM = functionLua(3, 'neumann', {xyz,nxyz})
+
+
+m = GModel()
+m:load("cylinder.geo")
+m:load("cylinder.msh")
+e = elasticitySolver(m,1)
+e:addElasticDomain (6, E, nu)
+
+print('---> creating elastic solver')
+e = elasticitySolver(m,1)
+e:addElasticDomain (1, E, nu)
+e:addDirichletBC (2,100,0,1)
+e:addDirichletBC (2,100,1,0)
+e:addDirichletBC (2,100,2,0)
+e:addDirichletBC (2,200,0,0)
+e:addDirichletBC (2,200,1,0)
+e:addDirichletBC (2,200,2,0)
+--e:addNeumannBC (2,300,{10,0,0})
+--e:addNeumannBC (2,300,{0,0,0})
+--e:addNeumannBCFct (2,300,NEUM)
+
+
+sys = linearSystemCSRTaucs()
+
+print('----------> assembling')
+e:assemble(sys)
+
+print('----------> solving')
+--sys:systemSolve()
+e:solve()
+
+print('----------> writing results')
+format = 5 --gmsh format
+--view = e:buildVonMisesView("vonMises")
+--view = e:buildDisplacementView("displacement")
+view = e:buildStressesView("stresses")
+view:write("post.msh",format,false)
+
+view2 = e:buildDisplacementView("displacements")
+view2:write("disp.msh",format,false)
+