Netgen meshing with PrePoMax

[wilzzer]

Hi all !
I’m trying to understand how PrePoMax does its meshing and replicate it with a python program. For this, i’m using the python solutions provided by NGSolve to launch Netgen Mesher. For each shell object, the simple part of the code looks like this :

from netgen.occ import *
from netgen.meshing import *
geo = OCCGeometry(stp_file) #stp_file being a file name as str
export_meshRefinement(freecad_shape_obj, ref, refine_val, MESHSIZEFILE)  #Homemade function that creates MESHSIZEFILE. This file has the exact same structure as *meshRefinement* file of PrePoMax (as I inspired myself from it). I compared their output results and they are indeed exactly the same.
param = MeshingParameters(curvaturesafety=2, segmentsperedge=2, maxh=1, meshsizefilename=MESHSIZEFILE)
mesh = geo.GenerateMesh(param)
export_mesh(mesh, "mesh.inp") #Homemade function to export a netgen mesh to abaqus .inp format, which can be interpreted by both FreeCAD and PrePoMax.

This small snippet of code creates a netgen mesh, with refinement along edges defined in ref. I don’t touch the other meshing parameters as their default values are seemingly the same as in PrePoMax.
I then input the resulting .inp mesh in PrePoMax to visualize the meshing and compare FEM results.
First of all, we can see that the mesh, even though really similar, are not perfectly the same. I thought that despite this, I could get similar results, as the differences are so small. Unfortunately, that’s not the case ! With the same simulation setup, the homemade mesh crashes whereas the PrePoMax mesh works perfectly fine.

I have found that the parameters used by PrePoMax are stored in the temp file meshParameters. For netgen, i found a list of editable parameters here, in the function Ng_SetMeshingParameters. Unfortunately, most parameters don’t have the same designation in PrePoMax’s file and in netgen’s git.
So the question is the following : what are the netgen parameters defined by PrePoMax ? is it even possible to recreate the same mesh through netgen in python or am I going for a dead end ?
As a new user, i cannot attach files, but here’s a swisstransfer link : https://www.swisstransfer.com/d/278bf3d2-ebf3-4906-ad43-25b7d94b625d. You’ll find the two .pmx files, the .inp of my homemade generated mesh and the initial .stp file.

Any hint/idea/help will be very much appreciated !

[synt]

hi,

as i know, NetGen library used by PrePoMax still an old one’s version 5.x thus probably make it lead to a little bit of different result with the latest. Another possible reason is in the mesh size file reference and values.

to get exactly the same result, the library and executable of NetGen mesher can be call by sub-process in Python script. Alternatively, Gmsh library with Python wrapper can be used, it’s more convenient way since mesh size is defined in script and further calculated internally not by reference file explicitly as in NetGen.

[wilzzer]

Hi Synt,

Thank you for your answer.

Unfortunately no, the mesh size file follow exactly the same logic in structure : a series of lines defined by a start and an end point’s coordinates, followed by a value indicating the refinement. I made sure as well that my curves were cut in exactly the same number of segments by finding out how curves are simplified in a suite of segment in PrePoMax and FreeCAD.

I don’t exactly understand what you mean by:

the library and executable of NetGen mesher can be call by sub-process in Python script.
Do you mean using the executable NetGenMesher.exe used by PrePoMax to compute the mesh ?
If so, as far as I’ve found out by looking at the file _output_obj.txt in the temp folder, the calling of NetGenMesher.exe should be the following:

Path\PrePoMax v1.4.1\NetGen\NetGenMesher.exe BREP_MESH "Path\PrePoMax v1.4.1\Temp\geometry.brep" "Path\PrePoMax v1.4.1\Temp\geometry.vol" "Path\PrePoMax v1.4.1\Temp\meshParameters" "Path\PrePoMax v1.4.1\Temp\meshRefinement"

Unfortunately, I automatically generate my files from FreeCAD scripting and FreeCAD cannot export .vol and .brep files.

