Issue with unstructured hex mesh

[Paras Kumar]

Dear Deal.II community,

I am working on homogenization of particulate nano-composites and thus have to solve the finite deformation quasistatics problem on different microstructures, such as one shown in the attachment pic-1.png.

In order to automate the generation of (2D & 3D) meshes for microstructures comprising of randomly distributed filler particles of arbitrary radius, we employ a self-written python code using GMSH for mesh generation. The setup works well for 2D but for 3D, the resulting mesh (left) culminates in irregularities in the stress contour as compared to that obtained with a structured mesh (right) as can be seen from the attachment pic-2.png. Due to the random spatial distribution and arbitrary size of particles, unstructured mesh appeared as the feasible option.

In order to verify that the cause of such irregularities is the mesh and not our solver (based on deal.ii), we tested the mesh in Abaqus and similar irregularities were observed, cf. pic-3.png. Refining the mesh does not help either.

Has anyone else working with unstructured Hex meshes, observed similar issues? Any clues on how to deal with this issue or on other tools (preferably open source) which one could use for generating better meshes in 3D would be of great help.

I understand that this is probably not a question directly linked to deal.II and thus my apologies for the inconvenience caused.

Best regards,

Paras Kumar

[Wolfgang Bangerth]

Paras,

> I am working on homogenization of particulate nano-composites and thus have to
> solve the finite deformation quasistatics problem on different
> microstructures, such as one shown in the attachment pic-1.png.
>
> In order to automate the generation of (2D & 3D) meshes for microstructures
> comprising of randomly distributed filler particles of arbitrary radius, we
> employ a self-written python code using GMSH for mesh generation. The setup
> works well for 2D but for 3D, the resulting mesh (left) culminates in
> irregularities in the stress contour as compared to that obtained with a
> structured mesh (right) as can be seen from the attachment pic-2.png. Due to
> the random spatial distribution and arbitrary size of particles, unstructured
> mesh appeared as the feasible option.
>
> In order to verify that the cause of such irregularities is the mesh and not
> our solver (based on deal.ii), we tested the mesh in Abaqus and similar
> irregularities were observed, cf. pic-3.png. Refining the mesh does not help
> either.

That suggests that you are doing something wrong. I don't know what exactly
the problem is, and in particular how you assign material properties to cells,
but we know from theory that the finite element solution must converge to the
exact solution if you just refine the mesh often enough. If it doesn't for
you, then something is amiss.

This does not say anything about the *absolute level of accuracy*. I would not
be surprised if for a mesh like the one you show (in which pretty much every
hexahedron is poorly shaped), the solution is less than for the corresponding
tetrahedal mesh. But if you refine it a couple of times, you should get a more
accurate solution. Is this not the case?

> Has anyone else working with unstructured Hex meshes, observed similar issues?
> Any clues on how to deal with this issue or on other tools (preferably open
> source) which one could use for generating better meshes in 3D would be of
> great help.

The quality of meshes matters for the absolute level of error. The hex meshes
you get by subdividing tets are generally quite poor, though I know of people
who are using these routinely and report that they nevertheless get good
accuracy. A better approach is certainly to see if you can find ways to
*directly* create hex meshes. gmsh can do this to some degree. For the sphere
you show, deal.II can also generate a high-quality mesh itself.

Best
W.

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

[Paras Kumar]

Dear Prof. Bangerth,

Thank you for your response.

On Wed, Sep 9, 2020 at 1:08 AM Wolfgang Bangerth <bang...@colostate.edu> wrote:

Paras,

> I am working on homogenization of particulate nano-composites and thus have to
> solve the finite deformation quasistatics problem on different
> microstructures, such as one shown in the attachment pic-1.png.
>
> In order to automate the generation of (2D & 3D) meshes for microstructures
> comprising of randomly distributed filler particles of arbitrary radius, we
> employ a self-written python code using GMSH for mesh generation. The setup
> works well for 2D but for 3D, the resulting mesh (left) culminates in
> irregularities in the stress contour as compared to that obtained with a
> structured mesh (right) as can be seen from the attachment pic-2.png. Due to
> the random spatial distribution and arbitrary size of particles, unstructured
> mesh appeared as the feasible option.
>
> In order to verify that the cause of such irregularities is the mesh and not
> our solver (based on deal.ii), we tested the mesh in Abaqus and similar
> irregularities were observed, cf. pic-3.png. Refining the mesh does not help
> either.
 

