Area of a closed curve mesh

dolfinx
1

Hi everyone,

I have constructed my mesh in this way:

# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 1.0e-3  # time step

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.

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

    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)
    model.mesh.setOrder(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.
    """
    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_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

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

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
P2 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)

ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))
V = functionspace(mesh, P2)

This is a manifold of geometrical dimension 2 but topological dimension 1. Is there a way to compute the area enclosed by this curve (in this case, it is just a circle, but the mesh evolves over time)?

Than you

2

Hey sdcardozof,

Tried to run your MWE but got an error at the following line due to the presence of gdim=2:

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)

Regardless, I believe the following line of code should compute the area of your mesh (although double check this):

area = fem.assemble_scalar(fem.form(fem.Constant(mesh, 1.0)*dx))

You may need to define your measure as well if you want to do more specific calculations with it.

I will also provide the code I used to generate a unit circle down below:

# Initialise the GMSH module
gmsh.initialize()

# Create the membrane and start the computations by the GMSH CAD kernel, to generate the relevant underlying data structures
# Create Circular Mesh with area of A0 = 1
membrane = gmsh.model.occ.addDisk(0, 0, 0, np.sqrt(1/np.pi), np.sqrt(1/np.pi)) # x, y, z, x-radius, y-radius
gmsh.model.occ.synchronize()

# Make membrane a physical surface for GMSH to recognise when generating the mesh
gdim = 2 # 2D Geometric Dimension of the Mesh
gmsh.model.addPhysicalGroup(gdim, [membrane], 0) # Dimension, Entity tag, Physical tag

# Generate 2D Mesh with uniform mesh size
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.023)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.023)
gmsh.model.mesh.generate(gdim)

# Create a mesh from the GMSH model
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)
domain.name = "initial_mesh"
gmsh.write("circle.msh")

# Finalize GMSH
gmsh.finalize()

Cheers,
Volkan