To give a bit more context: I am building a script that automatically iterates through designs and performs cyclic FEMs on them to optimize their geometry based on given criterions. I do my parametric design and FEM via FreeCAD. But FreeCAD hasn’t implemented yet 2D axisymmetric multi-material multi-contact, which was blocking for me. So I looked at how PrePoMax built its .inp file and wrote a script that adapts FreeCAD’s .inp file to one that is capable of running 2D axisymmetric FEAs.
But for some unknown reasons, I couldn’t get the same results. They were in the same order of magnitude but not similar enough to satisfy me. The only difference was in the mesh generation. FreeCAD can only refine meshes with Gmsh, which makes it kind of a “standard” when you want to perform these complicated simulations, as you need some refinement. But PrePoMax on its side uses netgen by default, and with refinements ! Therefore I want to find a way to mesh correctly with netgen only through python to link this part to the rest of my algorithm.
(This is way i’m not looking at scripting Gmsh, even though it’d have been way easier )
This brings us here.
Currently my simulation can run fairly well, but only in node to surface contact, instead of surface to surface when working completely from prepomax (which goes way faster).

Therefore, I’m still looking for PrePoMax’s standard settings/differences from the default of the netgen used by NGSolve. Are these results only different because of netgen’s version to you ?

Thank you again for your precious help !

[ANYS]

Could you create a Reference point at 0,0 to show it on the screen and check if both meshes are positioned the same way?
Be sure there are no elements in the y<0 plane.
Could you try the explode view to discard the different components are merged.?

Netgen can generate different meshes depending on the initial orientation of the model.

EDITED: Sorry I didn’t saw you posted the file.

[synt]

@wilzzer

Unfortunately no, the mesh size file follow exactly the same logic in structure : a series of lines defined by a start and an end point’s coordinates, followed by a value indicating the refinement.

it seems PrePoMax generates several internal points between its end when minimum and maximum is defined, or local refinement assigned to edges. Screenshot attached previously also shown some edge division is not the same, so i guessed it’s related to mesh size file. I’m not yet to check another edges since it is too large to count manually.

[Matej]

For me, the models do not work due to the Fixed support. For some reason, adding a boundary condition on the rotational degree of freedom does not work in CalculiX. So, replace the fixed boundary conditions with a displacement one.

Then you are right. One model can converge while the other does not. It seems the mesh density is to blame for that. I do not think there is a difference in how the mesh is prepared (using which mesher), but it is more of a contact-dependent problem. I would suggest using a mesh with a similar size of elements on the whole contact area, especially on the Sleeve part.

[wilzzer]

it seems PrePoMax generates several internal points between its end when minimum and maximum is defined, or local refinement assigned to edges. Screenshot attached previously also shown some edge division is not the same, so i guessed it’s related to mesh size file. I’m not yet to check another edges since it is too large to count manually.

Omg well noticed ! I definitely overlooked that. I dive again into the mesh size files to try and notice the differences. I was definitely convinced my mesh size files were the same

For me, the models do not work due to the Fixed support. For some reason, adding a boundary condition on the rotational degree of freedom does not work in CalculiX. So, replace the fixed boundary conditions with a displacement one.

As my first model, the one fully prepared with PrePoMax, does converge, will replacing fixed boundary with displacement change the output ? The displacement should be set as 0 in x, y, z and rY, right ? Or every displacement should be set to 0 ?

I do not think there is a difference in how the mesh is prepared (using which mesher)

But how come my two meshes are different ? Should it be theoretically to produce the exact same mesh using PrePoMax’s netgen and using Netgen from python directly ?

I would suggest using a mesh with a similar size of elements on the whole contact area, especially on the Sleeve part.

Most of the time, when a model doesn’t work I try to play with the mesh refinements along the edge to ensure a slave_element_site < master_element_size. In the end I generally can make models converge by tweaking this.

Be sure there are no elements in the y<0 plane.

Yes, i superposed the meshes and nothing is in the y<0 plane and both are in the same position (since they’re based on the same step files).

Could you try the explode view to discard the different components are merged.?

What is this for ? How does one do that ?

Netgen can generate different meshes depending on the initial orientation of the model.

Interesting, didn’t think of that either ! I’ll look into it.

Thank you all very much for your help !

[Matej]

@wilzzer

As my first model, the one fully prepared with PrePoMax, does converge, will replacing fixed boundary with displacement change the output ? The displacement should be set as 0 in x, y, z and rY, right ? Or every displacement should be set to 0 ?

The fixed support is the same as setting all the DOFs to 0. But setting UR3 to 0 somehow prevents the solver from running at all. And it will not impact your results if you do not constrain it.

@wilzzer

I do not think there is a difference in how the mesh is prepared (using which mesher)

