FEniCS shells mesh

Hi FEniCS family,

I have few queries about FEniCS shells.

  1. Importing shell mesh to FEniCS. I have nodes and element connectivity. I tried it importing using gmsh (converted .cgns file to .msh ). On converting to .xml using dolfin-convert, I am getting error " Unable to find cells of supported types".
    gmsh mesh_shell

Can you provide a convenient way of importing shell mesh?

  1. I want to use my own basis functions. How can I add it to UFL library so that installed FEniCS supports it? Can I define a UFL function in the code itself for new basis functions.

(My requirement was to use curvilinear coordinates (xi,eta) to interpolate the geometry as used by commercial software? Is there any existing basis function satisfying this requirement? The basis functions are for 8 noded shell element).

Any help is greatly appreciated.

dolfin-convert is outdated. Use meshio for conversion if you have to stick with legacy FEniCS.

If you want to use your own basis functions, I would strongly suggest that you use FEniCSx, as you can define your own finite elements there, see for instance @mscroggs brilliant demo: Creating TNT elements using Basix’s custom element interface — DOLFINx 0.7.2 documentation

2 Likes

Thanks @dokken for your response.

I would try to switch to meshio. Can you recommend the specific examples for meshio? Does it treat shell element as a quad element or any other specific element?

For FEniCSx, is there any reference guide other than available demos to follow changed basic UFL, Dolfin syntax like in legacy FENiCS (there was an entire available book) for getting better clarity. Any suggestion is helpful.

Legacy fenics doesn’t have proper support for quad meshes, so I would really advice to stay away from it.

I’ve created several guides on FEniCSx:
https://jsdokken.com/dolfinx-tutorial/
https://jsdokken.com/FEniCS23-tutorial/README.html
and also the actual demos of DOLFINx
https://docs.fenicsproject.org/dolfinx/v0.7.2/python/demos.html

1 Like

Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken (jsdokken.com)

In your tutorial, you define that for using meshio, we need to remove the z coordinate. But, my mesh is defined in 3D space using quad elements shown in above figure. I am not able to figure out with meshio.

$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
40
      1     -1.76776695           -0.25     -1.76776695
      2              0.           -0.25            -2.5
      3      1.76776695           -0.25     -1.76776695
      4             2.5           -0.25              0.
      5      1.76776695           -0.25      1.76776695
      6              0.           -0.25             2.5
      7     -1.76776695           -0.25      1.76776695
      8            -2.5           -0.25              0.
      9     -1.76776695            0.25     -1.76776695
     10              0.            0.25            -2.5
     11      1.76776695            0.25     -1.76776695
     12             2.5            0.25              0.
     13      1.76776695            0.25      1.76776695
     14              0.            0.25             2.5
     15     -1.76776695            0.25      1.76776695
     16            -2.5            0.25              0.
     17     -2.30969882           -0.25     -0.95670861
     18     -1.76776695              0.     -1.76776695
     19     -2.30969882            0.25     -0.95670861
     20            -2.5              0.              0.
     21     -0.95670861           -0.25     -2.30969882
     22              0.              0.            -2.5
     23     -0.95670861            0.25     -2.30969882
     24      0.95670861           -0.25     -2.30969882
     25      1.76776695              0.     -1.76776695
     26      0.95670861            0.25     -2.30969882
     27      2.30969882           -0.25     -0.95670861
     28             2.5              0.              0.
     29      2.30969882            0.25     -0.95670861
     30      2.30969882           -0.25      0.95670861
     31      1.76776695              0.      1.76776695
     32      2.30969882            0.25      0.95670861
     33      0.95670861           -0.25      2.30969882
     34              0.              0.             2.5
     35      0.95670861            0.25      2.30969882
     36     -0.95670861           -0.25      2.30969882
     37     -1.76776695              0.      1.76776695
     38     -0.95670861            0.25      2.30969882
     39     -2.30969882           -0.25      0.95670861
     40     -2.30969882            0.25      0.95670861
$EndNodes
$Elements
8
1     3 2 0 1	8     1     9    16    
2     3 2 0 1   1     2    10     9   
3     3 2 0 1 	2     3    11    10    
4     3 2 0 1	3     4    12    11    
5     3 2 0 1	4     5    13    12    
6     3 2 0 1	5     6    14    13    
7     3 2 0 1   6     7    15    14    
8     3 2 0 1	7     8    16    15   
$EndElements

