I am working on a test script (FEniCSx-real in docker) that combines elastodynamics (modified from Elastodynamics Example ) with a penalty method to prevent penetration of a rigid surface. In the test (the MWE included below), an elastic cube, starting at rest, is in free fall in a gravitational field over a rigid surface. The elastodynamics equations are iterated to solve the displacement of the cube, and the penalty method used at each step to prevent penetration of the rigid surface using the absolute position of the cube (displacement + reference configuration).
However, right now, if there is no initial penetration, then the object just remains in freefall and the center of mass (My in the code) follows a parabola. If there is some initial penetration, then the object follows an upward parabola (constant upward force). From this, it seems like the penalty term is staying at its initial value and not updating in the solver. If I print out values during the iterations, the values are changing.
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)?
The MWE of the elastodynamics/contact problem is included below:
import numpy as np
import ufl
import pyvista
import matplotlib.pyplot as plt
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx import *
#### Defining domain and function spaces
domain = mesh.create_rectangle(MPI.COMM_WORLD,[[0.0,0.0],[1,1]],[100,10],mesh.CellType.triangle)
V = fem.VectorFunctionSpace(domain,("CG",2)) # function space for the displacements and test functions
Vsig = fem.TensorFunctionSpace(domain,("DG",0)) # function space for the stress tensor
# gravity vector
b = -0.1
B = fem.Constant(domain,PETSc.ScalarType((0,b)))
#### Getting facets of the bottom edge that will come in contact ####
def bottom(x):
return np.isclose(x[1],0)
fdim = domain.topology.dim -1
bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
marked_facets = np.hstack([bottom_facets])
marked_values = np.hstack([np.full(len(bottom_facets),3,dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain,fdim,marked_facets[sorted_facets],marked_values[sorted_facets])
#### Defining elasticity constants ####
E = 1000.0
nu = 0.3
mu = fem.Constant(domain,E / (2.0*(1+nu)))
lmbda = fem.Constant(domain,E*nu / ((1.0+nu)*(1.0-2.0*nu)))
k_pen = fem.Constant(domain,100000.0)
rho = 1.0
#### Defining parameters associated with Generalized alpha method and timestepping ####
am = 0.2
af = 0.4
g = 0.5 + af - am
alpha_m = fem.Constant(domain,am)
alpha_f = fem.Constant(domain,af)
gamma = fem.Constant(domain,g)
beta = fem.Constant(domain,(g+0.5)**2 / 4.0)
T = 10
Nsteps = 100
dt = fem.Constant(domain,T/Nsteps)
#### Defining functions of interest ####
u_ = ufl.TestFunction(V)
u = fem.Function(V,name='Displacement')
u_old = fem.Function(V)
v_old = fem.Function(V)
a_old = fem.Function(V)
penetrate = fem.Function(V) # value of the penetration at each element
abs_penetrate = fem.Function(V)
X0 = fem.Function(V) # Reference Configuration of the domain
X0.interpolate(lambda x: [x[0],x[1]])
X = fem.Function(V) # Current Configuration of the domain
X.interpolate(X0)
NN = len(X.x.array)
# array indices of the x- and y-components of the vector-valued function arrays
x_inds = np.arange(0,NN,2)
y_inds = np.arange(1,NN+1,2)
#### Defining measures ####
metadata = {"quadrature_degree":4}
ds = ufl.Measure('ds',domain=domain,subdomain_data=facet_tag,metadata=metadata)
dx = ufl.Measure("dx",domain=domain,metadata=metadata)
#### defining elements of GA method ####
def sigma(r):
return lmbda*ufl.nabla_div(r)*ufl.Identity(len(r)) + 2*mu*ufl.sym(ufl.grad(r))
def m(u,u_):
return rho*ufl.inner(u,u_)*dx
def k(u,u_):
return ufl.inner(sigma(u),ufl.grad(u_))*dx
def penalty(u,u_):
X.x.array[:] = u.x.array[:] + X0.x.array[:]
penetrate.x.array[y_inds] = 0.0 - X.x.array[y_inds] # Boundary at y=0
abs_penetrate.x.array[:] = (penetrate.x.array[:] + np.abs(penetrate.x.array[:]))/2
return ufl.inner(k_pen*abs_penetrate,u_)*ds(3)
def update_a(u,u_old,v_old,a_old,ufl_=True):
if ufl_:
dt_ = dt
beta_ = beta
else:
dt_ = dt
beta_ = beta.value
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
def update_v(a,u_old,v_old,a_old,ufl_=True):
if ufl_:
dt_ = dt
gamma_ = gamma
else:
dt_ = dt
gamma_ = gamma.value
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(u,u_old,v_old,a_old):
u_vec,u0_vec = u.x.array[:],u_old.x.array[:]
v0_vec,a0_vec = v_old.x.array[:],a_old.x.array[:]
a_vec = update_a(u_vec,u0_vec,v0_vec,a0_vec,ufl_=False)
v_vec = update_v(a_vec,u0_vec,v0_vec,a0_vec,ufl_=False)
v_old.x.array[:] = v_vec
a_old.x.array[:] = a_vec
u_old.x.array[:] = u_vec
def avg(x_old,x_new,alpha):
return alpha*x_old + (1-alpha)*x_new
a_new = update_a(u,u_old,v_old,a_old,ufl_=True)
v_new = update_v(a_new,u_old,v_old,a_old,ufl_=True)
#### Defining the problem and solver ####
res = m(avg(a_old,a_new,alpha_m),u_) + k(avg(u_old,u,alpha_f),u_) - rho*ufl.inner(B,u_)*dx - penalty(u,u_) #residual to solve
problem = fem.petsc.NonlinearProblem(res,u)#,bcs=bcs)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
#### Quatities to track ####
time = np.linspace(0, T, Nsteps+1)
V0 = 1
# center of mass positions
Mx_ = np.zeros((Nsteps+1,))
My_ = np.zeros((Nsteps+1,))
Mx = fem.form((1/V0)*(u[0]+X0[0])*dx)
My = fem.form((1/V0)*(u[1]+X0[1])*dx)
Mx_[0] = fem.assemble_scalar(Mx)
My_[0] = fem.assemble_scalar(My)
#### Iterating through time ####
for (i, dt) in enumerate(np.diff(time)):
t = time[i+1]
num_its,converged = solver.solve(u) # solve the current time step
assert(converged)
u.x.scatter_forward()
# Update old fields with new quantities
update_fields(u, u_old, v_old, a_old)
Mx_[i+1] = fem.assemble_scalar(Mx)
My_[i+1] = fem.assemble_scalar(My)
# Plot center of mass evolution
plt.figure()
plt.plot(time,Mx_,time,My_,time,My_[0]+0.5*b*time**2)
plt.legend(("Mx","My","My_freefall"))
plt.xlabel("Time")
plt.ylabel("Center of Mass")
plt.show()
Thanks so much!