Triangle elements in 3D broken in dolfinx?

General
1

Hi all,

After a recent update to latest nightly dolfinx build, I realize that an example where I have “triangle3D” elements isn’t working anymore. MWE to reproduce:

#!/usr/bin/env python3
from mpi4py import MPI
from dolfinx import fem, io
import ufl

comm = MPI.COMM_WORLD

meshname='triangle3D.xdmf'

with io.XDMFFile(comm, meshname, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    mesh = infile.read_mesh(name="Grid")

V_s = fem.FunctionSpace(mesh, ("CG", 1))
V_v = fem.VectorFunctionSpace(mesh, ("CG", 1))
V_t = fem.TensorFunctionSpace(mesh, ("CG", 1))

where the mesh file “triangle3D.xdmf” is

<Xdmf Version="3.0">
  <Domain>
    <Grid Name="Grid">
      <Geometry GeometryType="XYZ">
        <DataItem DataType="Float" Dimensions="3 2" Format="XML" Precision="8">
          0.0 0.0 0.0
          1.0 0.0 0.0
          0.0 1.0 0.0
        </DataItem>
      </Geometry>
      <Topology NodesPerElement="3" NumberOfElements="1" TopologyType="triangle">
        <DataItem DataType="Int" Dimensions="1 3" Format="XML" Precision="4">
          0 1 2
        </DataItem>
      </Topology>
    </Grid>
  </Domain>
</Xdmf>

Problem: CG vector and tensor function space creations fail with the error ufl.log.UFLException: Non-matching cell of finite element and domain., whereas scalar funtion space creation works.

Remedy would be to use a 2D triangle specifying only “GeometryType=“XY”” in the xdmf file. But, I wonder if some functionality has been removed for the triangle3D case, or if this is a bug. Because I can confirm that this example worked for dolfinx from 28 Aug 2022 or earlier.

Thanks!

Best,
Marc

1 Like
2

Thanks for spotting this issue. I’ve reported it at: Manifold meshes from XDMF broken · Issue #2389 · FEniCS/dolfinx · GitHub

1 Like
3

Thanks for reporting this! I just updated dolfinx and can confirm the MWE works now. However, I had a case where I applied DBCs to the triangle’s edges’ Z direction, and this is yielding an assertion error. Consider the MWE with the same mesh as posted above:

#!/usr/bin/env python3
from mpi4py import MPI
from dolfinx import fem, io, mesh
import ufl

# function to mark z=0
def twodimZ(x):
    return np.isclose(x[2], 0.0)

comm = MPI.COMM_WORLD

meshname='triangle3D.xdmf'

with io.XDMFFile(comm, meshname, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    msh = infile.read_mesh(name="Grid")

V_s = fem.FunctionSpace(msh, ("CG", 1))
V_v = fem.VectorFunctionSpace(msh, ("CG", 1))
V_t = fem.TensorFunctionSpace(msh, ("CG", 1))

u = fem.Function(V_v)

dbc = fem.dirichletbc(u.sub(2), fem.locate_dofs_topological(V_v.sub(2), msh.topology.dim-1, mesh.locate_entities_boundary(msh, msh.topology.dim-1, twodimZ)))

Last line yields an AssertionError “assert self.ufl_element().num_sub_elements() > i”.

Is there anything wrong / has something changed in the syntax of applying the DBC in Z here?

Thanks,
Marc

1 Like
4

I reopened the issue

1 Like
5

Should hopefully be resolved with: Pass manifold cell information to base element of vector elements by jorgensd · Pull Request #617 · FEniCS/basix · GitHub
If you could test it I would be grateful.

6

I am facing the same problem:

from mpi4py import MPI
from dolfinx import fem, io
from dolfinx.mesh import create_mesh
import meshio
from ufl import FiniteElement, Cell, Mesh, VectorElement

points = [
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [1.0, 1.0, 0.0]
]
cells = [("triangle", [[0, 1, 2], [1, 3, 2]])]
mesh = meshio.Mesh(points, cells)

shape = "triangle"
degree = 1
cell = Cell(shape)
domain = Mesh(VectorElement("Lagrange", cell, degree))
mesh = create_mesh(MPI.COMM_WORLD, mesh.cells_dict['triangle'], mesh.points, domain)

vec_elem = VectorElement("CG", mesh.ufl_cell(), 2)
fin_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, vec_elem)
Q = fem.FunctionSpace(mesh, fin_elem)

@marchirschvogel I tried it with fenics-dolfinx=0.5.0 since this is the release before your mentioned date 28 Aug 2022. But the error:

ufl.log.UFLException: Non-matching cell of finite element and domain.
still occurs.
Is there any workaround until the bug is fixed? I tried multiple approaches but this bug seems to prevent all 3D simulations on custom triangle meshes.

Many thanks in advance.

7

The bug has been resolved since the 0.5.2 release (on the main branch).
The bug does not affect all simulations on custom triangle meshes, only those on manifolds (i.e. a triangle embedded in 3D space). If you want to work on manifolds you need to use the main branch (or the nightly docker image).

The example you supply can be fixed by calling:

mesh = meshio.Mesh(points[:,:2], cells)

instead of

8

@dokken Thank you for the quick reply. I am already using the latest release fenics-dolfinx=0.5.2, reproduced the bug there and on the latest dolfinx/dolfinx:nightly docker image. I just wanted to show that the bug also occurs in fenics-dolfinx=0.5.0.

Your suggested fix just reduces my code from 3D into 2D space. But my code is just a minimal example to show the bug. My actual use case is in 3D (not reducible to 2D).
Anyway I found a workaround for creating the mesh but this introduced another bug.
Minimal example to reproduce:

from mpi4py import MPI
from dolfinx import fem, io
from petsc4py import PETSc
import numpy as np
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import locate_dofs_topological, dirichletbc
import meshio
from ufl import FiniteElement, VectorElement


def boundary(x):
    return np.full(x.shape[1], True)


points = [
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [1.0, 1.0, 0.0]
]
cells = [("triangle", [[0, 1, 2], [1, 3, 2]])]
tri_mesh = meshio.Mesh(points, cells)
meshio.write("tri-bug-mesh.xdmf", tri_mesh)
with io.XDMFFile(MPI.COMM_WORLD, "tri-bug-mesh.xdmf", 'r') as infile:
    mesh = infile.read_mesh(name="Grid")

vec_elem = VectorElement("CG", mesh.ufl_cell(), 2)
fin_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
V = fem.FunctionSpace(mesh, vec_elem)
Q = fem.FunctionSpace(mesh, fin_elem)

# Boundary conditions
facetdim = mesh.topology.dim - 1
bndry_facets = locate_entities_boundary(mesh, facetdim, boundary)
wall_dofs = locate_dofs_topological(V, facetdim, bndry_facets)
u_noslip = np.array((0,) * mesh.topology.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

Output:
> Traceback (most recent call last):
>   File "/home/user/thrombectomics/is_it_a_bug.py", line 42, in <module>
>     bc_noslip = dirichletbc(u_noslip, wall_dofs, V)
>   File "/home/user/anaconda3/envs/thrombectomics/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 182, in dirichletbc
>     return formcls(value, dofs, V)
>   File "/home/user/anaconda3/envs/thrombectomics/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 127, in __init__
>     super().__init__(_value, dofs, V._cpp_object)  # type: ignore
> RuntimeError: Creating a DirichletBC using a Constant is not supported when the Constant size is not equal to the block size of the constrained (sub-)space. Use a fem::Function to create the fem::DirichletBC.

Maybe I misunderstand something about when to use triangles or tetrahedrons in 3D. I want to do CFD in a pipe created from a custom mesh which consists of points and triangle-faces. Seems like a common use case to me.

Many thanks in advance.

9

We are not backporting fixes, thus to use main branch of dolfin if you want to use a manifold.

So you want the facets of the mesh to be triangles? Then your mesh should consist of tetrahedrons.

This is still related to you trying to do something invalid.
As you send in your mesh as a 3D mesh (points are 3D, thus your mesh consists of triangular cells embedded in 3D). and you are sending in the topological dimension (which is equal to 2) to your DirichletBC

10

@dokken

So you want the facets of the mesh to be triangles? Then your mesh should consist of tetrahedrons.

I am confused about how to reorganize my 3D-triangle mesh to tetrahedrons. Since a tetrahedron has one additional point compared to a triangle, where I am supposed to put this 4th point so it works with fenicsx/dolfinx? I already tried pygalmesh.generate_volume_mesh_from_surface_mesh() but the resulting CFD looks like dolfinx is using all the tetrahedrons inside my pipe mesh as walls, although it should flow through the pipe.

Unfortunately I didn’t find any documentation about this.
Many thanks in advance.

11

Could you please provide a picture of the domain you would like to simulate something on? I really do not understand how you can simulate a flow in a pipe (assuming something cylindrical) with triangular elements. Or do you mean that you have a 2D slice of a pipe, similar to Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial?
Are the z-coordinates of your points anything else than a constant value?

12

@dokken Sure!

airways_coronal
airways_coronal863×939 155 KB

airways_inflow
airways_inflow1813×988 188 KB

Visualized with pyvista using:

verts, faces, normals, values = skimage.measure.marching_cubes(img, step_size=1)
verts /= max(img.shape)

# Meshio
cells = [("triangle", faces)]
mesh = meshio.Mesh(verts, cells)
vtu_path = 'triangle_surface.vtu'
meshio.write(vtu_path, mesh)
# faces = np.insert(faces, 0, np.ones(faces.shape[0]), axis=1)
mesh = pygalmesh.generate_volume_mesh_from_surface_mesh(
    "triangle_surface.vtu",
    max_facet_distance=0.001,
)

# Dolfinx
vec_elem = VectorElement("Lagrange", tetrahedron, 1)
domain = Mesh(vec_elem)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, mesh.cells_dict['tetra'], mesh.points, domain)

m_topology, m_cell_types, m_geometry = create_vtk_mesh(mesh)
m_grid = pyvista.UnstructuredGrid(m_topology, m_cell_types, m_geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(m_grid, show_edges=True, opacity=0.5)
plotter.show_bounds(
    grid='front',
    location='outer',
    all_edges=True)
plotter.view_xz()
plotter.show()

The CFD is supposed to be inside the “pipes”.

Let me know if I should create a new thread for this topic.

13

Then you need to mesh the interior or the pipe.

2 Likes
14

Sorry for the late reply on this. I can confirm that the MWE

#!/usr/bin/env python3
from mpi4py import MPI
from dolfinx import fem, io, mesh
import ufl

# function to mark z=0
def twodimZ(x):
    return np.isclose(x[2], 0.0)

comm = MPI.COMM_WORLD

meshname='triangle3D.xdmf'

with io.XDMFFile(comm, meshname, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    msh = infile.read_mesh(name="Grid")

still fails with latest nightly dolfinx Docker image, with error python3: /src/dolfinx/cpp/dolfinx/io/xdmf_mesh.cpp:264: std::pair<std::vector, std::array<long unsigned int, 2> > dolfinx::io::xdmf_mesh::read_geometry_data(MPI_Comm, hid_t, const pugi::xml_node&): Assertion `gdims[1] == (int)gdim’ failed.

Mesh file from opening post.

Best,
Marc

15

There is a typo in your initial mesh,

The dimensions here should be 3 by 3 right? not 3 by 2.

1 Like
16

Right, well spotted. This seems to fix the problem. Thanks!

1 Like
17

I would like to reopen this post as I think there might be a new issue with 3D triangles. Consider the following MWE failing with latest dolfinx nightly Docker image (pulled 13 May):

#!/usr/bin/env python3
from mpi4py import MPI
from dolfinx import fem, io
import ufl

comm = MPI.COMM_WORLD

meshname='triangle3D.xdmf'

with io.XDMFFile(comm, meshname, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    msh = infile.read_mesh(name="Grid")

V = fem.VectorFunctionSpace(msh, ("CG", 1))
u = fem.Function(V)

print(u.vector.getBlockSize())

gr = ufl.sym(ufl.grad(u))

The symmetric part operator yields ValueError: Cannot take symmetric part of rectangular matrix with dimensions (2, 3).

The block size of vector u is now 2, whereas I believe it should be 3.

The mesh file triangle3D.xdmf is


<Xdmf Version="3.0">
  <Domain>
    <Grid Name="Grid">
      <Geometry GeometryType="XYZ">
        <DataItem DataType="Float" Dimensions="3 3" Format="XML" Precision="8">
          0.0 0.0 0.0
          1.0 0.0 0.0
          0.0 1.0 0.0
        </DataItem>
      </Geometry>
      <Topology NodesPerElement="3" NumberOfElements="1" TopologyType="triangle">
        <DataItem DataType="Int" Dimensions="1 3" Format="XML" Precision="4">
          0 1 2
        </DataItem>
      </Topology>
    </Grid>
  </Domain>
</Xdmf>

Is there any new syntax/something I miss here?

Best,
Marc

18

Thanks for reporting this, see: Fix 3D manifolds from vector and tensorfunctionspace by jorgensd · Pull Request #2655 · FEniCS/dolfinx · GitHub

19

Hello,

I am trying to model 2D plate in 3D space (so that each node of the 2D plate mesh have 6 dofs).
I base my code on this tutorial, just simplifying the geometry to unit square Linear shell model — Numerical tours of continuum mechanics using FEniCS master documentation
I managed to have it working in dolfin (just had to regenerate the mesh in 3D, instead using 2D) and now, I am retrying it using dolfinx.
This is my code to generate 3D mesh:

Blockquote
import dolfinx
import ufl
from mpi4py import MPI

n_div = 2
mesh_2D = dolfinx.mesh.create_unit_square(
MPI.COMM_SELF, n_div, n_div, dolfinx.mesh.CellType.triangle
)
gdim, shape, degree = 3, “triangle”, 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(
ufl.MixedElement(
[
ufl.VectorElement(“Lagrange”, cell, degree),
ufl.VectorElement(“CR”, cell, degree),
]
)
)
cells = mesh_2D.geometry.dofmap
x = mesh_2D.geometry.x
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)

Ue = (
ufl.VectorElement(
“Lagrange”,
ufl.triangle,
2,
dim=3,
),
)
Te = ufl.VectorElement(“CR”, ufl.triangle, 1, dim=3)

V = dolfinx.fem.FunctionSpace(mesh, ufl.MixedElement([Ue, Te]))

v = dolfinx.fem.Function(V)
u, theta = ufl.split(v)

ufl.grad(u)

I receive the follwoing error:
Traceback (most recent call last):

  • File “/home/misko/fenicsx_start/to_del.py”, line 31, in *
  • V = dolfinx.fem.FunctionSpace(mesh, ufl.MixedElement([Ue, Te]))*
  •    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
  • File “/home/misko/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py”, line 581, in FunctionSpace*
  • return functionspace(mesh, element, form_compiler_options, jit_options)*
  •       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*
  • File “/home/misko/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py”, line 528, in functionspace*
  • raise ValueError(“Non-matching UFL cell and mesh cell shapes.”)*
    ValueError: Non-matching UFL cell and mesh cell shapes.

if I do not specify the geometrical_dimention for ufl.Cell, then I get the follwoing error:
Traceback (most recent call last):

  • File “/home/misko/fenicsx_start/to_del.py”, line 36, in *
  • ufl.grad(u)*
  • File “/home/misko/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/operators.py”, line 346, in grad*
  • return Grad(f)*
  •       ^^^^^^^*
  • File “/home/misko/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/differentiation.py”, line 260, in init*
  • self._dim = find_geometric_dimension(f)*
  •            ^^^^^^^^^^^^^^^^^^^^^^^^^^^*
  • File “/home/misko/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/domain.py”, line 396, in find_geometric_dimension*
  • raise ValueError(“Cannot determine geometric dimension from expression.”)*
    ValueError: Cannot determine geometric dimension from expression.

I would be grateful for any help.

20

Please format the code with 3x`, ie

```python
# add code here
```