But how come my two meshes are different ? Should it be theoretically to produce the exact same mesh using PrePoMax’s netgen and using Netgen from python directly ?

What I meant is that, in this case, creating exactly the same mesh will help, but in any other case, this is not the solution. Probably using PrePoMax mesh gives you a solution by chance and not while the PrePoMax mesher is different. So you need to find mesh settings that will work in most cases (make it more robust).

You should get identical mesh if you use exactly the same Netgen version and the same settings. Probably. It might be that the API you are using interprets some of the things differently.

[wilzzer]

What I meant is that, in this case, creating exactly the same mesh will help, but in any other case, this is not the solution. Probably using PrePoMax mesh gives you a solution by chance and not while the PrePoMax mesher is different. So you need to find mesh settings that will work in most cases (make it more robust).

Ooooh okay I get it now. I tried to replicate perfectly PrePoMax as its behavior seemed to be the best but you’re right that in my case I should just find a way to make things work all the time. As long as my results are proportional to reality depending on the design, I should be able anyways to extrapolate my most optimized design !

One more thing:
I noticed that settings node to surface contact or surface to surface contact would result in different results in the elongation of my mandrel. Currently, the pulling of the mandrel creates a slight material bump at the entry of the sleeve, therefore creating elongation earlier than wanted.
In the surface to surface contact simulation, the mandrel elongates a loooot, which I believe in reality would cause a break of the material. Whereas in the node to surface simulation, the elongation seems much more reasonable.
Therefore here are my questions:

• From the very technical point of view, how does node to surface behave ? How does surface to surface behave ? Is there any documentation of the maths/working principle behind these definitions ? I couldn’t find online any really technical explanations regarding contact definitions in FEAs.
• Which contact is the most “trustworthy” ?
• And finally, as node to surface would be instinctively the way to go here, are there any other part of an FEA that can be optimized to reduce runtime (apart from the meshsize i guess) ? Is process parallelization possible to your knowledge ? I know that this last question is not specific to PrePoMax so I understand if it’s a bit too far-fetch by I’m guessing your guys are quite used to FEAs processes

Thanks a lot !

[synt]

below i generate mesh using NetGen GUI, mesh size files taken from PrePoMax setting. Indeed, mesh results still have a bit different at curved contact regions. Maybe using previous versions of 5.x will give an insight.

[FEAnalyst]

@wilzzer

• From the very technical point of view, how does node to surface behave ? How does surface to surface behave ? Is there any documentation of the maths/working principle behind these definitions ? I couldn’t find online any really technical explanations regarding contact definitions in FEAs.
• Which contact is the most “trustworthy” ?

From the CalculiX User’s Manual:

Node to surface:

For each node on the dependent surface, a face on the independent surface is localized such that it contains the orthogonal projection of the node. If such is face is found a nonlinear spring element is generated consisting of the dependent node and all vertex nodes belonging to the independent face.

Surface to surface:

The spring element which was described in the previous section is now applied between an integration point of a slave face and a master face. The contact force at the integration point is subsequently transferred to the nodes of the slave face. This results in contact spring elements connecting a slave face with a master face. The integration points in the slave faces are not the ordinary Gauss points. Instead, the master and slave mesh are put on top of each other, the common areas, which are polygons (sides of quadratic elements are approximated by piecewise linear lines), are identified and triangulated. For each triangle a 7-node scheme is used. This can result to up to 100 or more integration points within one slave face. It usually leads to a very smooth pressure distribution. Furthermore, it is now irrelevant which side is defined as master and which as slave.

Generally speaking, surface to surface is recommended because of a better accuracy/distribution of the contact stresses, reduced risk of snagging and penetrations as well as reduced sensitivity with respect to the master/slave assignment when mesh densities are similar. Pretty much only disadvantage of the surface to surface contact is increased solution cost per iteration (since there are more nodes involved in constraints). On the other hand, this formulation may reduce the number of iterations needed to achieve convergence.

[Matej]

@wilzzer

Is process parallelization possible to your knowledge ? I know that this last question is not specific to PrePoMax so I understand if it’s a bit too far-fetch by I’m guessing your guys are quite used to FEAs processes

You can turn on parallel solving in Calculix in the PrePoMax settings: Tools → Settings → Calculi → Number of processors. Or you can prepare multiple models at the same time and run them in parallel.

[wilzzer]

