Unable to solve my Lagrange multiplier problem with either dolfinx or dolfin

Legacy DOLFIN
, ,
1

I have a variational formulation that uses the Lagrange multiplier method. I need to use a real element for this multiplier. I’m a bit confused by real elements being in FEniCSx but still not implemented: calling

lagrange_elem = basix.ufl.element("Lagrange", domain.ufl_cell().cellname(), 1)
real_elem = basix.ufl.real_element(domain.ufl_cell().cellname(), (1,))
mixed_elem = basix.ufl.mixed_element([lagrange_elem, real_elem])

is fine, but

V = fem.functionspace(domain, mixed_elem)

is not for some reason and raises NotImplementedError. I’m aware of this issue but it’s contradicted by this one so I’m a bit confused.

So I decided to try my luck with legacy Dolfin. My problem now is I can’t find how to generate a circular mesh with holes in legacy Dolfin. It was easy in 1.3 with CSG, or in Dolfinx with Gmsh, but I can’t find anything in the 1.9 documentation about that.

How can I create the mesh without relying on an external xdmf file? I need to loop on different meshes that I generate on the fly and the performance will be hit too hard if I need to handle file reading and writing.

My problem is as follows, with \eta a scalar function, \lambda a real number and i\in\{1,2\}:

\int_\Omega \nabla \eta \cdot \nabla \varphi d\Omega + \int_\Omega \eta\varphi d\Omega + \int_{\Gamma_0}\eta\varphi d\Gamma-\sum_i \tan\theta_i \int_{\Gamma_i} \varphi d\Gamma\\+2\lambda\sum_i\int_{\Gamma_i}\nabla\eta\cdot\nabla\varphi d\Gamma + \mu\sum_i\int_{\Gamma_i}\left[\nabla\eta^2-\tan\theta_i^2\right]d\Gamma_i=0

image

The differential equation is (\Delta-1)\eta=0 and the conditions are

\nabla \eta \cdot \mathbf{t} = 0 \text{ on } \Gamma_i\\ \nabla \eta\cdot\mathbf{n}=\tan\theta_i \text{ on } \Gamma_i\\ \nabla \eta \cdot \mathbf{n} = -\eta \text{ on } \Gamma_0

Maybe there is a different way to solve this?

Here is my Dolfinx code:

import basix.ufl
from mpi4py import MPI

import numpy as np

import ufl
import basix
import dolfinx
from dolfinx.fem.petsc import LinearProblem

tan_theta_1 = 1
tan_theta_2 = 1
interchannel_distance = 4
channel_radius = 1
domain_radius = 20

### MESH ###

import gmsh
gmsh.initialize()

# Geometry
ext_boundary = gmsh.model.occ.addCircle(0, 0, 0, domain_radius)
ext_boundary_loop = gmsh.model.occ.addCurveLoop([ext_boundary])

channel1 = gmsh.model.occ.addCircle(-interchannel_distance/2, 0, 0, channel_radius)
channel1_loop = gmsh.model.occ.addCurveLoop([channel1])

channel2 = gmsh.model.occ.addCircle(interchannel_distance/2, 0, 0, channel_radius)
channel2_loop = gmsh.model.occ.addCurveLoop([channel2])

membrane = gmsh.model.occ.addPlaneSurface([ext_boundary_loop, channel1_loop, channel2_loop])
gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)

# Mesh generation
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 100)

gmsh.model.mesh.generate(gdim)

# DOLFINx interface
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

### Boundaries ###
from dolfinx.mesh import locate_entities, meshtags

# Markers
def is_on_ext_boundary(x) -> bool:
    return np.isclose(x[0]**2 + x[1]**2, domain_radius**2)

def is_on_channel1_boundary(x) -> bool:
    return np.isclose((x[0] + interchannel_distance/2)**2 + x[1]**2, channel_radius**2)

def is_on_channel2_boundary(x) -> bool:
    return np.isclose((x[0] - interchannel_distance/2)**2 + x[1]**2, channel_radius**2)

def is_on_boundary(x) -> bool:
    return np.logical_or(is_on_ext_boundary(x), is_on_channel1_boundary(x), is_on_channel2_boundary(x))

boundaries = list(enumerate([is_on_ext_boundary, is_on_channel1_boundary, is_on_channel2_boundary]))

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Debugging
from dolfinx.io import XDMFFile
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags.xdmf", "w", encoding=XDMFFile.Encoding.ASCII) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)

### Variational formulation ###

# Function space
from dolfinx import fem
V = fem.functionspace(domain, ("Lagrange", 1))

# Build function space with Lagrange multiplier
lagrange_elem = basix.ufl.element("Lagrange", domain.ufl_cell().cellname(), 1)
real_elem = basix.ufl.real_element(domain.ufl_cell().cellname(), (1,))



lagrange_V = fem.functionspace(domain, lagrange_elem)
real_V = fem.functionspace(domain, real_elem) # NotImplementedError
mixed_V = fem.functionspace(domain, lagrange_V * real_V)
2

IMHO, you are mixing up too many questions.

You start saying that you would like a “real” element in dolfinx. But, it wasn’t available, so you moved back to dolfin. Still, midway through the post the issue seems to be that you can’t mesh your domain. Cherry on top, when showing us the code, you actually post a dolfinx code, while you said you had moved to dolfin. Overall, the reader may be confused.

  • on the real element being unavaible: that is still the case. Search the forum with the word nullspace to find alternative solutions. Your problem may also benefit from using the dolfinx-mpc addon library.
  • on creating a hole in gmsh: see gmsh documentation, or Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken, where the author uses gmsh.model.occ.cut.
