Are shell quad element T-junction topologies unsupported for dolfinx meshes?

I am hoping to use FEniCS-Shells to perform analysis on a somewhat complex shell geometry. I also hope to use quad elements if possible and, based on the discussion here, I opted to use dolfinx as it seems quad elements will be better supported. I installed dolfinx via Docker using simply:

docker run dolfinx/dolfinx

To demonstrate my issue, I made a simple T-junction mesh (I would attach the .xdmf file but the file extension seems not to be allowed…?):

I am trying to load the file with the code:

import dolfinx
from dolfinx import io
from dolfinx.io import XDMFFile
from mpi4py import MPI

infile = XDMFFile(MPI.COMM_WORLD, "tTopology.xdmf", "r", encoding=dolfinx.cpp.io.XDMFFile.Encoding.ASCII)
mesh = infile.read_mesh(name="Grid")
infile.close()

I then get the following error:

2021-03-23 15:08:14.799 ( 0.291s) [main ] graphbuild.cpp:93 ERR| Found three identical facets in mesh (local)
Traceback (most recent call last):
File “xdmfToPvd.py”, line 11, in
mesh = infile.read_mesh(name=“Grid”)
File “/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/io.py”, line 59, in read_mesh
mesh = cpp.mesh.create_mesh(self.comm(), cpp.graph.AdjacencyList_int64(cells),
RuntimeError: Inconsistent mesh data in GraphBuilder: three identical facets

The “three identical facets” comment seems to indicate that such a topology is not currently supported. Is there a particular reason for this, or might it be supported in the future? If it is unsupported, does anybody have suggestions for how this issue might be overcome?

Hi,

I don’t think fenics-shells will allow you to do T-junctions since their approach is based on an explicit mapping to a reference 2D surface. I am working on an alternative formulation which is not based on such a direct mapping but uses directly the non-manifold mesh pretty much like this demo for spatial beams. It is almost ready but written in old dolfin at the moment. It does not use quads but very simple locking-free triangular plate elements such as those discussed here.
What kind of analysis are you planning to do ? Small displacements in linear elasticity or finite strains ?

Thanks for the response! Having not yet dug into the details of fenics-shells it would not surprise me if there are various aspects that would have to be modified/extended for my purposes (e.g. a new mapping procedure). I fully expect to have to do some work in that regard and will work through those details as needed. Ultimately I hope to do a variety of analyses including linear and nonlinear/large deformation elasticity. I appreciate you linking some resources, I will take a closer look at those!

All of that said, these high-level details are moot if I FEniCS inherently cannot handle the types of mesh topologies I am dealing. Note in the example above that I have not done anything with fenics-shells…I am simply trying to read my mesh into the program. Resolving that issue is my immediate concern here.

I don’t know what is the status of non-manifold meshes in dolfinx, maybe you should file an issue. What’s sure is that it can be handled for continuous interpolations in legacy dolfin. However, quad elements are not well supported in legacy dolfin.

1 Like

Makes sense–I can and will file an issue as well. Helpful to know that such topologies can theoretically be handled in legacy dolfin…I suppose adopting triangular elements could be an alternative, albeit less attractive, solution in this case. Thanks for your input!

This test was added recently, to detect meshing errors. It could easily be disabled, if the mesh s known to be OK.

1 Like

Hi aherrema

For your mesh question:- if you want to use old dolfin then you can use a mixed dimensional mesh with MeshView or a 2d mesh in a 3d space. I personally wouldn’t bother with quad meshes, get it working with tris first.
Here is some sample code I used to make the mesh for my project - it won’t work directly for you because it relies on other functions as well but you can get the idea from it
Note the lines which create a VectorFunctionSpace and dolfin.File(self.Name+’.xml.gz’) <<self._dolfin_mesh are important - they are doing nothing but for some reason the code would crash without them

def dolfin_meshes(self, *, meshType='mixedDimensional',gdim=3): 
    cellTypes={1:df.CellType.Type.interval.__repr__().split('.')[-1] , 2: df.CellType.Type.triangle.__repr__().split('.')[-1] ,3:df.CellType.Type.tetrahedron.__repr__().split('.')[-1]}
    
    if 'separate' in meshType.lower():
        ret=[None]*3
        for dim in [1, 2, 3]:
            elts, nodes=self.elts_and_nodes_of_dim(dim)
            if len(elts):
                model=self.tri_model(only2d=1) if dim==2 else self 
                gdim1= dim if 'dims' in  meshType.lower() else gdim
                mesh = df.Mesh()
                mesh_editor = df.MeshEditor()        
                mesh_editor.open(mesh,cellTypes[dim],   dim, gdim1) 
                mesh_editor.init_vertices(len(nodes))   
                mesh_editor.init_cells(len(elts ))#
                nodedict={nd.Id:i for i, nd in enumerate(nodes)}
                for i, nd in enumerate(nodes):mesh_editor.add_vertex(i,nd.x[:gdim1])          
                for i, el in enumerate(elts): mesh_editor.add_cell(i,numpy.array([nodedict[n.Id]  for n in el.nodes], dtype=numpy.uintp)) 
                mesh_editor.close()
                VectorSpace1= df.VectorFunctionSpace(mesh, 'Lagrange', 1, dim=gdim)#for some reason this seems needed to make fenics initiate the facets 
                ret[dim-1]=mesh
        return ret
    if 'mixed' in meshType.lower():
        if self._dolfin_mesh==None:
            self._dolfin_mesh = df.Mesh()
            mesh_editor = df.MeshEditor()      
            gdim=3
            cellTypes={1:df.CellType.Type.interval.__repr__().split('.')[-1] , 2: df.CellType.Type.triangle.__repr__().split('.')[-1] ,3:df.CellType.Type.tetrahedron.__repr__().split('.')[-1]}
            dim=2 if len(self.shell_elts()[0]) else 1 #TODO deal with adding different dims
            model=self
            maxdim=3-[len(model.elts_of_dim(dim)[0])>0 for dim in [1, 2, 3]][::-1].index(1)
            if maxdim==2:
                elts=self.tri_model(only2d=1).elts_of_dim(maxdim)[0] 
            else:
                elts=model.elts_of_dim(maxdim)[0]
            mesh_editor.open(self._dolfin_mesh,cellTypes[maxdim],   maxdim, gdim) 
            mesh_editor.init_vertices(len(model.nodes))   
            mesh_editor.init_cells(len(elts ))#TODO this is probably wrong for mixed dim meshes
            for i, nd in enumerate(model.nodes):mesh_editor.add_vertex(i,nd.x)
            if len(elts):
                for i, el in  enumerate(elts): 
                    mesh_editor.add_cell(i,numpy.array([model.nodedict[n.Id]  for n in el.nodes], dtype=numpy.uintp)) 
            mesh_editor.close()
            dolfin.File(self.Name+'.xml.gz') <<self._dolfin_mesh 
            self._dolfin_meshes =[None]*3
            for dim in [1, 2, 3]:
                if dim==maxdim:
                    self._dolfin_meshes[dim-1]=self._dolfin_mesh
                elif len(self.elts_of_dim(dim)[0]) :
                    marker =dolfin.MeshFunction("size_t", self._dolfin_mesh, dim, 0) # mesh.topology().dim() - 1, 0)
                    for f in dolfin.facets(self._dolfin_mesh): marker[f] = 1
                    self._dolfin_meshes[dim-1] = dolfin.MeshView.create(marker, 1)#dim)
                    
        return self._dolfin_meshes    

Jermey Bleyer is right about fenics shells (its 3D example is a projection to 2D) but you could extend fenics shells with the same approach Jeremy used for beams, I did that a couple of years ago (i.e. use ideas from Elastic 3D beam structures — Numerical tours of continuum mechanics using FEniCS master documentation to update the fenics shells example. It worked in principle, even for large deformations, but it had some strange issues when the node order of adjacent elements changed , my maths was not strong enough to resolve it.

Using Jeremy’s new shell code will probably be a better starting point - his 3D beam code is well laid out, its easy to adapt his 3D beam example to work an arbitrary beam mesh, I guess his new shell code will also be clear. You could probably add the fenics shell implementations to it later. If you do add the fenics shells then I would recommend to start with the PSRI element, it performs well and is fairly easy to understand, the MITC elements are a bit more complex.

There is some IGA based code from David Kamensky which handles 3d shells - tIGAr , I haven’t used it myself but the code looks good and that also might be useful to you.

Thanks for the helpful tips, Thomas. On the dolfinx side, as previously suggested, I submitted an issue and am hopeful that the issue can be resolved there:

If I do decide to go with legacy dolfin + tris I will certainly take the approach you suggest into consideration. I am also familiar with David Kamensky’s work with tIGAr and may use that as a reference in the future.

Hi! I am one of the developers of FEniCS-Shells and also DOLFINx. Everything said above by Thomas and Jeremy said above is entirely correct, particularly the limitation regarding geometry.

Corrado Maurini and I are planning on porting FEniCS-Shells to DOLFINx. This will include the existing functionality (e.g. the existing models, MITC elements, reduced integration elements) and some new features. In particular, we plan to make use tangential differential calculus (TDC) for handling general discrete geometry (including junctions!) and various new finite element methods based on Regge elements for solving shear and membrane locking issues.

Thanks @jackhale for confirming the points above and for providing some more details about your plans for FEniCS-Shells in DOLFINx! It sounds like you are moving in a direction that will be beneficial both in general and for me specifically. I’m sure it will take time, but I will be sure to keep an eye out for this!

@bleyerj Hi Jeremy, I found this thread and wonder if you ever published your alternative shell element formulation for old dolfin, as it would be interesting to see? If not did you maybe encounter any specific issues with the implementation?

Hi, it just got published here:
https://comet-fenics.readthedocs.io/en/latest/demo/linear_shell/linear_shell.html

1 Like

Thank you, all your examples are really great and help a lot :+1:

Hi the reference Kirchhoff defection corresponds to clamped boundary conditions.
Here you enforce homogeneous BC on V.sub(0) which contains displacement only, you are thus modelling simple supports which should explain the approximate factor 4 that you have.

Hi, yes there was a typo remaining in the demo. The constitutive law should read

N = thick * plane_stress_elasticity(eps)
M = thick ** 3 / 12 * plane_stress_elasticity(kappa)
Q = mu * thick * gamma

Now it should work.
It has been corrected online, thanks for pointing this out.

1 Like