Using the above .msh file I tried meshio for the above quad shell element in 3D space. It is reading the shell mesh but, I need to ask how to proceed to get subdomains and mesh used by PBC and functionspace respectively.

import meshio
import numpy as np

# Read in mesh
msh = meshio.read("a.msh")

# Helper function to extract the cells and physical data
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

print(np.unique(msh.cell_data_dict['gmsh:physical']['quad']))
#print(np.unique(msh.cell_data_dict['gmsh:physical']['hexahedron']))

# Convert mesh to XDMF
quad_mesh = create_mesh(msh, "quad", prune_z=False)
#meshio.write("3Dmesh/3Dblock_facet_mesh.xdmf", quad_mesh)

#hexahedron_mesh = create_mesh(msh, "hexahedron", prune_z=False)
#meshio.write("3Dmesh/3Dblock_mesh.xdmf", hexahedron_mesh)
quad_mesh
<meshio mesh object>
  Number of points: 40
  Number of cells:
    quad: 8
  Cell data: name_to_read

My query is: what is the meshio alternative to get mesh and physical space.

mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")

Legacy dolfin does not support quadrilateral elements read from file. You would have to use DOLFINx for this.

I do not understand the last bit of your question. Reading in cell markers or facet markers has to be done by adding cell data to the corresponding meshes, as done in the tutorials.

Sorry, I was unclear in describing.
1.

after getting quad_mesh, how can I use this mesh data to create a Function Space.

Q = FunctionSpace(quad_mesh, ("CG", 1))
AttributeError                            Traceback (most recent call last)
Cell In[26], line 1
----> 1 Q = FunctionSpace(quad_mesh, ("CG", 1))

File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py:581, in FunctionSpace(mesh, element, form_compiler_options, jit_options)
    561 def FunctionSpace(mesh: Mesh,
    562                   element: typing.Union[ufl.FiniteElementBase, ElementMetaData,
    563                                         typing.Tuple[str, int, typing.Tuple, bool]],
    564                   form_compiler_options: typing.Optional[dict[str, typing.Any]] = None,
    565                   jit_options: typing.Optional[dict[str, typing.Any]] = None) -> FunctionSpaceBase:
    566     """Create a finite element function space.
    567 
    568     .. deprecated:: 0.7
   (...)
    579 
    580     """
--> 581     return functionspace(mesh, element, form_compiler_options, jit_options)

File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py:521, in functionspace(mesh, element, form_compiler_options, jit_options)
    519 try:
    520     e = ElementMetaData(*element)
--> 521     ufl_e = basix.ufl.element(e.family, mesh.basix_cell(), e.degree, shape=e.shape,
    522                               symmetry=e.symmetry, gdim=mesh.ufl_cell().geometric_dimension())
    523 except TypeError:
    524     ufl_e = element  # type: ignore

AttributeError: 'Mesh' object has no attribute 'basix_cell'

(there is no example of FEniCSx where you used external imported mesh and solved the variational problem). Can you give MWE for imported mesh to create FE formulation?
2.

# Convert mesh to XDMF
quad_mesh = create_mesh(msh, "quad", prune_z=False)
meshio.write("a_phy.xdmf", quad_mesh)
WARNING:py.warnings:/home/bagla0/.local/lib/python3.10/site-packages/h5py/__init__.py:36: UserWarning: h5py is running against HDF5 1.10.7 when it was built against 1.14.2, this may cause problems
  _warn(("h5py is running against HDF5 {0} when it was built against {1}, "

And with this the kernel dies everytime.

In legacy-fenics, we use the below code (MeshFunction) to get different subdomain using MeshFunction. What is MeshFunction alternative in fenicsx to define different subdomains?

subdomains = MeshFunction("size_t", quad_mesh, fname + "_physical_region.xml")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[27], line 1
----> 1 subdomains = MeshFunction("size_t", quad_mesh, fname + "_physical_region.xml")

NameError: name 'MeshFunction' is not defined

I am new to FEniCSx and facing issues to run a simple problem using imported quad mesh. Any example with imported mesh would be helpful.

https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html?highlight=gmsh#interfacing-with-gmsh-in-dolfinx
and Defining subdomains for different materials — FEniCSx tutorial
are two of many examples using external meshes.

This is related to how you have installed h5py. See for instance Importing mesh gives error - #8 by dokken

For final question about meshfunctions, see the demos ive linked to earlier in this post

1 Like