Custom made non uniform mesh

mesh
1

Hello everyone,

I am using fenics-dolfinx 0.9.0 with Linux (Ubuntu 24.04.2 LTS).

I am trying to create a custom mesh to solve the mixed Poisson problem, following the approach described here: Mixed formulation for the Poisson equation.

I followed the steps outlined in this tutorial: Mesh generation, but when I try to plot_mesh the domain, the process takes too long and never finishes.

Here is my code:

from mpi4py import MPI

import matplotlib.pyplot as plt
import numpy as np
import pyvista

import basix.ufl
import dolfinx
import ufl

nodes = np.array([[1.0, 0.0], [2.0, 0.0], [3.0, 2.0], [1, 3]], dtype=np.float64)

connectivity = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)

c_el = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(nodes.shape[1],)))

domain = dolfinx.mesh.create_mesh(MPI.COMM_SELF, connectivity, nodes, c_el)

def plot_mesh(mesh: dolfinx.mesh.Mesh, values = None):
    """
    Given a DOLFINx mesh, create a `pyvista.UnstructuredGrid`,
    and plot it and the mesh nodes.

    Args:
        mesh: The mesh we want to visualize
        values: List of values indicating a marker for each cell in the mesh

    Note:
        If `values` are given as input, they are assumed to be a marker
        for each cell in the domain.
    """
    # We create a pyvista plotter instance
    plotter = pyvista.Plotter()

    # Since the meshes might be created with higher order elements,
    # we start by creating a linearized mesh for nicely inspecting the triangulation.
    V_linear = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    linear_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(V_linear))

    # If the mesh is higher order, we plot the nodes on the exterior boundaries,
    # as well as the mesh itself (with filled in cell markers)
    if mesh.geometry.cmap.degree > 1:
        ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
        if values is not None:
            ugrid.cell_data["Marker"] = values
        plotter.add_mesh(ugrid, style="points", color="b", point_size=10)
        ugrid = ugrid.tessellate()
        plotter.add_mesh(ugrid, show_edges=False)
        plotter.add_mesh(linear_grid,style="wireframe", color="black")
    else:
        # If the mesh is linear we add in the cell markers
        if values is not None:
            linear_grid.cell_data["Marker"] = values
        plotter.add_mesh(linear_grid,show_edges=True)

    # We plot the coordinate axis and align it with the xy-plane
    plotter.show_axes()
    plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        plotter.show()

pyvista.start_xvfb(1.0)
plot_mesh(domain)

Am I missing something? Do I need to make any modifications?

Thanks in advance!

2

You say you are running DOLFINx on linux.
How did you install DOLFINx, is it with conda, spack, apt or through docker?
Secondly, what happens if you remove the line:

1 Like
3

@dokken thank you very much. Removing the line pyvista.start_xvfb(1.0) it worked.

Let me ask one more question. How can i refine my mesh ? i found the function

refined_domain_uniform = dolfinx.mesh.refine(domain)

using that

nodes = np.array([[1.0, 0.0], [2.0, 0.0], [3.0, 2.0], [1, 3]], dtype=np.float64)
connectivity = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(nodes.shape[1],)))
domain = dolfinx.mesh.create_mesh(MPI.COMM_SELF, connectivity, nodes, c_el)

plot_mesh(domain)

refined_domain_uniform = dolfinx.mesh.refine(domain)
plot_mesh(refined_domain_uniform)

i get the following

Traceback (most recent call last):
  File "/home/vboxuser/Desktop/mixpoisson2.py", line 68, in <module>
    refined_domain_uniform = dolfinx.mesh.refine(domain)
                             ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/mesh.py", line 571, in refine
    mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
                                       ^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Edges must be initialised

Also i install DOLFINx with apt.

4

As the error message states, you need to initialize the edges of the mesh to be able to compute a refinement.
Edges are those topological entities of dimension 1.
See table 1 in: https://dl.acm.org/doi/pdf/10.1145/3524456
You can initialize these by calling

domain.topology.create_entities(1)

prior to calling refine.

1 Like
5

@dokken thank you for your response but i still getting an error

this is the code

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
import dolfinx
import basix
import pyvista
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
                 div, exp, inner)

def plot_mesh(mesh: dolfinx.mesh.Mesh, values = None):
    """
    Given a DOLFINx mesh, create a `pyvista.UnstructuredGrid`,
    and plot it and the mesh nodes.

    Args:
        mesh: The mesh we want to visualize
        values: List of values indicating a marker for each cell in the mesh

    Note:
        If `values` are given as input, they are assumed to be a marker
        for each cell in the domain.
    """
    # We create a pyvista plotter instance
    plotter = pyvista.Plotter()

    # Since the meshes might be created with higher order elements,
    # we start by creating a linearized mesh for nicely inspecting the triangulation.
    V_linear = dolfinx.fem.functionspace(mesh, ("CG", 1))
    linear_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(V_linear))

    # If the mesh is higher order, we plot the nodes on the exterior boundaries,
    # as well as the mesh itself (with filled in cell markers)
    if mesh.geometry.cmap.degree > 1:
        ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
        if values is not None:
            ugrid.cell_data["Marker"] = values
        plotter.add_mesh(ugrid, style="points", color="b", point_size=10)
        ugrid = ugrid.tessellate()
        plotter.add_mesh(ugrid, show_edges=False)
        plotter.add_mesh(linear_grid,style="wireframe", color="black")
    else:
        # If the mesh is linear we add in the cell markers
        if values is not None:
            linear_grid.cell_data["Marker"] = values
        plotter.add_mesh(linear_grid,show_edges=True)

    # We plot the coordinate axis and align it with the xy-plane
    plotter.show_axes()
    plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        plotter.show()

nodes = np.array([[1.0, 0.0], [2.0, 0.0], [3.0, 2.0], [1, 3]], dtype=np.float64)
connectivity = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)
c_el = ufl.Mesh(basix.ufl.element("CG", "triangle", 1, shape=(nodes.shape[1],)))
domain = dolfinx.mesh.create_mesh(MPI.COMM_SELF, connectivity, nodes, c_el)

plot_mesh(domain)

domain.topology.create_entities(1)

refined_domain_uniform = dolfinx.mesh.refine(domain)
plot_mesh(refined_domain_uniform)

and this is the error

Traceback (most recent call last):
  File "/home/vboxuser/Desktop/mixpoisson2.py", line 71, in <module>
    plot_mesh(refined_domain_uniform)
  File "/home/vboxuser/Desktop/mixpoisson2.py", line 36, in plot_mesh
    V_linear = dolfinx.fem.functionspace(mesh, ("CG", 1))
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py", line 607, in functionspace
    dtype = mesh.geometry.x.dtype
            ^^^^^^^^^^^^^
AttributeError: 'tuple' object has no attribute 'geometry'

Thank you again for your help

6

Well, now you need to check the output type of refine.

Refine returns the mesh, plus two optional arrays (used to map meshtags from one mesh to the other).

1 Like
7

Thank you again for your help.