MixedElement formulation leads to unintuitive behaviour of the API

dolfinx
1

Hello all,

We are currently trying to solve a coupled PDE with 4 degrees of freedom (3 displacement and 1 pressure) preferably with a MixedElement approach.
Right now we have several issues with the code (the snippet omits the definition of the weak forms but it should be clear anyway):

  • after definition of the MixedElement we can’t access the underlying Functions for u and lambda (pressure). [Line 58 and following]
  • definition of the boundary conditions seems to be way to complex for such an easy problem setup (especially for all the single components of u). [Line 68 - 135]
  • the PDE contains the time derivative of the displacement. For the approximation we try to incorporate a simple backwards differentiation scheme based on the currently known solution of u. However, this approach does not work well with the MixedElement formulation (due to the fact that we cannot extract the underlying functions of the MixedElement) [Line 141 - 147]
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import mesh, fem, plot, io, nls

# ---------------------------------------------------------------------------- #
#                                  Parameters                                  #
# ---------------------------------------------------------------------------- #
L, W, H = 1, 1, 1
mesh_size = 0.2

# Simulation parameter
dt = 0.1

# Elasticity parameters
E = 1e7
nu = 0.3

# Environmental influences
b = [0, 0, 0]

# 3D box with hexahedron cells
n_elements = [max(int(x / mesh_size), 1) for x in [L, W, H]]
domain = mesh.create_box(
    comm=MPI.COMM_WORLD,
    points=[np.array([0, 0, 0]), np.array([L, W, H])],
    n=n_elements,
    cell_type=mesh.CellType.hexahedron,
)


# ---------------------------------------------------------------------------- #
#                                FiniteElements                                #
# ---------------------------------------------------------------------------- #
# u  = displacement (dim=3)
# lm = pressure (dim=1)

element_u = ufl.VectorElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=2,
)
element_lm = ufl.FiniteElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=1,
)

element = ufl.MixedElement(element_u, element_lm)
V = fem.FunctionSpace(domain, element)

# The trial/test functions get defined on the (mixed) function space
dudlm = ufl.TestFunction(V)
du, dlm = ufl.split(dudlm)

ulm = fem.Function(V)
(u, lm) = ufl.split(ulm)

print(type(lm))  # ufl.indexed.Indexed      <<<< why no Function
print(type(u))  # ufl.tensors.ListTensor   <<<< from split?

# we use dx as the measure on the domain volume and ds on the boundary
dx = ufl.dx
ds = ufl.ds


# ---------------------------------------------------------------------------- #
#                              Boundary Conditions                             #
# ---------------------------------------------------------------------------- #
# is there an easier way for the definition of the BCs in mixed elements?


def clamped_boundary(x, axis, offset):
    return np.isclose(x[axis], offset)


clamp_u = lambda x: clamped_boundary(x, 2, 0)
clamp_lm_left = lambda x: clamped_boundary(x, 0, 0)
clamp_lm_right = lambda x: clamped_boundary(x, 0, L)

facet_dim = domain.topology.dim - 1

# get the mesh entities, where the clamped boundary shall be applied
boundary_facets_u = mesh.locate_entities_boundary(
    domain,
    facet_dim,
    clamp_u,
)
boundary_facets_lm_left = mesh.locate_entities_boundary(
    domain,
    facet_dim,
    clamp_lm_left,
)
boundary_facets_lm_right = mesh.locate_entities_boundary(
    domain,
    facet_dim,
    clamp_lm_right,
)

# get the topological (element) dofs that correspond to the clamped boundary
topological_dofs_x = fem.locate_dofs_topological(
    V.sub(0).sub(0),
    facet_dim,
    boundary_facets_u,
)
topological_dofs_y = fem.locate_dofs_topological(
    V.sub(0).sub(1),
    facet_dim,
    boundary_facets_u,
)
topological_dofs_z = fem.locate_dofs_topological(
    V.sub(0).sub(2),
    facet_dim,
    boundary_facets_u,
)
topological_dofs_lm_left = fem.locate_dofs_topological(
    V.sub(1),
    facet_dim,
    boundary_facets_lm_left,
)
topological_dofs_lm_right = fem.locate_dofs_topological(
    V.sub(1),
    facet_dim,
    boundary_facets_lm_right,
)

