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!