Long-range coupling

dolfinx
1

Hi,

I need to implement a long-range coupling between disconnected subdomains, that represent springs between them. My problem is similar to this simple one:

\Omega = \Omega_1 \cup \Omega_2,
\int_\Omega \nabla u : \nabla \delta u = \int_\Omega \delta u q,
with the spring force q(x) = -k(u(x) - u(T(x))),
k a parameter,
and T(x) a mapping from \Omega_1 to \Omega_2 and from \Omega_2 to \Omega_1.

The way I am solving it now and I’d like to improve involves working with an extended function space that contains copies of all subdomain fields in each subdomain. Then, I use dolfinx_mpc to express that all these are copies, and I solve the problem involving only the master dofs.

So, for the above equations, that would be:

elt = basix.ufl.element(...)
N = 2
# V = fem.functionspace(domain, elt)  # that would be the "natural" function space
V_big = fem.functionspace(domain, basix.ufl.mixed_element([elt] * N))  # here each subdomain is copied on all other subdomains
u = ufl.TrialFunctions(V_big)  # len(...) = N -> u_n=u[n]
du = ufl.TestFunctions(V_big)

# K = ufl.inner(ufl.grad(u), ufl.grad(du)) * dx  # if dealing with V
K = ufl.inner(ufl.grad(u[0]), ufl.grad(du[0])) * dx(0) + ufl.inner(ufl.grad(u[1]), ufl.grad(du[1])) * dx(1)

# long-range coupling
q0 = -k * (u[0] - u[1])
q1 = -k * (u[1] - u[0])
W_ext = ufl.inner(q0, du[0]) * dx(0) + ufl.inner(q1, du[1]) * dx(1)

K_total = K - W_ext

and then comes a dolfinx_mpc.MultiPointConstraint(V_big) to express that u[0] on tag 1 is mapped to u[0] on tag 0 through a given mapping T, and similar for u[1].

Now, this works, but I don’t find it very elegant so I am wondering whether there could be a better way. Drawbacks: the slave dofs grow as N^2 for N subdomains; the “classical” part of the equation (here K) does not express very nicely; and mostly I am forced to treat this as a multiphysics problem if there are other subdomains that are not long-range coupled.
I think ideally I would like to work with two function spaces simultaneously, and define and mapping between them. Say, something like:

V = fem.functionspace(domain, elt)
u = ufl.TrialFunction(V)
du = ufl.TestFunction(V)
K = ufl.inner(ufl.grad(u), ufl.grad(du)) * dx

# long-range coupling
fictitious_subdomain = ...
V_not_bigger = fem.functionspace(fictitious_subdomain, basix.ufl.mixed_element([elt] * N))
u_f = ufl.TrialFunctions(V)
du_f = ufl.TestFunctions(V)

q0 = -k * (u_f[0] - u_f[1])
q1 = -k * (u_f[1] - u_f[0])
W_ext = ufl.inner(q0, du_f[0]) * dx + ufl.inner(q1, du_f[1]) * dx

# and then a way to define the mapping T(x) between u{tag n} and u_f[n]

Is this doable?

Nonlinear Poisson with periodic coupling
2

Is the mapping T(x) explicit, that in a sense, given an x in \Omega_i it gives you y=T(x) in \Omega_j, i\neq j?

I think this could be done with usage of some global interpolation operators, similar to what one does for non-conforming 3D-1D coupling in fenics_ii (GitHub - MiroK/fenics_ii) and fenicsx_ii (GitHub - scientificcomputing/fenicsx_ii: A reimplementation of some of the functionality from MiroK's FEniCS_ii (FEniCS trace))

1 Like
3

Yes it is explicit. Thank you for the links, I’ll have a look.

4

I’ve made a small demo illustrating this kind of coupling for coefficients (on nonmatching grids with non-matching cell type) at:

It should also work for trial and or test functions, which would make it possible to set up your problem in a UFL syntax (see: A coupled 3D-1D manufactured solution — fenicsx_ii for further details on bilinear and linear forms).

2 Likes
5

This looks great! I’ll try.

6

Hi again,

I just tried to replace the fem.Function in your example with ufl.Trial / TestFunction and I can’t make it work.

Is this usage what you suggested?

# instead of
# u_0 = dolfinx.fem.Function(V0, name="u0")
# I define
u0, v0 = ufl.TrialFunction(V0), ufl.TestFunction(V0)

# and then
Q1 = dolfinx.fem.functionspace(mesh1, quadrature)
u1, v1 = ufl.TrialFunction(Q1), ufl.TestFunction(Q1)

u0_on_omega1 = Average(u0, restriction, Q1)
v0_on_omega1 = Average(v0, restriction, Q1)

dx1 = ufl.Measure("dx", domain=mesh1)
x1 = ufl.SpatialCoordinate(mesh1)

# This...
a = ufl.inner(u0_on_omega1, v1) * dx1
# ...  does not work
aform = dolfinx.fem.form(a)
>>> ..... fenicsx_ii/ufl_operations.py:141, in Average.__init__ .....
>>> TypeError: Can only average arguments or coefficients, got <class 'ufl.referencevalue.ReferenceValue'>

# Neither...
a = ufl.inner(u0_on_omega1, v0_on_omega1) * dx1
# ... does this
aform = dolfinx.fem.form(a)
>>> same error here
7

You can’t use dolfinx.fem.form directly, as there are some intermediate steps encoded into fenicsx_ii to handle this.

If you try calling fenicsx_ii.assemble_matrix, it will work:

from typing import Callable

from mpi4py import MPI

import basix.ufl
import dolfinx
import numpy as np
import numpy.typing as npt
import ufl

from fenicsx_ii import Average,  assemble_matrix
from fenicsx_ii.quadrature import Quadrature

# -

# Next we create the two domains

# +
mesh0 = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, cell_type=dolfinx.mesh.CellType.triangle
)

