diff --git a/benchmarks/elasticity/hole.geo b/benchmarks/elasticity/hole.geo
new file mode 100644
index 0000000000000000000000000000000000000000..9187cf9f3d8dfe1374a78c50649bc6334deff7f8
--- /dev/null
+++ b/benchmarks/elasticity/hole.geo
@@ -0,0 +1,39 @@
+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
new file mode 100644
index 0000000000000000000000000000000000000000..db59983d63844f9bbd4add9b3cc168b192f0cd48
--- /dev/null
+++ b/benchmarks/elasticity/hole.lua
@@ -0,0 +1,121 @@
+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))