Why FEniCS values are so different than theoretical ones?

Legacy DOLFIN
3

Thanks for the useful suggestion. It did work increasing mesh density, for 30, the FEM values increase a bit as follows -

Maximum pressure FEM:  1.04893 Hertz:  1.39916
Applied force    FEM:  0.28177 Hertz:  0.02930

But it didn’t get close to the calculated values on further increment.

I also changed the newton solver by adding the Custom Newton Solver by @nate and tried changing different linear solvers and preconditioners but it didn’t help much and the values were same as above.

Any further pointers are highly encouraged.

Here is the custom newton solver integrated MWE -

from dolfin import *
from mshr import *
import numpy as np

domain = Rectangle(Point(-1, -1),Point(1, 1))
mesh = generate_mesh(domain, 30)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1) and on_boundary
        
def bottom(x, on_boundary):
        return near(x[1], -1) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d+(pow(x[0],2))/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom)]

class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)

class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("pc_type", "ilu")
        self.linear_solver().set_from_options()

E = Constant(10.)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u)

problem = Problem(J, form, bc)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

p.assign(-project(sigma(u)[1, 1], V0))
gap.assign(project(obstacle-u[1], V2))

a = sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*pi*a**2)
print("Maximum pressure FEM: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FEM: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))