translation_vector = np.array([4, 3])

mesh1 = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD,
    [translation_vector, translation_vector + np.ones(2)],
    [20, 20],
    cell_type=dolfinx.mesh.CellType.quadrilateral,
)
# -

# Furthermore, we implement the translation-restriction from $\Omega_1$ to $\Omega_0$
# as a general mapping restriction, using the {py:class}`fenicsx_ii.ReductionOperator`
# protocol, where we have to implement a mapping that takes a set of reference points
# in $\mathbb{R}^2$ and a set of cells in $\Omega_1$ and computes the mapping
# $T(F(x_{ref}))$ where $F$ is the push forward operation from $K_{ref}$
# (reference element) to $K$ (element in physical space) and $T$ the mapping operator.


class MappedRestriction:
    quadrature_degree: int

    def __init__(
        self,
        mesh,
        operator: Callable[[npt.NDArray[np.floating]], npt.NDArray[np.floating]],
    ):
        self._mesh = mesh
        self._operator = operator

    def compute_quadrature(
        self, cells: npt.NDArray[np.int32], reference_points: npt.NDArray[np.floating]
    ) -> Quadrature:
        phys_points = np.zeros(
            (len(cells) * reference_points.shape[0], self._mesh.geometry.dim)
        )
        x_geom = self._mesh.geometry.x[
            self._mesh.geometry.dofmap[cells], : self._mesh.geometry.dim
        ]
        for i, x_i in enumerate(x_geom):
            phys_points[
                i * reference_points.shape[0] : (i + 1) * reference_points.shape[0], :
            ] = self._mesh.geometry.cmap.push_forward(reference_points, x_i)
        translated_points = self._operator(phys_points.T).T
        return Quadrature(
            name="Translation",
            points=translated_points,
            weights=np.ones((phys_points.shape[0], 1)),
            scales=np.ones(phys_points.shape[0]),
        )

    @property
    def num_points(self) -> int:
        return 1


# Next we specify the actual translation operator from $\Omega_1$ to $\Omega_0$ and
# create the restriction operator.


# +
def translate_1_to_0(x):
    x_out = x.copy()
    for i, ti in enumerate(translation_vector):
        x_out[i] -= ti
    return x_out


restriction = MappedRestriction(mesh1, translate_1_to_0)
# -

# Next, we are ready to test the implementation.
# We start by defining $u_0$ as a known function.


# +
def f(x):
    return x[0] + 2 * x[1] * x[1]


V0 = dolfinx.fem.functionspace(mesh0, ("Lagrange", 2))
u0, v0 = ufl.TrialFunction(V0), ufl.TestFunction(V0)

# and then
q_degree = 2
quadrature = basix.ufl.quadrature_element(
    mesh1.basix_cell(), value_shape=(), degree=q_degree
)
Q1 = dolfinx.fem.functionspace(mesh1, quadrature)
u1, v1 = ufl.TrialFunction(Q1), ufl.TestFunction(Q1)

u0_on_omega1 = Average(u0, restriction, Q1)
v0_on_omega1 = Average(v0, restriction, Q1)

dx1 = ufl.Measure("dx", domain=mesh1)
x1 = ufl.SpatialCoordinate(mesh1)

# This...
a = ufl.inner(u0_on_omega1, v1) * dx1
# ...  does not work
A = assemble_matrix(a)
A.assemble()
# # Neither...
a = ufl.inner(u0_on_omega1, v0_on_omega1) * dx1
A11 = assemble_matrix(a)
A11.assemble()

If you want to create the matrix once, you can use:

8

Oh, I should have seen it, sorry. It works, thanks a lot.