Operator on an integral

1

Hi all, I want to compute the maximum von Mises stress over all the elements and compute the derivative in the future.

In order to use the automatic differentiation, I want to express the maximum von Mises stress with an integration via p-norm (not sure if we have any better choices?) as ( \int{ v^p(\mathbf{x}) \text{d}\mathbf{x}} )^{1/p} where v is the von-Mises stress. Here is the minimal code.

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.cpp.mesh import CellType
import ufl

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([4, 2, 0])], \
    [4,2], CellType.quadrilateral)
V = dolfinx.VectorFunctionSpace(mesh, ('CG', 1))

def clamped_boundary(x):
    return np.isclose(x[0], 0)
fdim = mesh.topology.dim - 1
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, clamped_boundary)
u_D = dolfinx.Function(V)
with u_D.vector.localForm() as loc:
    loc.set(0)
bc = dolfinx.DirichletBC(u_D, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets))

T = dolfinx.Constant(mesh, (0, 1)) # traction
f = dolfinx.Constant(mesh, (0, 0)) # body force

boundaries = [(1, lambda x: np.isclose(x[0], 4))]

# loop through all the boundary conditions to identify the facets
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_tag)

lambda_ = 1 # first Lame constant
mu = 1 # second Lame constant (shear modulus)
epsilon = lambda u: ufl.sym(ufl.grad(u))
sigma = lambda u: lambda_*ufl.nabla_div(u)*ufl.Identity(len(u)) + 2*mu*epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(1)

problem = dolfinx.fem.LinearProblem(a, L, bcs=[bc], petsc_options={'ksp_type': 'preonly', 'pc_type': 'lu'})
uh = problem.solve()

s = lambda u: sigma(u) - 1/3*ufl.tr(sigma(u))*ufl.Identity(len(u)) # deviatoric stress
vonStress = lambda u: ufl.sqrt(3/2*ufl.inner(s(u), s(u))) # von Mises stress

# we want to construct this integration with p-norm:
von_stress = vonStress(u)
p = 12.0
integral = (von_stress**p * ufl.dx) ** (1/p)

However, an error message appears:

Traceback (most recent call last):
  File "/shared/test.py", line 56, in <module>
    integral = (von_stress**p * ufl.dx) ** (1/p)
TypeError: unsupported operand type(s) for ** or pow(): 'Form' and 'float'

Although we can compute the integrand first and then consider the power later, it might induce the overflow problem. So do we have any remedies? Thanks!

2

You need to assemble the integral first, as you are creating a scalar integral, i.e.
(\int_{Omega} stress^p~\mathrm{d}x)^{1/p}=(\sum_{cell} \int_{cell} stress^p \mathrm{d}x)^{1/p}.
There is no way to simplify this root without assembling over all cells first.
However, you could normalize the stresses or the functional in general to avoid overflow.