diff --git a/benchmarks/elasticity/conge.geo b/benchmarks/elasticity/conge.geo
deleted file mode 100644
index ca3c052d2856abd298f7f2a674a619e5f0f3df6c..0000000000000000000000000000000000000000
--- a/benchmarks/elasticity/conge.geo
+++ /dev/null
@@ -1,94 +0,0 @@
-
-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) = {14,15,16};
-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) = { 3, 4, 5};
-Line(13) = { 3, 2};
-Line(14) = { 2, 1};
-
-Line(15) = {18,19};
-Circle(16) = {21,20,24};
-Circle(17) = {24,20,19};
-Circle(18) = {18,23,25};
-Circle(19) = {25,23,22};
-Line(20) = {21,22};
-
-Line Loop(21) = {17,-15,18,19,-20,16};
-//Plane Surface(22) = {21};
-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/conge.lua b/benchmarks/elasticity/conge.lua
deleted file mode 100644
index ee8860733889f1b9ae6c418719a50b7902707235..0000000000000000000000000000000000000000
--- a/benchmarks/elasticity/conge.lua
+++ /dev/null
@@ -1,53 +0,0 @@
-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")
-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 = linearSystemCSRTaucs()
-e:assemble(sys)
-sys:systemSolve()
-view = e:buildVonMisesView("vonMises")
-
-view = e:buildDisplacementView("displacement")
-view:write("displacement.msh",5,false)
-
-errorOnView(view,3)
diff --git a/benchmarks/elasticity/hole.geo b/benchmarks/elasticity/hole.geo
deleted file mode 100644
index 9187cf9f3d8dfe1374a78c50649bc6334deff7f8..0000000000000000000000000000000000000000
--- a/benchmarks/elasticity/hole.geo
+++ /dev/null
@@ -1,39 +0,0 @@
-lcExt = 0.1;
-lcIn = 0.01;
-a=0.05;
-L=3;
-
-Point(1)={-L/2,-L/2,0,lcExt};
-Point(2)={L/2,-L/2,0,lcExt};
-Point(3)={L/2,L/2,0,lcExt};
-Point(4)={-L/2,L/2,0,lcExt};
-
-Line(1)={1,2};
-Line(2)={2,3};
-Line(3)={3,4};
-Line(4)={4,1};
-
-Point(10) = {0,0,0,lcIn};
-Point(11) = {a,0,0,lcIn};
-Point(12) = {0,a,0,lcIn};
-Point(13) = {-a,0,0,lcIn};
-Point(14) = {0,-a,0,lcIn};
-Circle (10)= {11,10,12};
-Circle (11)= {12,10,13};
-Circle (12)= {13,10,14};
-Circle (13)= {14,10,11};
-
-Line Loop(1)={1,2,3,4};
-Line Loop(2)={10,11,12,13};
-Plane Surface(1) = {1,2};
-
-Physical Line("Bottom")={1};
-Physical Line("Right")={2};
-Physical Line("Top")={3};
-Physical Line("Left")={4};
-Physical Line("Cicrcle")={10,11,12,13};
-Physical Surface("Surface")={1};
-Physical Point("Center") = {11};
-
-Mesh.CharacteristicLengthExtendFromBoundary=1;
-Mesh.CharacteristicLengthFromPoints=1;
diff --git a/benchmarks/elasticity/hole.lua b/benchmarks/elasticity/hole.lua
deleted file mode 100644
index db59983d63844f9bbd4add9b3cc168b192f0cd48..0000000000000000000000000000000000000000
--- a/benchmarks/elasticity/hole.lua
+++ /dev/null
@@ -1,121 +0,0 @@
-a = 0.05
-T = 10e9
-E = 200e9
-nu = 0.3
--- --[[
---pull on sides
-function analyticalSolution(x,y,z)
- local r = math.sqrt(x*x+y*y)
- local theta = math.atan2(y,x)
- local ar = a/r
- local ur = T*r/(2*E) * (1-nu+(1+nu)*ar^2 + math.cos(2*theta)*((1+nu)*(1-ar^4)+4*ar^2))
- local utheta = - T*r/(2*E) * math.sin(2*theta)*((1+nu)*(1+ar^4)+2*(1-nu)*ar^2)
- return {ur*math.cos(theta)-utheta*math.sin(theta), ur*math.sin(theta)+utheta*math.cos(theta),0}
-end
--- --]]
-
- --[[
---pull on corner
-function analyticalSolution(x,y,z)
- local r = math.sqrt(x*x+y*y)
- local theta = math.atan2(y,x)
- local ar = a/r
- local ur = math.sin(2*theta)*T*r/E * (1+nu+4*ar^2-ar^4)
- local utheta = math.cos(2*theta)*T*r/E * (1+nu + 2*(1-nu)*ar^2+ar^4)
- return {ur*math.cos(theta)-utheta*math.sin(theta), ur*math.sin(theta)+utheta*math.cos(theta),0}
-end
- --]]
-
-function analyticalSolutionX (x,y,z)
- return analyticalSolution(x,y,z)[1]
-end
-
-function analyticalSolutionY (x,y,z)
- return analyticalSolution(x,y,z)[2]
-end
-
-
-function errorOnView (view, dim, solutionFunction)
- local basis = polynomialBasis.find(2)
- local gaussPoints = fullMatrix(1,1)
- local gaussWeights = fullMatrix(1,1)
- gaussIntegration.getTriangle(3,gaussPoints,gaussWeights)
-
- local nnodes = 3 -- only triangles
- local viewData = view:getData()
- local values = fullMatrix(nnodes,dim)
- local nodes = fullMatrix(nnodes,3)
- local gaussXYZ = fullMatrix(gaussPoints:size1(),3)
- local gaussValues = fullMatrix(gaussPoints:size1(),dim)
- local nEntities = viewData:getNumEntities(-1)
- local psi = fullMatrix(gaussPoints:size1(),3)
-
- local face = m:getFaceByTag(1)
-
- local errorViewFactory = PViewFactory("error", "ElementNodeData",m, 0, 1)
-
- basis:f(gaussPoints,psi)
- local nGaussPoints = gaussPoints:size1()
- local totalError = 0
- local errMat = fullMatrix(3,1)
- for i=0,nEntities-1 do
- local nElements = viewData:getNumElements(-1,i)
- for j=0,nElements-1 do
- local element = face:getMeshElement(j)
- if(viewData:getNumNodes(0,i,j)==nnodes) then
- viewData:getAllValuesForElement(0,i,j,values)
- viewData:getAllNodesForElement(0,i,j,nodes)
- gaussXYZ:gemm(psi,nodes,1,0)
- gaussValues:gemm(psi,values,1,0)
- for k=0,nGaussPoints-1 do
- local sol = solutionFunction(gaussXYZ:get(k,0),gaussXYZ:get(k,1),gaussXYZ:get(k,2))
- local err = 0
- for l=0,dim-1 do
- err = err + (sol[l+1]-gaussValues:get(k,l))^2
- end
- totalError = totalError + err*gaussWeights:get(k,0)*element:getJacobianDeterminant(gaussPoints:get(k,0),gaussPoints:get(k,1),gaussPoints:get(k,2))
- end
- local numVertices = element:getNumVertices()
- for k=0,numVertices-1 do
- local vertex = element:getVertex(k)
- local sol = solutionFunction(vertex:x(), vertex:y(), vertex:z())
- local err = 0
- for l=0,dim-1 do
- err = err + (sol[l+1]-values:get(k,l))^2
- end
- errMat:set(k,0,err)
- end
- errorViewFactory:setEntry(element:getNum(),errMat)
- end
- end
- end
- errorViewFactory:createView()
- return math.sqrt(totalError)
-end
-
-m = GModel()
-m:load("hole.geo")
-m:load("hole.msh")
-e = elasticitySolver(m,1)
-e:addElasticDomain (6, E, nu)
-
---e:addDirichletBCLua(1,5,0,"analyticalSolutionX")
---e:addDirichletBCLua(1,5,1,"analyticalSolutionY")
-e:addNeumannBC (1,5,{0,0,0})
-e:addDirichletBCLua(1,1,0,"analyticalSolutionX")
-e:addDirichletBCLua(1,1,1,"analyticalSolutionY")
-e:addDirichletBCLua(1,2,0,"analyticalSolutionX")
-e:addDirichletBCLua(1,2,1,"analyticalSolutionY")
-e:addDirichletBCLua(1,3,0,"analyticalSolutionX")
-e:addDirichletBCLua(1,3,1,"analyticalSolutionY")
-e:addDirichletBCLua(1,4,0,"analyticalSolutionX")
-e:addDirichletBCLua(1,4,1,"analyticalSolutionY")
-
-sys = linearSystemCSRTaucs()
-e:assemble(sys)
-sys:systemSolve()
-e:postSolve()
---view = e:buildVonMisesView("vonMises")
-view = e:buildDisplacementView("displacement")
-
-print(errorOnView(view,3,analyticalSolution))