Hi guys,

First of all thanks a lot for your advices and pointers. It really helps in my current work/reflection.

I decided to follow Matej’s advice:

What I meant is that, in this case, creating exactly the same mesh will help, but in any other case, this is not the solution. Probably using PrePoMax mesh gives you a solution by chance and not while the PrePoMax mesher is different. So you need to find mesh settings that will work in most cases (make it more robust).

And work with my current method, not focusing on making the “perfect mesh”. Unfortunately I noticed that there might be an issue with how I refine my mesh at the moment. Currently, If I use my meshRefinement file, my simulation runs like this and of course comes to a crash:

We can see that at the contact point the mesh gets completely distorted.
In a second try, I use, for the refinement of this piece, the meshRefinement file generated by PrePoMax. The meshRefinement file of the sleeve is the only thing that changes between bot simulations! Not only does the simulation not crash but there’s no more weird distortion in my result (picture taken roughly at the same displacement of my mandrel):

At this point, the meshRefinement files are REALLY simple as the refinement happens only on 3 edges and 1 curve (hence 6 refinement lines). Here they are:
Homemade meshRefinement file

0
6
2.967 2.31 0.0 2.967 6.91 0.0 0.1
2.1395 5.01 0.0 2.1395 2.523397 0.0 0.1
2.146199 2.498397 0.0 2.142515 2.506296 0.0 0.1
2.142515 2.506296 0.0 2.14026 2.514715 0.0 0.1
2.14026 2.514715 0.0 2.1395 2.523397 0.0 0.1
2.146199 2.498397 0.0 2.48591 1.91 0.0 0.1

PrePoMax meshRefinement file

0
6
2.967 2.31 0 2.967 6.91 0 0.1 
2.1395 5.01 0 2.1395 2.523397 0 0.1 
2.1395 2.523397 0 2.14026 2.514715 0 0.1 
2.14026 2.514715 0 2.142515 2.506296 0 0.1 
2.142515 2.506296 0 2.146199 2.498397 0 0.1 
2.48591 1.91 0 2.146199 2.498397 0 0.1 

The two first edges are described in the same way: firstVertex → lastVertex → refinementVal. Then comes the curve (discretized as 3 segments) and the last edge. For the curve, we can see that i get the exact same discretization but not displayed in the same order ! Same for the last edge: right values but not in the right order.

As this seems to impact so significantly the meshing and the results: how does PrePoMax determine which is the first vertex and which is the last one to place them correctly in its meshRefinement file ?

Thanks in advance !

All the best.

[Matej]

Not looking at the code, I will try to answer the question. So you have some connected line elements with 2 nodes each, and you want to sort them and their nodes in the same way as they are connected?

First, you find the endpoints. Those are the points that are only connected to one element. The easiest way to find them is to go through all line elements and increase the node counter for each single node as you encounter it.

Then you know where to start. Select one endpoint A, and find the element (element E) connected to it. Go through the element to the other side - the other node B. Then start from that node B and find all elements connected to it. One of the connected elements will be the element E, and the other element will be the next element on your path. Then repeat.

[Matej]

But in my tests of your model, it was much better if I used elements of the same size much deeper into the deformed material - mesh grading can be used to achieve this.

[wilzzer]

Hi Matej,

I’m not sure I completely understood your description. I believe maybe my first explanation was a bit too vague as I’m not looking for my mesh elements. I made a little illustration of the situation to make it clearer.
We know that a meshRefinement file has the following structure:

nP      ###Number of refined points
###Blablabla list of points and their refinement value
nS      ###Number of refined segments
x1 y1 z1 x2 y2 z2 refVal
x3 y3 z3 x4 y4 z4 refVal
...

Let P1, P2, P3, P4, … be the short name of point such that P1:(x1, y1, z1).
I’ve noticed that in the meshRefinement file, the order in which the two points describing a segment are written has an impact on the final mesh and its behavior during simulation (c.f. previous pictures). Therefore

0
1
x1 y1 z1 x2 y2 z2 refVal

Will not create the same mesh/meshing behavior as

0
1
x2 y2 z2 x1 y1 z1 refVal

I don’t have an issue finding the end points of the segments for my refinement but rather knowing how to arrange them. Which endpoint should be written first and which endpoint second ?
On to this issue in a specific case to illustrate it (c.f. last post). My refinements are applied during the meshing of the following piece. The two meshRefinement files previously cited are used in my script during the meshing of this piece and create the two separate results shown before.