1 Like
3

Thank you for your answer, I will look into “null space”.
I didn’t post any dolfin code because I don’t even have a mesh (which I do in Dolfinx thanks to the Gmsh interface). The tutorial you linked to uses Dolfinx to interact with Gmsh, can I use a Dolfinx mesh in legacy Dolfin? If so, then my problem is solved.

4

See Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio to use a gmsh mesh in dolfin (without requiring dolfinx at all). You may want to start reading it from the bottom, since there have been many replies.

1 Like
5

It works with MPC!

image
image751×215 61.7 KB

from typing import Callable
from mpi4py import MPI

import numpy as np

import gmsh

import ufl
# import dolfinx
from dolfinx import fem, plot
from dolfinx.io import gmshio, XDMFFile
from dolfinx.mesh import locate_entities, locate_entities_boundary, meshtags
from dolfinx_mpc import LinearProblem, MultiPointConstraint

import pyvista

tan_theta_1 = 1
tan_theta_2 = 1
interchannel_distance = 4
channel_radius = 1
domain_radius = 20

### MESH ###

gmsh.initialize()

# Geometry
ext_boundary = gmsh.model.occ.addCircle(0, 0, 0, domain_radius)
ext_boundary_loop = gmsh.model.occ.addCurveLoop([ext_boundary])

channel1 = gmsh.model.occ.addCircle(-interchannel_distance/2, 0, 0, channel_radius)
channel1_loop = gmsh.model.occ.addCurveLoop([channel1])

channel2 = gmsh.model.occ.addCircle(interchannel_distance/2, 0, 0, channel_radius)
channel2_loop = gmsh.model.occ.addCurveLoop([channel2])

membrane = gmsh.model.occ.addPlaneSurface([ext_boundary_loop, channel1_loop, channel2_loop])
gmsh.model.occ.synchronize()

gdim = 2
gmsh.model.addPhysicalGroup(gdim, [membrane], 1)

# Mesh generation
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 100)

gmsh.model.mesh.generate(gdim)

# DOLFINx interface

gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

### Boundaries ###

# Markers
def is_on_ext_boundary(x) -> bool:
    return np.isclose(x[0]**2 + x[1]**2, domain_radius**2)

def is_on_channel1_boundary(x) -> bool:
    return np.isclose((x[0] + interchannel_distance/2)**2 + x[1]**2, channel_radius**2)

def is_on_channel2_boundary(x) -> bool:
    return np.isclose((x[0] - interchannel_distance/2)**2 + x[1]**2, channel_radius**2)

def is_on_boundary(x) -> bool:
    return np.logical_or(is_on_ext_boundary(x), is_on_channel1_boundary(x), is_on_channel2_boundary(x))

boundaries = list(enumerate([is_on_ext_boundary, is_on_channel1_boundary, is_on_channel2_boundary]))

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# Debugging
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, "facet_tags.xdmf", "w", encoding=XDMFFile.Encoding.ASCII) as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(facet_tag, domain.geometry)

### Variational formulation ###

# Function space
V = fem.functionspace(domain, ("Lagrange", 1))

facets = locate_entities_boundary(
    domain,
    dim=(domain.topology.dim - 1),
    marker=is_on_boundary,
)

dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)

# Multi-point constraint
master_x_channel1 = (-interchannel_distance/2, channel_radius)
master_x_channel2 = (interchannel_distance/2, channel_radius)

def slave_x_channel1(x_slave:np.ndarray) -> bool:
    return np.logical_and(
        is_on_channel1_boundary(x_slave), 
        np.logical_or(
            np.logical_not(np.isclose(x_slave[0], master_x_channel1[0])), 
            np.logical_not(np.isclose(x_slave[1], master_x_channel1[1]))
            )
        )

def slave_x_channel2(x_slave:np.ndarray) -> bool:
    return np.logical_and(
        is_on_channel2_boundary(x_slave), 
        np.logical_or(
            np.logical_not(np.isclose(x_slave[0], master_x_channel2[0])), 
            np.logical_not(np.isclose(x_slave[1], master_x_channel2[1]))
            )
        )

def circular_periodic_condition(x_master:tuple[float,float]) -> Callable[[np.ndarray], np.ndarray]:
    def condition(x):
        vals = np.zeros(x.shape)
        vals[0] = x_master[0]
        vals[1] = x_master[1]
        return vals
    return condition

mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(
    V,
    slave_x_channel1,
    circular_periodic_condition(master_x_channel1),
    bcs=[]
    )
mpc.create_periodic_constraint_geometrical(
    V,
    slave_x_channel2,
    circular_periodic_condition(master_x_channel2),
    bcs=[]
    )
mpc.finalize()

# Measure on the boundary
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)

a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(u, v) * ufl.dx + ufl.dot(u, v) * ds(0)
L = tan_theta_1 * v * ds(1) + tan_theta_2 * v * ds(2)

problem = LinearProblem(a, L, mpc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

with XDMFFile(domain.comm, "out_membrane.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh)



cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
if pyvista.OFF_SCREEN:
    pyvista.start_xvfb(wait=0.1)
    plotter.screenshot("uh_membrane.png")
else:
    plotter.show()