u_D = PETSc.ScalarType(0.0)
bcx = fem.dirichletbc(u_D, topological_dofs_x, V.sub(0).sub(0))
bcy = fem.dirichletbc(u_D, topological_dofs_y, V.sub(0).sub(1))
bcz = fem.dirichletbc(u_D, topological_dofs_z, V.sub(0).sub(2))
bclml = fem.dirichletbc(PETSc.ScalarType(0), topological_dofs_lm_left, V.sub(1))
bclmr = fem.dirichletbc(PETSc.ScalarType(50), topological_dofs_lm_right, V.sub(1))

bcs = [bcx, bcy, bcz, bclml, bclmr]


# ---------------------------------------------------------------------------- #
#                  Definition of a time differentiation method                 #
# ---------------------------------------------------------------------------- #
# required: du/dt

u_prev = fem.Function(V.sub(0))
u_prev.interpolate(lambda x: np.zeros_like(x[:3]))  # initialize with 0

u_dot = (u - u_prev) / dt
# print(u_dot.x.array)  # <<< does not work. Is there an alternative?

# ---------------------------------------------------------------------------- #
#                               Define Weak Form                               #
# ---------------------------------------------------------------------------- #
F = get_weak_form()  # (dummy)

# ---------------------------------------------------------------------------- #
#                                     Solve                                    #
# ---------------------------------------------------------------------------- #

problem = fem.petsc.NonlinearProblem(F, ulm, bcs=bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
# ... some more options

for i in range(1):
    uh = solver.solve(ulm)
    u_prev.x.array[:] = u.x.array[:] # <<< does not work either (u is ListTensor)

Any input is highly appreciated.
Thanks in advance and all the best

2

If you consider the signature of ufl.split:

In [1]: ufl.split?
Signature: ufl.split(v)
Docstring:
UFL operator: If v is a Coefficient or Argument in a mixed space, returns
a tuple with the function components corresponding to the subelements.
File:      ~/shared/ufl/ufl/split_functions.py
Type:      function

it states that it returns the components of the function, not a DOLFINx function itself. To get a variable that is a dolfinx.fem.Function in the appropriate function space, you need to call ulm.split().
You can also see: fenics-project / DOLFIN / issues / #194 - split(u) versus u.split() — Bitbucket

You could of course use for-loops to create these boundary conditions instead of declaring them one by one.

For assigning data to sub-spaces, see: fenics-project / DOLFIN / issues / #194 - split(u) versus u.split() — Bitbucket

Please note that this does not work, as it will throw an error message:

RuntimeError                              Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 u_prev = fem.Function(V.sub(0))

File ~/shared/dolfinx/python/dolfinx/fem/function.py:246, in Function.__init__(self, V, x, name, dtype)
    244     self._cpp_object = functiontype(dtype)(V._cpp_object, x)
    245 else:
--> 246     self._cpp_object = functiontype(dtype)(V._cpp_object)
    248 # Initialize the ufl.FunctionSpace
    249 super().__init__(V.ufl_function_space())

RuntimeError: Cannot create Function from subspace. Consider collapsing the function space
3

Hello @dokken and thanks for your comments on our questions.

With the proposed solution ulm.split() and the extracted sub-functions u and lm the formulation of the weak form was not possible anymore and the whole Jupyter-Server crashed silently when trying to solve. Thus we still need to use ufl.split(ulm) for the formulation and ulm.split() for the extraction of subfunctions? This seems redundant. Is there a better way?

Furthermore: The actual reason we want the subfunctions (to get access to the values of the array of one variable is not solved with this approach. After the split, the subfunctions do not contain only the values of the respective variable but the values of the whole (mixed) space. The following example should explain what I mean:

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import mesh, fem, plot, io, nls

# 3D box with 1 hexahedron cell
domain = mesh.create_box(
    comm=MPI.COMM_WORLD,
    points=[np.array([0, 0, 0]), np.array([1, 1, 1])],
    n=[1, 1, 1],
    cell_type=mesh.CellType.hexahedron,
)


# ---------------------------------------------------------------------------- #
#                                FiniteElements                                #
# ---------------------------------------------------------------------------- #
# u  = displacement (dim=3)
# lm = pressure (dim=1)

element_u = ufl.VectorElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=1,
)
element_lm = ufl.FiniteElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=1,
)

element = ufl.MixedElement(element_u, element_lm)
V = fem.FunctionSpace(domain, element)

# The trial/test functions get defined on the (mixed) function space
dudlm = ufl.TestFunction(V)
du, dlm = ufl.split(dudlm)

ulm = fem.Function(V)

(u, lm) = ufl.split(ulm)
(u_, lm_) = ulm.split()

# we use dx as the measure on the domain volume and ds on the boundary
dx = ufl.dx
ds = ufl.ds