That suggests that you are doing something wrong. I don't know what exactly
the problem is, and in particular how you assign material properties to cells,
but we know from theory that the finite element solution must converge to the
exact solution if you just refine the mesh often enough. If it doesn't for
you, then something is amiss.

I also had this doubt that there was some bug in my code and thus I tried with the Abaqus solver where not much could go wrong except the mesh.

With regards to material property assignment, it is done by means of *material option within the .inp mesh file (as was also used for the PR https://github.com/dealii/dealii/pull/9466.) and the material IDs appear to be correct when we export the mesh as vtk using deal.II functions and visualize it in Paraview.

This does not say anything about the *absolute level of accuracy*. I would not
be surprised if for a mesh like the one you show (in which pretty much every
hexahedron is poorly shaped), the solution is less than for the corresponding
tetrahedal mesh. But if you refine it a couple of times, you should get a more
accurate solution. Is this not the case?

We also tried with a refined mesh (resulting in 360K elements as compared to 120K elements for the mesh in pic-2.png), but the irregularities still remain as  can be seen in pic-4.png.Considering that fact that we wish to model multi-particle systems, such large number of elements would render such a mehsing approach computationally infeasible.

 

> Has anyone else working with unstructured Hex meshes, observed similar issues?
> Any clues on how to deal with this issue or on other tools (preferably open
> source) which one could use for generating better meshes in 3D would be of
> great help.

The quality of meshes matters for the absolute level of error. The hex meshes
you get by subdividing tets are generally quite poor, though I know of people
who are using these routinely and report that they nevertheless get good
accuracy. A better approach is certainly to see if you can find ways to
*directly* create hex meshes. gmsh can do this to some degree. For the sphere
you show, deal.II can also generate a high-quality mesh itself.

We have already tried several options available in GMSH, but it did not seem to work. Could you please suggest some tool for automatic "direct" Hex mesh generation. 

Is it not a good idea to use the tet-to-hex approach for generating Hex meshes, since this appears to be the only feasible option with GMSH considering the complexity of the geometry.

The actual geometry in the current case comprises a spherical particle of material-2 embedded in a cube of material-1 (later, it could be 50 such particles embedded in the matrix). I went through the grid generation functions of deal.II but could not find any possibilities of  doing an "intersection" between triangulations as would be necessary for the matrix material. Did I miss something?

Best regards,

Paras Kumar

[blais...@gmail.com]

Dear Paras,

I have faced a lot of similar issues in the past and maybe these suggestions can help (or maybe not).

1. Generally visual rendering software (at least paraview) renders data on quad/hex meshes using triangular polygons. Consequently, when you look at the color map within a square, paraview is actually rendering the square using two triangles. Consequently, the color map often appears wrong because visually you are not seing the bi-linear interpolation, but some sort of "cut in half quad" where each cut is rendered using a linear simplex element. Do you think there is another metric you can look at other than the visualisation? What you are seeing right now might in fact be an artifact of the GPU rendering of your solution.

2. One issue with sub-divided tets into hex is that the added node on the triangular face does not match the topology of your geometry. In CFD this can introduce very large errors, so I guess it does the same thing for solid mechanics. What you can do is generate that mesh in GMSH, then "inflate it" to match the manifold of the real geometry you want to simulate. This means displacing a few nodes. You can achieve this robustly by solving a linear elasticity equation with dirichlet boundary conditoin to ensure that your displacement deforms your mesh everywhere. I have found this to work very well when the objects are spherical.

3. If the two solutions above don't work, you can try alternative meshing software (like cubit). However, they are harder to automatize in a GMSH-like fashion.

Feel free to reach out if you need any help!

Best

Bruno

[Wolfgang Bangerth]

On 9/9/20 12:31 PM, Paras Kumar wrote:
>
> This does not say anything about the *absolute level of accuracy*. I would
> not
> be surprised if for a mesh like the one you show (in which pretty much every
> hexahedron is poorly shaped), the solution is less than for the corresponding
> tetrahedal mesh. But if you refine it a couple of times, you should get a
> more
> accurate solution. Is this not the case?
>
> We also tried with a refined mesh (resulting in 360K elements as compared to
> 120K elements for the mesh in pic-2.png), but the irregularities still remain
> as  can be seen in pic-4.png.Considering that fact that we wish to model
> multi-particle systems, such large number of elements would render such a
> mehsing approach computationally infeasible.

Right. But it helps illuminate where the problem may be.

From the look of it, the solution is at least smoother. That might suggest
that the issue is not a bug in the code, but that the solution really does
converge, though maybe slowly. Do you have the resources to do one more
refinement, just to see whether the solution continues to get better? That's
not going to be a production run, but a one-time computation that might take a
day or two if necessary.

> The quality of meshes matters for the absolute level of error. The hex meshes
> you get by subdividing tets are generally quite poor, though I know of people
> who are using these routinely and report that they nevertheless get good
> accuracy. A better approach is certainly to see if you can find ways to
> *directly* create hex meshes. gmsh can do this to some degree. For the sphere
> you show, deal.II can also generate a high-quality mesh itself.
>
> We have already tried several options available in GMSH, but it did not seem
> to work. Could you please suggest some tool for automatic "direct" Hex mesh
> generation.

I'm no expert in using gmsh, so if you say that you couldn't get it to work,
then I have no other suggestion either.

> Is it not a good idea to use the tet-to-hex approach for generating Hex
> meshes, since this appears to be the only feasible option with GMSH
> considering the complexity of the geometry.
> The actual geometry in the current case comprises a spherical particle of
> material-2 embedded in a cube of material-1 (later, it could be 50 such
> particles embedded in the matrix). I went through the grid generation
> functions of deal.II but could not find any possibilities of  doing an
> "intersection" between triangulations as would be necessary for the matrix
> material. Did I miss something?

If it's just one sphere in a ball, then take a look at the various functions
in namespace GridGenerator that build geometries with holes. If you take a
look at these functions and how they are implemented, you will probably
understand how you can generate a mesh for a single hole in a cube, and then
combine it with a mesh for the spherical hole via GridGenerator::hyper_sphere.

That won't work if you have many inclusions. We generally deal with many
inclusions by just using a regular mesh for the domain and at each quadrature
point asking whether it is inside an inclusion or in the background material.
Here is an example of how this could be done:
https://github.com/geodynamics/aspect/blob/master/benchmarks/nsinker/nsinker.cc
The value() returns the coefficient at a particular point. The mesh does not
actually to to mesh these inclusions, we just leave it to each quadrature
point. The benchmark this file implements is like this:
https://geodynamics.org/cig/news/research-highlights/november-2019/
You can find more information about this "nsinker" benchmark here:
https://tigerprints.clemson.edu/cgi/viewcontent.cgi?article=3528&context=all_dissertations

[Paras Kumar]

Dear Prof. Bangerth,

From the look of it, the solution is at least smoother. That might suggest
that the issue is not a bug in the code, but that the solution really does
converge, though maybe slowly. Do you have the resources to do one more
refinement, just to see whether the solution continues to get better? That's
not going to be a production run, but a one-time computation that might take a
day or two if necessary.
 

As can be seen in figures 7 & 8 of the attached report, a further refinement results in relatively smoother contours, but the irregularities are still not eliminated.

 

That won't work if you have many inclusions. We generally deal with many
inclusions by just using a regular mesh for the domain and at each quadrature
point asking whether it is inside an inclusion or in the background material.
Here is an example of how this could be done:
https://github.com/geodynamics/aspect/blob/master/benchmarks/nsinker/nsinker.cc
The value() returns the coefficient at a particular point. The mesh does not
actually to to mesh these inclusions, we just leave it to each quadrature
point. The benchmark this file implements is like this:
https://geodynamics.org/cig/news/research-highlights/november-2019/
You can find more information about this "nsinker" benchmark here:
https://tigerprints.clemson.edu/cgi/viewcontent.cgi?article=3528&context=all_dissertations 

This approach would probably not be applicable for our use-case since, we  need to model the interface between the filler and the matrix regions accurately.

Best,

Paras

[Paras Kumar]

Dear Bruno,

Thank you for the suggestions. Here are a few remarks/queries to the same.

1. Generally visual rendering software (at least paraview) renders data on quad/hex meshes using triangular polygons. Consequently, when you look at the color map within a square, paraview is actually rendering the square using two triangles. Consequently, the color map often appears wrong because visually you are not seing the bi-linear interpolation, but some sort of "cut in half quad" where each cut is rendered using a linear simplex element. Do you think there is another metric you can look at other than the visualisation? What you are seeing right now might in fact be an artifact of the GPU rendering of your solution.

Stresses and strains (both based on the derivative of the solution field and thus suffering from the irregularities) are the commonly chosen metrics for visualization. In order to verify if this is just a software induced artefact, I also tired with VISIT, but got the same problem there as well. As an alternative to nodal stress values (computed by L2 projection of the QP stress values), I also tried to plot the cell-wise average stresses but in that case, the contour was even more irregular. Is there any visualization tool, which can render date directly using quads/hexes?  Further, it is noteworthy that the solution i.e. the displacement, does not show such irregularities.

 

2. One issue with sub-divided tets into hex is that the added node on the triangular face does not match the topology of your geometry. In CFD this can introduce very large errors, so I guess it does the same thing for solid mechanics. What you can do is generate that mesh in GMSH, then "inflate it" to match the manifold of the real geometry you want to simulate. This means displacing a few nodes. You can achieve this robustly by solving a linear elasticity equation with dirichlet boundary conditoin to ensure that your displacement deforms your mesh everywhere. I have found this to work very well when the objects are spherical.

In order to verify if the "improper" modeling of the curved particle surface due to the reason you explained above, is the cause for this irregularities, we tried with a cubic particle, thereby eliminating the above issue, but still observe irregularities in the stress contour on the particle surface, cf. Figures 3-5 of the attachment.  With this observation, we could probably conclude that improper modeling of the curved geometry is not the only reason for the irregularities.

 

3. If the two solutions above don't work, you can try alternative meshing software (like cubit). However, they are harder to automatize in a GMSH-like fashion.

Trelis [https://coreform.com/products/trelis/]  seems to be the next thing to try out. Could you please comment on in what sense would it be "harder to automate". The approach now would be to generate the RVE geometry within our python code in .step format, and then call Trelis from within our code for mesh generation.  

Since, we plan to make the RVE generation tool open source at some point in the future, an open source meshing backend would have been great. Are you aware of any open source alternatives which could probably solve this problem?

Thanks again for your help and best regards,

Paras

 

[Wolfgang Bangerth]

On 9/22/20 10:41 AM, Paras Kumar wrote:
> As can be seen in figures 7 & 8 of the attached report, a further refinement
> results in relatively smoother contours, but the irregularities are still not
> eliminated.

OK, that's good to know at least as a check that your program is doing
something right.

> This approach would probably not be applicable for our use-case since, we
> need to model the interface between the filler and the matrix regions accurately.

Then I'm afraid I'm out of ideas short of using a program that can actually
create hex meshes :-(

[Wolfgang Bangerth]

On 9/22/20 11:09 AM, Paras Kumar wrote:
> In order to verify if the "improper" modeling of the curved particle surface
> due to the reason you explained above, is the cause for this irregularities,
> we tried with a cubic particle, thereby eliminating the above issue, but still
> observe irregularities in the stress contour on the particle surface, cf.
> Figures 3-5 of the attachment.  With this observation, we could probably
> conclude that improper modeling of the curved geometry is not the only reason
> for the irregularities.

I know this isn't helpful from the perspective of actually solving your
problem, but: Excellent work testing out different hypotheses and validating
every step of your work. What you show here is a model of how one builds
numerical solutions: If there is something that seems odd, one systematically
simplifies the simulation (sphere -> cube), and then cross checks different
approaches (tets vs hexes vs uniform discretization) until it's clear what the
issue is. I wished every student project were as comprehensively explored! :-)

[Wolfgang Bangerth]

On 9/22/20 4:27 PM, Wolfgang Bangerth wrote:
> I wished every student project were as comprehensively explored! :-)

I tacitly assumed that Maurice Rohracker is your student. I didn't mean to
imply that the work you're doing is at the "student level" -- quite the
contrary: What your reports show is just really impressive!

[Timo Heister]

I am coming late to the conversation (and I only looked at it quickly,
sorry), but I have two additional suggestions:
1. Depending on the deformation of cells and the PDE in question, it
might be important to check if the quadrature degree is high enough to
not contribute to the overall error. This is easy to check by
increasing the quadrature degree in all places but you can of course
go back to your PDE and check what you need.
2. Did you play with linear solver tolerances to see if this matters?

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/5f6253d8-11a6-a530-ab57-5e6f020a650f%40colostate.edu.



--
Timo Heister
http://www.math.clemson.edu/~heister/

[Paras Kumar]

Dear Prof. Bangerth,

Thank you for the appreciation.

In the quest of alternative solutions, I came across the Mesquite library and was curious to know if that could be helpful to improve the quality of the Delaunay based Hex mesh so as to avoid the issue with stress irregularities.

It also seems that Mequite's integration within deal.II is under progress[https://github.com/dealii/dealii/pull/5971]. \

Since, I have no prior experience with this library, some guidance wrt the current issue would be of great help.

Best,

Paras

[Wolfgang Bangerth]

On 9/25/20 4:58 AM, Paras Kumar wrote:
>
> In the quest of alternative solutions, I came across the Mesquite library and
> was curious to know if that could be helpful to improve the quality of the
> Delaunay based Hex mesh so as to avoid the issue with stress irregularities.
> It also seems that Mequite's integration within deal.II is under
> progress[https://github.com/dealii/dealii/pull/5971]. \
> Since, I have no prior experience with this library, some guidance wrt the
> current issue would be of great help.

You could try that patch (though I don't know the current status of it -- the
issue is that Mesquite is no longer maintained).

But I suspect that it will not help very much. You can understand this in 2d:
Even if you have a perfect triangular mesh and then subdivide each triangle
into 3 quadrilaterals, you end up with quadrilaterals shows shape can not be
improved by moving around nodes. Furthermore, and that seems to be the key
issue for theoretical reasons too complicated to discuss here, you really want
2d meshes where at the majority of nodes 4 quadrilaterals come together. But
the quadrilateral mesh you get out of subdividing triangles has nodes where
either 3, 4, or 6 quadrilaterals come together (corresponding to triangle
centers, triangle edge midpoints, and triangle vertices). You end up with much
higher errors on these kinds of meshes, and I think that is what you are
seeing. Libraries like Mesquite optimize the mesh by moving around nodes, and
that has no effect on the valence of nodes.

(I'm sure Mesquite can also do other operations, but they are far more
complicated on quad and hex meshes and I don't think the interface enabled that.)

[Jean-Paul Pelteret]

Since you have access to Cubit/Trellis, “that patch” is probably not necessary. Cubit uses the Mesquite library for mesh quality improvement, so if you could try to increase the mesh quality that way. Although I’d have to agree with Wolfgang that it might not help too much in this case, but that’s an educated guess.

Best,
Jean-Paul

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f4460ec5-52b8-243c-3e34-63cc4171454c%40colostate.edu.

[Jean-Paul Pelteret]

Oh, and and just in case you’ve not heard of it before and it might be useful, there’s also this tool that does automated hex meshing:

https://www.openfoam.com/documentation/guides/latest/doc/guide-meshing-snappyhexmesh.html

I’ve never used it though, and I have no idea how you would get the mesh in and out of OpenFOAM, and whether one can prevent the creation of meshes with hanging nodes. Nevertheless, maybe its worth looking into.

Best,

Jean-Paul

[blais...@gmail.com]

Be careful about SnappyHexMesh. You can use it to make hex-only meshes, but it will generate a tesselated mesh (which does not really conform to the boundary, it is by all means a voxel mesh).

If you fully use Snappy, it will generate a polyhedral mesh with as many as 20 faces per element. This is all fine and dandy for FVM (like OpenFOAM), but it will not work within deal.II.