Solving mixed poisson in circunference

Hello everyone,

I am implementing a mixed Poisson problem on the circumference (exterior of the circle, which means 1D elements). Here is my code:

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    Args:
        model: Gmsh model to add the mesh to.
        name: Name (identifier) of the mesh to add.

    Returns:
        Gmsh model with a circle mesh added.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)

    circle= model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=2)
    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.

    Args:
        comm: MPI communicator top create the mesh on.
        model: Gmsh model.
        name: Name (identifier) of the mesh to add.
        filename: XDMF filename.
        mode: XDMF file mode. "w" (write) or "a" (append).

    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "Circle")
model.setCurrent("Circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCirc_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCirc_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

gmsh.finalize

# Define approximation spaces and basis functions
k = 1
Q_el = element("Lagrange", mesh.basix_cell(), k)
P_el = element("Lagrange", mesh.basix_cell(), k)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(mesh, V_el)


However, I encounter this error:

ValueError: Non-matching UFL cell and mesh cell shapes.

for the last line.

I have already changed some parameters of the mesh generation, but I still haven’t found a solution.

Does anyone know how to fix this?

Thanks.

Should the circumferencr be embedded in 2D or 3D space? You are omittimg the gdim kwarg, which by default is 3.

Depending on your version of dolfinx, you Also need to specify gdim in the Basix element and mixed element construction (this has been removed on the main branch).

The circumference is embedded in 2D. Did you mean to change this line?

msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0, gdim=1)

How should I define gdim in the Basix element?

Also, I tried to define explicitly:

k = 1

Q_el = element("Lagrange", "interval", k)
P_el = element("Lagrange", "interval", k)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(mesh, V_el)

And I have the same error

Gdim should be 2 if you are embedding in 2D

You can send in gdim as a keyword argument to Basix.ufl.element and Basix.ufl.mixed_element.

Note that this keyword argument will be removed on 0.8.x

Thanks for your response Jorgen. Since this keyword will be removed, how should I define this type element?

I follow your suggestions and it doesn’t get this error anymore, but, now with the solution I encountered this error

AttributeError: ‘LinearProblem’ object has no attribute ‘_solver’

I am practically following this tutorial with some changes.

Anyway, here is my code

# Define approximation spaces and basis functions
k = 1

Q_el = element("Lagrange", mesh.basix_cell(), k, gdim=2)
P_el = element("DG", mesh.basix_cell(), k-1, gdim=2)
R_el = element("DG", mesh.basix_cell(),k-1, gdim=2)
V_el = mixed_element([Q_el, P_el, R_el],gdim=2)
V = fem.functionspace(mesh, V_el)

(sigma, u, r) = TrialFunctions(V)
(tau, v, t) = TestFunctions(V)

x = ufl.SpatialCoordinate(mesh)

f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx + inner(r,v) * dx + inner(t,u)* dx
L = inner(f, v) * dx


problem = LinearProblem(a, L)
uh = problem.solve()

I already reviewed the documentation and I can’t figure out the reason.

On the main branch the geometrical dimension is interpreted from the mesh, so you simply dont need to pass anything if mesh.geometry.dim is correct

For your new issue, i need a reproducible code; ie, include mesh generation, imports and loading of the mesh in a single script

Ok, thank you,

Certainly! Below is my full code

from mpi4py import MPI
from petsc4py.PETSc import ScalarType  # type: ignore

import numpy as np
from pathlib import Path

import ufl

from basix.ufl import element, mixed_element
from dolfinx.io import XDMFFile, gmshio
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner, dx)


try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    Args:
        model: Gmsh model to add the mesh to.
        name: Name (identifier) of the mesh to add.

    Returns:
        Gmsh model with a circle mesh added.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)

    circle= model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.

    Args:
        comm: MPI communicator top create the mesh on.
        model: Gmsh model.
        name: Name (identifier) of the mesh to add.
        filename: XDMF filename.
        mode: XDMF file mode. "w" (write) or "a" (append).

    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0, gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "Circle")
model.setCurrent("Circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCirc_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCirc_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

gmsh.finalize

# Define approximation spaces and basis functions
k = 1

Q_el = element("Lagrange", mesh.basix_cell(), k, gdim=2)
P_el = element("DG", mesh.basix_cell(), k-1, gdim=2)
R_el = element("DG", mesh.basix_cell(),k-1, gdim=2)
V_el = mixed_element([Q_el, P_el, R_el],gdim=2)
V = fem.functionspace(mesh, V_el)

(sigma, u, r) = TrialFunctions(V)
(tau, v, t) = TestFunctions(V)

x = ufl.SpatialCoordinate(mesh)

f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)

a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx + inner(r,v) * dx + inner(t,u)* dx
L = inner(f, v) * dx

problem = LinearProblem(a, L)
uh = problem.solve()

with io.XDMFFile(mesh.comm, "out_poisson/poisson.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)

The issue is + inner(div(sigma), v) as sigma is a scalar value`, the divergence is undefined.

1 Like