Can I equate DoFs from different function spaces?

variational formulation
1

I have two or more function spaces with degrees of freedoms that coincide spatially and I want to equate the DoFs before solving the linear system.

In my case, I have three 1D meshes in a 2D plane that meet at a vertex. I want to replace the three DoFs at the vertex with a single one, effectively equating the value of that DoF beteen the function spaces. Is there a mechanism for this in Fenics?

I would expect it to be similar to the application of boundary conditions:

A = assemble(a)
b = assemble(L)
bc.apply(A, b)

but I have only found the possibility to set “u = g” for u in each function space individually.

2

This can be treated as a multipoint-constraint problem, and has an implementation in dolfinx, in the extension dolfinx_mpc.
However, is there a particular reason not for merging the meshes to be single 1D structure in a 2D space?

3

Thanks. I never tried dolfinx. I am interested in the recent mixed-dimensional additions to Fenics. I am not sure they are available in dolfinx?

I actually tried creating a single mesh and a CG function space on the forked 1D-mesh consisting of the three 1D-submeshes. My variational form consists of integrals over these 1D-subdomains of functions living in this 1D space in product with the trace of functions on 2D bulk domains for which these 1D subdomains are parts of their boundaries.

In principle, it should work fine. I am able to express that as a form, but I get the error below when trying to assemble the form, so I am a bit stuck on that route.

*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /home/fredrik/local/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  946dbd3e268dc20c64778eb5b734941ca5c343e5
4

Dolfinx does not have the mixed-dimensional assembly.
To be of further assistance, please provide a minimal code example showing how far you have gotten.
For instance consider a unitsquare with two 1D internal boundaries meeting at a point. I guess what you should do is to create the 2D mesh, then a single meshview for the union of all the internal boundaries that meet, and then use mesh-markers to subdivide the MeshView in assembly.

5

In my attempt to create a minimal example to reproduce the error I think I found the cause of the error that I pasted from petsc.

When there are two forms defined for the same pair of function spaces (block) but on different meshes, the sparsity pattern for the “GenericTensor” will be constructed based on only the mesh from the first form. When the other form (integrated over a different mesh) is about to be assembled, dolfin tries to add nonzeros to indices in the petsc matrix not present in the sparsity pattern.

Below is a minimal example that causes that error:

from dolfin import *

# Generate the meshes
square = UnitSquareMesh(8, 8)

# Mark the right boundary
facet_marker = MeshFunction("size_t", square, square.topology().dim()-1, 0)
for f in facets(square):
    facet_marker[f] = 1 - DOLFIN_EPS < f.midpoint().x() < 1 + DOLFIN_EPS
right_boundary = MeshView.create(facet_marker, 1)

# Mark the facets of the lower and upper part of that boundary
facet_marker = MeshFunction("size_t", square, square.topology().dim()-1, 0)
for f in facets(square):
    facet_marker[f] = ((1 - DOLFIN_EPS < f.midpoint().x() < 1 + DOLFIN_EPS)*
                       (1 + (0.5 - DOLFIN_EPS < f.midpoint().y())))
upper_right_boundary = MeshView.create(facet_marker, 1)
lower_right_boundary = MeshView.create(facet_marker, 2)

# Initialize function spaces for full domain on the full right boundary
V, V_boundary = [FunctionSpace(mesh, "CG", 1) for mesh in [square, right_boundary]]
W = MixedFunctionSpace(V, V_boundary)
u = TrialFunctions(W)
v = TestFunctions(W)

# Prepare to integrate over only parts of the right boundary
dsU = Measure("dx", domain=upper_right_boundary)
dsL = Measure("dx", domain=lower_right_boundary)

# Form
a = u[0]*v[1]*dsU + u[0]*v[1]*dsL
L = v[0]*dx

sol = Function(W)
assemble_mixed_system(a == L, sol)

If changing the form to, for instance

a = u[0]*v[0]*dsU + u[0]*v[1]*dsL

the assembly works fine. since there is at most 1 form per block.

The block (0, 1) in the matrix to the form in the first example seems to be assembled in two steps. On the line https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/dolfin/fem/MixedLinearVariationalSolver.cpp#lines-168:

assemble_mixed(*(As[i*u.size() + j]), *(a[i*u.size() + j][k]), bool(k>0));

Here i = 0, j = 1 and k = 0 for the first form and then k = 1 for the second. This is when it crashes.