Neumann bcs in mixed space

dolfinx
1

Hi,

how can I define a Neumann bc for a subdomain in mixed space? Before, I defined the bcs differently, but now I would like to use the meshtags. I guess there is a problem in Measure , but I cant figure out how to do it correctly.

My MWE is :

from dolfinx import  default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from basix.ufl import mixed_element

# Define mesh
L = 1.0
p = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)
V_el = mixed_element([Ue, Te]) # u, P

# Define function spaces
V = fem.functionspace(domain, V_el)
V0, _ = V.sub(0).collapse()

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

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

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
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])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V0, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V0)]

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Define functions
w = fem.Function(V)
u, P = ufl.split(w)
δw = ufl.TestFunction(V)
δu, δP = ufl.split(δw)

#Define physics
I = ufl.Identity(3)  
F = I + ufl.grad(u) 
E = 1 / 2 * (F.T * F - I)  
E_el = E - P  
S =  E_el + ufl.tr(E_el) * I  
δE = ufl.derivative( ufl.grad(u) , w, δw)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
form = ( 
        ufl.inner(δE, S) * dx
    - ufl.inner(δu, T) * ds(2)
    )

problem = NonlinearProblem(form, w, bcs)

and the error message is:

Traceback (most recent call last):
  File "/home/bjoern/examples/plasticity/mwe.py", line 66, in <module>
    problem = NonlinearProblem(form, w, bcs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 922, in __init__
    self._L = _create_form(
              ^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 337, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 331, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 254, in _form
    ufcx_form, module, code = jit.ffcx_jit(
                              ^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 62, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 212, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms
    raise e
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
    ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/ir/representation.py", line 135, in compute_ir
    _compute_integral_ir(
  File "/usr/lib/python3/dist-packages/ffcx/ir/representation.py", line 396, in _compute_integral_ir
    integral_ir = compute_integral_ir(
                  ^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/ir/integral.py", line 91, in compute_integral_ir
    mt_table_reference = build_optimized_tables(
                         ^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/ir/elementtables.py", line 444, in build_optimized_tables
    t = get_ffcx_table_values(
        ^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/ir/elementtables.py", line 147, in get_ffcx_table_values
    entity_points = map_integral_points(points, integral_type, cell, entity)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffcx/ir/representationutils.py", line 118, in map_integral_points
    assert points.shape[1] == tdim - 1
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError

It’s Fenicx 0.9 on WSL Ubuntu / Windows. Many thanks!

2

The issue doesn’t seem related to the facet_tags: the same error message is printed with ds = ufl.Measure('ds', domain=domain)

Rather, something awkward is triggered by the the quadrature_element, which should not even affect the ds part of the form. If you replace the quadrature element He with He = basix.ufl.element("CG", domain.basix_cell(), p) in your MWE, then the script runs…

3

Thank you! Now I can build NonlinearProblem. May I still ask if I defined the bc correctly? It is supposed to just forcing one side of the cube. I am asking, since I still have convergence problems for the full problem. I believe they stem from the bc’s, because if I impose a displacement instead of Neumann bc, it works (I actually wonder a bit about that, because I thought I would need the quadrature elements for plasticity).

Just for clarity, I post the full problem:

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from dolfinx.fem import locate_dofs_geometrical
from basix.ufl import mixed_element

name = "plasticity_beam"
o, p = 1, 1

# Define mesh
L = 1.0
p = 2
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.element("CG", domain.basix_cell(), p)
# He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)
V_el = mixed_element([Ue, Te, He, Te]) # u, P, h, B
V = fem.functionspace(domain, V_el)

# Define function spaces
Uf = fem.functionspace(domain, Ue)
Tf = fem.functionspace(domain, Te)
Hf = fem.functionspace(domain, He)

V0, _ = V.sub(0).collapse()

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

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

dofs_sol = locate_dofs_geometrical(V0, lambda x: np.isclose(x[0], L) & np.isclose(x[1], 1)& np.isclose(x[2], 1))

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

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])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)

left_dofs = fem.locate_dofs_topological(V0, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V0)]

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Define functions
w = fem.Function(V)
w0 = fem.Function(V)
u, P, h, B = ufl.split(w)
u0, P0, h0, B0 = ufl.split(w0)
S0 = fem.Function(Tf, name="S0")

δw = ufl.TestFunction(V)
δu, δP, δh, δB = ufl.split(δw)

#Define physics
mu = fem.Constant(domain, 100.0)
la = fem.Constant(domain, 10.)
Sy = fem.Constant(domain, 0.300)  
bh = fem.Constant(domain, 20.00)  
qh = fem.Constant(domain, 0.100) 
bb = fem.Constant(domain, 250.0) 
qb = fem.Constant(domain, 0.100) 

def rJ2(A):
    """Square root of J2 invariant of tensor A"""
    J2 = 1 / 2 * ufl.inner(A, A)
    rJ2 = ufl.sqrt(J2)
    return ufl.conditional(rJ2 < 1.0e-12, 0.0, rJ2)

I = ufl.Identity(3)  
F = I + ufl.grad(u)  
E = 1 / 2 * (F.T * F - I)  
E_el = E - P  
S = 2 * mu * E_el + la * ufl.tr(E_el) * I
S, B, h = ufl.variable(S), ufl.variable(B), ufl.variable(h)
f = ufl.sqrt(3) * rJ2(ufl.dev(S - B)) - (Sy + h)  
g = f
dgdS = ufl.diff(g, S)
S, B, h = S.expression(), B.expression(), h.expression()
δE = ufl.derivative(E, w, δw)
dλ = ufl.max_value(f, 0) 

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
form = ( 
        ufl.inner(δE, S) * dx
    + ufl.inner(δP, (P - P0) - dλ * dgdS) * dx
    + ufl.inner(δh, (h - h0) - dλ * bh * (qh * 1.00 - h)) * dx
    + ufl.inner(δB, (B - B0) - dλ * bb * (qb * dgdS - B)) * dx
    - ufl.inner(δu, T) * ds(2)
    )

problem = NonlinearProblem(form, w, bcs)

solver = NewtonSolver(domain.comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

log.set_log_level(log.LogLevel.INFO)
tval0 = 0.1
for n in range(1, 2):
    T.value[0] = n * tval0
    T.value[1] = n * tval0
    T.value[2] = n * tval0

    num_its, converged = solver.solve(w)
    assert (converged)
    w.x.scatter_forward()
    print(n)
4

This could certainly be true, my ‘fix’ was just meant to illustrate that the bug was somewhere else in the code, and seems to be a dolfinx issue ( @dokken ).

It has been too long since I worked on something plasticity related, so I can’t comment on the correctness of the full form, but the Neumann boundary terms - ufl.inner(δu, T) * ds(2) certainly does look correct to me.

Just to be clear, currently T is a zero vector. Even then this issue occurs? (Of course, if T remains a zero vector, then you don’t need this integral…)

1 Like
5

Including a quadrature element in a ds integral does not make sense, if the quadrature element is defined in the cell itself.

A Quadrature element is a set of values at the quadrature points (in physical space). There are no underlying basis functions, that makes it possible to tabulate them at any other point within the physical/reference cell.

6

Right, but unless I am overlooking something, the performed integral does not involve functions living in the quadrature space. So that seems like a bug to me?

7

Ah, so the issue lies in:

which returns true for a mixed element.

It is not so clear how one should circumvent this (without using ufl.MixedFunctionSpace for each of the spaces, and get a block form).

I.e. do something like:

from dolfinx import  default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
import basix
from basix.ufl import mixed_element

# Define mesh
L = 1.0
p = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

Ue = basix.ufl.element("P", domain.basix_cell(), p, shape=(domain.geometry.dim,))
He = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), degree=p)
Te = basix.ufl.blocked_element(He, shape=(domain.geometry.dim, domain.geometry.dim), symmetry=True)

Us = fem.functionspace(domain, Ue)
Te = fem.functionspace(domain, Te)

# Define function spaces
V = ufl.MixedFunctionSpace(Us, Te)

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

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

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
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])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(Us, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, Us)]

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Define functions
u = fem.Function(Us)
P = fem.Function(Te)
δu, δP = δw = ufl.TestFunctions(V)

#Define physics
I = ufl.Identity(3)  
F = I + ufl.grad(u) 
E = 1 / 2 * (F.T * F - I)  
E_el = E - P  
S =  E_el + ufl.tr(E_el) * I  
δE = ufl.derivative( ufl.grad(u) , [u, P], δw)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
form = ( 
        ufl.inner(δE, S) * dx
    - ufl.inner(δu, T) * ds(2)
    )

F_compiled = fem.form(ufl.extract_blocks(form))
J_compiled = fem.form(ufl.extract_blocks(ufl.derivative(form, [u, P], ufl.TrialFunctions(V))))

and then for instance use the BlockedNewtonSolver from scifem:

or soon the PETSc SNES implementation in DOLFINx, ref: SNES solver interface for standard, block and nest systems by jorgensd · Pull Request #3648 · FEniCS/dolfinx · GitHub

2 Likes
8

@Stein @dokken many thanks for the explanation!

9

Just to be clear, currently T is a zero vector. Even then this issue occurs? (Of course, if T remains a zero vector, then you don’t need this integral…)

I am updating T in the loop with T.value[0] = n * tval0. I will try to solve it now with the suggested steps using the BlockedNewtonSolver. Thanks again for the help!