Enforcing Average Prescribed Displacement

General
1

Hello,

I am brand new to FEniCSx and I am trying to solve a problem where I specify the displacement u_D in average sense on a boundary Gamma_D. The weak form for this problem looks like this:

image
image1066×268 26.4 KB

where Psi_,u is the first variation of some hyperelastic material model. The displacement is enforced in an average sense, rather than pointwise, by letting q be discretized by a scalar (if u_D is enforced in a single direction) or a vector (if u_D is enforced in more than one direction).
At present, I am trying to enforce u_D only in one direction, such that q can be discretized by a scalar, i.e. it is essentially just a single Lagrange multiplier.

I found this tutorial on hyperelasticity:

tutorial link

And tried to change the enforced traction to an enforced average displacement, but I have not been able to do so. Do you have suggestions on how to add q as a variable to the system?

Best,
Andreas

2

Your Lagrange multiplier is a single scalar, so you’d need the ‘real_space’ from scifem. See: scifem/tests/test_real_functionspace.py at main · scientificcomputing/scifem · GitHub

The snippet test_singular_poisson also illustrates how you would use that inside a weak form. You’d have to replace the lines:

a01 = ufl.inner(c, v) * dx
a10 = ufl.inner(u, d) * dx

with the two boundary integrals that involve the Lagrange multiplier.

2 Likes
3

Thanks, Stein

I am closer to having it work now. I added:

R = scifem.create_real_functionspace(domain, (1,))
disp = fem.Constant(domain, default_scalar_type((0)))
q = ufl.TrialFunction(R)
test_q = ufl.TestFunction(R)

to the example in the tutorial, and changed

F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

to

F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx + ufl.inner(test_q[0],(u[1]-disp)) * ds(2) + ufl.inner(q[0],v[1]) * ds(2)

But now I don’t know how to add the variable q to the system in this line:

problem = NonlinearProblem(F, u, bcs)

Rather than only ‘u’, it should be ‘u’ and ‘q’. Otherwise, you get:

ValueError: Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_0
  v_1
  v_0

As far as I understand, the way to go is making a mixed formulation, but I can’t get it to work with

R = scifem.create_real_functionspace(domain, (1,))

Is it the right approach, or is there something that I am missing?

4

See for instance

2 Likes
5

Thanks a lot @dokken now it works :smiley:

For the sake of completeness, here is the edited version of the tutorial on Hyperelasticty with enforced average displacement condition.

I only started using FEniCS and FEniCSx yesterday, so this code is very open to improvements.

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
#from dolfinx.nls.petsc import NewtonSolver
from scifem import NewtonSolver, assemble_scalar, BlockedNewtonSolver
import pyvista
import numpy as np
import ufl
import scifem

from mpi4py import MPI
from dolfinx import fem, mesh, plot

from petsc4py import PETSc
ftype = PETSc.ScalarType

NX, NY = 20, 10         # number of elements in horizontal direction and vertical direction
LX, LY = 10., 5.       # length in repective directions

# Create a 2D rectangular mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [LX, LY]], [NX, NY], mesh.CellType.quadrilateral)

# Define the function space for 2D
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))
R = scifem.create_real_functionspace(domain, (1,)) # Adding a single Lagrange multiplier to the global system

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[1], LY) * (x[0]>0.8*LX-1e-5) # The boundary where the displacement is prescribed

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# Fixed boundary condition on the left edge
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

# Enforved displacement
disp = fem.Constant(domain, default_scalar_type((0)))

# Define function space including the lagrange multiplier
W = ufl.MixedFunctionSpace(V, R)
v, test_q = ufl.TestFunctions(W)
du, dq = ufl.TrialFunctions(W)

# Set zero-valued initial guess
u = fem.Function(V, dtype=ftype)
#u.x.array[:] = 0.0 # edit this for non-zero initial guess
q = fem.Function(R, dtype=ftype)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0)
nu = default_scalar_type(0.4)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F0 = ufl.inner(ufl.grad(v), P) * dx + ufl.inner(q[0],v[1]) * ds(2)
F1 = ufl.inner(test_q[0],(u[1]-disp)) * ds(2) 
F_ = F0 + F1
F = list(ufl.extract_blocks(F_))
J = ufl.extract_blocks(ufl.derivative(F_, u, du) + ufl.derivative(F_, q, dq))

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
solver = BlockedNewtonSolver(F, [u, q], bcs=bcs, J=None, petsc_options=petsc_options)
#solver = NewtonSolver(F, J, [u, q], bcs=bcs, max_iterations=25, petsc_options=petsc_options)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

pyvista.start_xvfb()
plotter = pyvista.Plotter()
plotter.open_gif("deformation.gif", fps=3)

topology, cells, geometry = plot.vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)

values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid["u"] = values
function_grid.set_active_vectors("u")

# Warp mesh by deformation
warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")

# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])

# Compute magnitude of displacement to visualize in GIF
Vs = fem.functionspace(domain, ("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(ufl.sqrt(sum([u[i]**2 for i in range(len(u))])), Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array

log.set_log_level(log.LogLevel.INFO)
tval0 = -1.0
for n in range(1, 10):
    disp.value = n * tval0
    num_its = solver.solve()
    u.x.scatter_forward()
    q.x.scatter_forward()
    print(f"Time step {n}, Number of iterations {num_its}, Disp {disp.value}, Load {q.x.array[:]}")
    function_grid["u"][:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
    magnitude.interpolate(us)
    warped.set_active_scalars("mag")
    warped_n = function_grid.warp_by_vector(factor=1)
    warped.points[:, :] = warped_n.points
    warped.point_data["mag"][:] = magnitude.x.array
    plotter.update_scalar_bar_range([0, 1])
    plotter.write_frame()
plotter.close()