diff --git a/benchmarks/elasticity/conge2.geo b/benchmarks/elasticity/conge2.geo
new file mode 100644
index 0000000000000000000000000000000000000000..6edcb31d88af71ecdd0d1f0238aeedf4532d900c
--- /dev/null
+++ b/benchmarks/elasticity/conge2.geo
@@ -0,0 +1,93 @@
+
+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};
+
diff --git a/benchmarks/elasticity/conge2.lua b/benchmarks/elasticity/conge2.lua
new file mode 100644
index 0000000000000000000000000000000000000000..d159df2bc882bb12d5da4d26a401ec730affdbc2
--- /dev/null
+++ b/benchmarks/elasticity/conge2.lua
@@ -0,0 +1,64 @@
+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)
diff --git a/benchmarks/elasticity/cylinder.geo b/benchmarks/elasticity/cylinder.geo
new file mode 100644
index 0000000000000000000000000000000000000000..86109a2e55a6136db726955d903ad773ef02a47c
--- /dev/null
+++ b/benchmarks/elasticity/cylinder.geo
@@ -0,0 +1,33 @@
+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};
diff --git a/benchmarks/elasticity/cylinder.lua b/benchmarks/elasticity/cylinder.lua
new file mode 100644
index 0000000000000000000000000000000000000000..c4192d904accaa9dc3f2d120e897a73957d37ac9
--- /dev/null
+++ b/benchmarks/elasticity/cylinder.lua
@@ -0,0 +1,59 @@
+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)
+