3D Mixed Simplex Mesh

[Alex Quinlan]

Hello everyone,

I have been reading the messages in this user group for a while, but this is my first post!  I am hoping that you can help with my problem. 

I am coming from a background of solid mechanics FEA and I've been using CalculiX as an openSource program for some time.  My group has decided to switch to deal.ii to implement some advanced functionality that the other programs don't have.  I'm an engineer with some programming experience, though fairly new to C++.

We are working with some complex geometries that are not practical to mesh with only hex elements.  The meshes are generated externally, so I am trying to import the meshes into deal.ii.  I have tried to do this in two general ways and run into problems with each, so I will split my post into 2 parts.

Part 1:  Abaqus .inp format

I had success importing hex meshes from the abaqus <.inp> format, but it seems that the UCD reader does not accept tet meshes.  I am basing this on the following bit of code in grid_in::read_ucd around Line 948

      if (((cell_type == "line") && (dim == 1)) ||
          ((cell_type == "quad") && (dim == 2)) ||
          ((cell_type == "hex") && (dim == 3)))
        // found a cell
        {

The absence of

       ((cell_type == "tet") && (dim == 3))

suggests that simplex elements are not supported when importing from .inp.  Am I correct in this conclusion.  How difficult would it be for me to add this functionality in?

Part 2:  GMSH format

While looking at grid_in.cc, I saw that the MSH reader supported tet elements.  I found the FAQ page supplied these examples of simplex meshes: https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html

I first ran the 2D simplex mesh example, and then I modified it to work for 3D without much problem.  Then, I switched to the mixed mesh, and had success running the 2D example.  However, I started having problems when I switched to 3D.  The mesh seems to import, but then I have issues operating with it. I found several places within fe_simplex_p.cc (L714, L759) where the dimension is asserted to be 2:

       AssertDimension(dim, 2);

 

This seems like a problem for a 3D analysis, so I am wondering if I need to do something different so that this function doesn't get called.

 

The problem seems to arrive when I call  dof_handler.distribute_dofs(fe);  I've attached my modified version of the example program along with my simple mixed-element 3D mesh.

Is there something that I am missing regarding the dof_handler when trying to use a 3D mixed-element mesh?

Thanks very much,

Alex

[Alex Quinlan]

I realized that the mesh file I provided had an issue.  The corner nodes of the tetrahedral elements coincided with a mid-face node of the hexahedron element.  I've fixed that and uploaded a new mesh.

Thanks,

Alex

[Wolfgang Bangerth]

Alex:

> I am coming from a background of solid mechanics FEA and I've been using
> CalculiX as an openSource program for some time.  My group has decided
> to switch to deal.ii

Good choice -- we applaud you :-)

> Part 1:  Abaqus .inp format
> I had success importing hex meshes from the abaqus <.inp> format, but it
> seems that the UCD reader does not accept tet meshes.  I am basing this
> on the following bit of code in grid_in::read_ucd around Line 948
>
>       if (((cell_type == "line") && (dim == 1)) ||
>           ((cell_type == "quad") && (dim == 2)) ||
>           ((cell_type == "hex") && (dim == 3)))
>         // found a cell
>         {
>
> The absence of
>        ((cell_type == "tet") && (dim == 3))
> suggests that simplex elements are not supported when importing from
> .inp.  Am I correct in this conclusion.  How difficult would it be for
> me to add this functionality in?

Not very difficult in all likelihood. The key difficulty in adding
support for other cells in this function is that it uses the
GeometryInfo<dim> class, which only provides information about hypercube
cells. For example, we read vertex indices like this in that function:
for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
in >> cells.back().vertices[GeometryInfo<dim>::ucd_to_deal[i]];

This needs to be changed to account for the number of vertices a cell
has, given the cell type, not the dimension. The right approach will be
to use the ReferenceCell objects defined in namespace ReferenceCells
(plural), and the other readers in that function will serve as guidance.
You might want to take a look at the gmsh reader, for example.

> Part 2:  GMSH format
> While looking at grid_in.cc, I saw that the MSH reader supported tet
> elements.  I found the FAQ page supplied these examples of simplex
> meshes: https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html
>
> I first ran the 2D simplex mesh example, and then I modified it to work
> for 3D without much problem.  Then, I switched to the mixed mesh, and
> had success running the 2D example.  However, I started having problems
> when I switched to 3D.  The mesh seems to import, but then I have issues
> operating with it. I found several places within fe_simplex_p.cc (L714,
> L759) where the dimension is asserted to be 2:
>
>        AssertDimension(dim, 2);
> This seems like a problem for a 3D analysis, so I am wondering if I need
> to do something different so that this function doesn't get called.

No, I think this is just not implemented. The function is called on
mixed meshes when the library wants to know how elements of type
FE_SimplexP on tetrahedra connect to another element types on other cell
shapes (namely, on wedges and pyramids). The comparison is simply not
implemented, but if you want to compute on mixed 3d meshes, you will
have to implement the missing parts of this function. As long as you
want to use the same polynomial degree on all cells, this should not
even be particularly difficult and I would be happy to walk you through
the necessary steps if you were interested in writing up a patch!

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Alex Quinlan]

Dear Wolfgang,

Thanks very much for your reply.  I'm happy to hear that you think my problems can be resolved without too much difficulty. I'd like to try writing a patch for the mixed 3D meshes.  I have done some development work before, but I have never sent it out of our group; your advice on writing the steps to write a patch would be most welcome!

I will re-review the relevant code and see if I can put together a general plan on how to approach this.

Best regards,

Alex

[Wolfgang Bangerth]

On 11/11/22 07:06, Alex Quinlan wrote:
>
> Thanks very much for your reply.  I'm happy to hear that you think my problems
> can be resolved without too much difficulty. I'd like to try writing a patch
> for the mixed 3D meshes.  I have done some development work before, but I have
> never sent it out of our group; your advice on writing the steps to write a
> patch would be most welcome!

Alex:
At the most basic level, the steps to contribute patches is outlined here:
https://dealii.org/participate.html
Let me know if you have questions and I'll be happy to walk you through
whatever is missing in these documents.

> I will re-review the relevant code and see if I can put together a general
> plan on how to approach this.

And here too: Feel free to discuss it with me or here on this list.

For the specific code where you ended up with an exception: This is the place
where a finite element describes how the degrees of freedom on one of its
faces should interact with the degrees of freedom of the finite element on a
neighboring cell. For example, if you have a P1 and a P2 element on
neighboring cells (let's assume triangles for the moment), then the mid-edge
degree of freedom of the P2 element needs to have a specific value to ensure
that the resulting finite element field is continuous. In your situation, the
function has simply not been taught what these constraints should be if you
have a FE_SimplexP(1) element on a tetrahedron's triangular face meet the
degrees of freedom of a FE_WedgeP(1) or FE_PyramidP(1) on that common
triangular face. The answer is that the degrees of freedom should all be equal
to each other, as in some of the other cases handled by the function.