# ---------------------------------------------------------------------------- #
#                                  Unintuitive                                 #
# ---------------------------------------------------------------------------- #
print(ulm.x.array.shape) # >> (32,) = ( 8 x 3 ) + ( 8 x 1 )
print(u_.x.array.shape)  # >> (32,) (expected 24)
print(lm_.x.array.shape) # >> (32,) (expected 8)

print(ulm.sub(0).collapse().x.array.shape) # >> (24,)
print(ulm.sub(1).collapse().x.array.shape) # >> (8,)

Regarding assignment of data to subspaces: Could you please check, if the link you provided is the correct one? It points to the same page as the split-Issue

The part with collapsing is still quite confusing. Why would you want to remove the “connection” between the source subspace and the one used to define a new function. Shouldn’t they be always related? Same thing for the collapsing of subfunctions in the example above: If I collapse it and my original function updates (e.g. during solution), my collapsed function will not update?

Thanks in advance and all the best

4

I meant to link to: How to correctly assign values to mixed elements in FEniCS-X? - #3 by dokken
which gives you the following workflow:

from IPython import embed
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io

L, W, H = 1, 1, 1
mesh_size = 0.2

# 3D box with hexahedron cells
n_elements = [max(int(x / mesh_size), 1) for x in [L, W, H]]
domain = mesh.create_box(
    comm=MPI.COMM_WORLD,
    points=[np.array([0, 0, 0]), np.array([L, W, H])],
    n=n_elements,
    cell_type=mesh.CellType.hexahedron,
)

element_u = ufl.VectorElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=2,
)
element_lm = ufl.FiniteElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=1,
)

element = ufl.MixedElement(element_u, element_lm)
V = fem.FunctionSpace(domain, element)

# The trial/test functions get defined on the (mixed) function space
dudlm = ufl.TestFunction(V)
du, dlm = ufl.split(dudlm)

ulm = fem.Function(V)
(u, lm) = ulm.split()
# Fill ulm with dummy data
ulm.x.array[:] = np.arange(len(ulm.x.array))
# Create subspace and map from sub space to parent
V0, V0_to_V = V.sub(0).collapse()


u_prev = fem.Function(V0)
u_prev.x.array[:] = u.x.array[V0_to_V]
u_prev.name = "Collapsed space"
with io.XDMFFile(domain.comm, "u_prev.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u_prev)

with io.XDMFFile(domain.comm, "u0.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(ulm.sub(0))

In many cases, the users do not need the map from the sub space to the parent space (for instance when defining Dirichlet conditions, as we send in both the collapsed and non-collapsed space), and therefore we do not store is as part of the Function object (as it increases memory footprint).

When you want to write variational forms, you should use ufl.split(ulm) to get each component of the function.
If you want to access data from these functions as a dolfinx.fem.Function object, you should use ulm.split() or ulm.sub(i) where i is the index of the sub space.

Here is a minimal example illustrating the use-cases:

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, nls, log

L, W, H = 1, 1, 1
mesh_size = 0.2

# 3D box with hexahedron cells
n_elements = [max(int(x / mesh_size), 1) for x in [L, W, H]]
domain = mesh.create_box(
    comm=MPI.COMM_WORLD,
    points=[np.array([0, 0, 0]), np.array([L, W, H])],
    n=n_elements,
    cell_type=mesh.CellType.hexahedron,
)

element_u = ufl.VectorElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=2,
)
element_lm = ufl.FiniteElement(
    family="Lagrange",
    cell=domain.ufl_cell(),
    degree=1,
)

element = ufl.MixedElement(element_u, element_lm)
V = fem.FunctionSpace(domain, element)

# The trial/test functions get defined on the (mixed) function space
dudlm = ufl.TestFunction(V)
du, dlm = ufl.split(dudlm)

ulm = fem.Function(V)
(u, lm) = ufl.split(ulm)

x = ufl.SpatialCoordinate(domain)
F = ufl.inner(u, du)*ufl.dx + ufl.inner(lm, dlm)*ufl.dx - \
    ufl.inner(x, du)*ufl.dx + ufl.inner(x[2], dlm)*ufl.dx

log.set_log_level(log.LogLevel.INFO)
problem = fem.petsc.NonlinearProblem(F, ulm, bcs=[])
solver = nls.petsc.NewtonSolver(domain.comm, problem)
solver.solve(ulm)


(u, lm) = ulm.split()
with io.XDMFFile(domain.comm, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u)

with io.XDMFFile(domain.comm, "lm.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(lm)