In my case, i’m interested in refining Edge1, Edge2, Edge3 (which is a curve and not a line) and Edge4. Edges are numbered in this order based on the numbering done by FreeCAD. During my application, I want the user to type in Edge n°s based on how they’re named in FreeCAD. My script then runs through the edges in an ordered way (lowest edge nr to highest edge nr) and, among many other things, writes the meshRefinement file prior to meshing.
Note that to write the refinement of Edge3, which is a curve, I discretize it, resulting in 3 segments. And i get the exact same coordinates as PrePoMax does, which makes me think this is the way to go.With this in mind, I’ve detailed my previous homemade meshRefinement file:

And as a comparison, the details of the previously uploaded PrePoMax’s meshRefinement:

Therefore, we see that PrePoMax writes Edge3 and Edge4 in a reverse manner.
Currently, to chose the order in which an edge is written (i.e P2->P1 instead of P1->P2), I use the fact that FreeCAD stores edgeX.firstVertex and edgeX.lastVertex variables.

But so how does PrePoMax decide in which order to write its refined points ? Because it clearly has something that I don’t understand as my mesh ends up by being distorted completely.

The only thing that changes between the two simulations (ran from PrePoMax’s interface after importing the mesh) are these 4 little lines in the meshRefinement file.

But in my tests of your model, it was much better if I used elements of the same size much deeper into the deformed material - mesh grading can be used to achieve this.

I’ll get a look at this in the meantime! Thank you so much.

All the best!

[Matej]

I think that now I understand your question. I would have to look into the code to determine from where the edge directions are defined in PrePoMax. But again, this should not affect the resulting simulation. I dont think that your mesh is an any way different than PrePoMax mesh except in sizing.

If you change the mesh size in both approaches, does it affect the restults in the same way? PrePoMax mesh still works and yours mesh still doesn’t work? Try a couple of sizes.

Is the contact surface defined in the same way for both meshes?

[wilzzer]

Hi Matej!

But again, this should not affect the resulting simulation.

I thought so as well ! In my script I made a manual modification that simply changes the edge directions of these specific edges (i.e : cheating ahah). Then, I had the same directions as the one made by PrePoMax. The resulting mesh was still different BUT at contact point, it wouldn’t behave the same way as in my mesh_distortion.png file, uploaded earlier in this discussion. It’d behave normally, in a similar manner as mesh_no_distortion.png.

Try a couple of sizes.

I’ll try a few mesh sizes/mesh refinement sizes. I noticed that in my custom solution if reach a refinement 10 times smaller than in PrePoMax and define node to surface contacts instead, then the simulation runs fully but the mesh at the contact points still ends up with these weird extreme distortions.

Is the contact surface defined in the same way for both meshes?

Contact in PrePoMax is defined the same way. I’ll have a deeper look at the resulting .inp files to see if I can spot any differences.

Thanks again !

[wilzzer]

Hi Matej,

I dug into how you create your meshRefinement file. In Controller.cs, line 4010, the function CreateMeshRefinementFile is created. In there you use the list allLines and linesList to write meshRefinement.
These lists are built using the function GetVertexAndEdgeCoorFromGeometryIds from FeMesh.cs (line 6024). Itself retrieves the coordinates from GetCellsFromGeometryIds and GetCellsFromGeometryId.
But seemingly all these functions do is retrieve a cell/vertex corresponding to an ID (as explicitly stated by their names lol).
So it doesn’t seem like your edges are built in a specific direction or ordered in a specific manner.

You don’t seem to change their orientation at all or do anything particular ! Just taking the vertexes composing an edge (as one would do !)

This explains even less why so slight changes in the meshing results in a different behavior near contact points.

There’s definitely something I am missing here.

Back to digging!

[Matej]

At the time the mesh refinement is prepared no edge direction changes are applied. The edge directions are defined by the geometry, which is imported into PrePoMax via the geometry.vis file. You can find it in the temp folder. It contains the triangulated representation of the CAD model with the definition of the surfaces and edges, which is done using OpenCascade in the Netgen (PrePoMax uses Netgen to get the geometry). The points of the edges are taken from the OpenCascade directly, but I did not check if the edge orientation of the triangulation matches the orientation in the CAD kernel.