Hello everyone,
I am currently trying to transfer the " Hertzian contact with a rigid indenter using a penalty approach" example (Comet-FEniCS Example) into FEniCSx but I have difficulties formulating the penalty-formulation.
So far the code looks like this and fails in the last line:
import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl
N=30
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N//2)
def bottom(x):
return np.isclose(x[1],0)
def top(x):
return np.isclose(x[1], N//2)
fdim = domain.topology.dim -1
bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,fdim,top)
marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
metadata = {"quadrature_degree":4}
ds = ufl.Measure('ds',domain=domain,subdomain_data=facet_tag.find(2),metadata=metadata)
V = fem.VectorFunctionSpace(domain, ("CG",1))
u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
tdim = domain.topology.dim
fdim = tdim - 1
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)
R = 0.5
d = 0.02
x = ufl.SpatialCoordinate(domain)
circle = fem.Function(V)
circle.interpolate(lambda x: np.stack((x[0], x[1], -d+(pow(x[0],2)+pow(x[1], 2))/2/R)))
E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return lmbda*ufl.tr(eps(v))*ufl.Identity(3) + 2.0*mu*eps(v)
def ppos(x):
return (x+np.abs(x))/2.
pen = fem.Constant(domain, ScalarType(1e4))
form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen*ufl.dot(u_[2], ppos(u[2]-circle))*ds
The error I get is the following:
ufl.log.UFLException: Can't add expressions with different shapes.
Since it is no longer possible to use the Expression class, I interpolated an anonymous in 3D. I am sure there is a proper way of implementing this expression but I have not found something comparable.
This Post comes quite close to my problem, but there was no answer to:
“Are there particular ways that penalty methods have to be set up in FEniCSx if the level of penetration is a function of absolute position as opposed to the displacement (as is done in the Hertzian contact example)?”
I use FEniCSx 0.6.0 inside a Docker container.
Many thanks in advance!