Benchmarking elasticity example against Comsol solution provides different results

Question

I'm modeling square of size 1 with small hole of radius r = 0.05 in the center of it. The left side of the square is fixed (zero displacement) while the right side of the square is being pulled.

from __future__ import print_function
from dolfin import *
from dolfin_adjoint import *

n = 100
mesh  = UnitSquareMesh(n, n)
p = Constant(2)        
eps = Constant(1.0e-8) 

class HoleExpression(Expression):
    def eval(self, values, x):
        values[0] = 1.0
        if (x[0]-0.5)**2 + (x[1]-0.5)**2 <  (0.05)**2:
            values[0] = 0
        return values

left  =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]

A = FunctionSpace(mesh, "CG", 3)  
PP = VectorFunctionSpace(mesh, "CG", 3)

g = interpolate(Expression(('1.0','0.0'), degree = 2), PP) 

bc1 = DirichletBC(PP, Constant((0.0, 0.0)), left)
bcs = [bc1]

parameters["std_out_all_processes"] = False

def k(a):
    return eps + (1 - eps) * a**p

E, nu = 1.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 - nu**2)))


def sigma(rho, v):
    return 2.0*mu*rho*sym(grad(v)) + lmbda*rho*tr(sym(grad(v)))*Identity(len(v))

def get_u(a):
    u = TrialFunction(PP)
    v = TestFunction(PP)
    frm = inner(sigma(k(a), u), grad(v))*dx
    rhs = dot(g, v)*ds(2)
    u = Function(PP)
    problem = LinearVariationalProblem(frm, rhs, u, bcs = bcs)
    solver = LinearVariationalSolver(problem)
    solver.parameters["linear_solver"] ="lu"
    solver.solve()
    return u
if __name__ == "__main__":
    a = interpolate(HoleExpression(degree = 2), A, name="Control")
    u = get_u(a)
    File("u.pvd") << u

From comsol I'm getting the following result for displacement among x:

and among y:

The scale of values on both pictures are in range 1e-5.

While the code above provide a different solution:
a displacement among x:

and a displacement among y:

The pictures are remotely similar, but what bothers me is the fact, that for material with E=1, mu=0.3 I'm getting displacement in range of 1e0, while comsol provides me with displacements in range 1e-5. (The comsol file could be obtained from https://yadi.sk/d/ZIFOPoZ73JZdu9 or used original https://www.comsol.com/community/exchange/56/ (but it needs modifications in material parameters)).

Comments

I have never used comsol but linear elasticity is, of course, a linear problem, and so it scales linearly. If the solutions are so similar, maybe your problem is in the units. Please check that E, mu and F (applied force) are being set within the same scale. If not, no idea.

---

You have a 1x1 mesh, with E=1, and thickness = 1, plane stress, unit load per lenght, right? If yes, then the stiffness, the load and the displacement will be in the order of 1. Thus, Dolfin's solution is right, and Comsol solution is wrong. Or, more probably, you used different data, as already suggested by @nabarnaf.
That said, you should avoid doing dirty trick in Dolfin to define the hole. Simple use a mesh with an hole, as you are doing in Comsol.

---

Yes, I think it was a scale issue, because I was able to replicate stress field with good accuracy.

Answer 1

Print the solution at nodes don't blindly believe the legend in FEniCS.

To print solution at a node point. use:

point = (1,1)
print 'numerical u at the corner point:', u(point)