Purpose of the Meshing Applications
The Meshing application combines the features that do not belong to any specific physical module (be it fluid, thermal or structural applications), but can be used in conjunction with either of them. Basically Meshing application contains interfaces with meshers (Triangle, Tetgen, Nestor) and Mapping utilities, that permit data transfer from one mesh to another.
Meshing application is of particular importance for the problems that use moving meshes, specially the Particle Finite Element Method (PFEM), where mesh re-creation is obligatory at each time step.
Different versions of interfaces accounting for different problem requirements are listed in the next section.
The Meshing Application is using the following libraries:
1) Triangle (for triangulation): http://www.cs.cmu.edu/~quake/triangle.html
2) Tetgen 1.4.0 (for tetrahedrization): http://tetgen.berlios.de/
4) Msuit (alternative mesher for both 2D and 3D) - WEBADRESS MISSING!!!!
3) ANN library (for nearest neighbor/spatial search) (soon to be changed to kd-tree of Kratos): http://www.cs.umd.edu/~mount/ANN/
NOTE: it was noticed that the Tetgen 1.4.2 results in an error when used together with adaptive mesh refinement
Trigen PFEM refine
Variables Needed and their meaning
The following 4 variables need to be defined in order to perform the re-meshing using "Trigen PFEM refine":
By default, only the elements of the fluid are re-created. Therefore, you have to mark all the nodes of the elements you are passing to the mesher as "IS_FLUID". All the nodes that might be added during adaptive the refinement (see add_nodes argument of the function ReGenerateMesh) have IS_FLUID=1.0 by default.
IS_STRUCTURE switch is necessary to distinguish such parts of the domain as walls etc. Mesher treats the elements whose nodes have non-zero IS_STRUCTURE flag in a different way. The elements that have all nodes with IS_STRUCTURE=1.0 are automatically removed from the mesh.
IS_BOUNDARY is necessary to determine the boundaries, that include free surfaces, and boundaries between the fluid anb the structure (i.e. between the fluid and the wall).
IS_FREE_SURFACE is a flag that is set on the free-surfaces to 1.0 and 0.0 elsewhere.
if any of these variables are not added in your application, an error will be thrown.
Input of the function
For re-meshing using "Trigen PFEM refine" you need to execute one single function ReGenerateMesh. The function is defined as follows:
ReGenerateMesh (ModelPart& ThisModelPart , Element const& rReferenceElement, Condition const& rReferenceBoundaryCondition, NodeEraseProcess& node_erase, bool rem_nodes = true, bool add_nodes=true, double my_alpha = 1.4, double h_factor=0.5)
We shall have a closer look at each argument one has to pass to the function.
1. ModelPart& ThisModelPart:
this is the model part you would like to re-mesh. see (Link to model_part for its meaning)
2. Element const& rReferenceElement:
is the name (of type string) of the element that has to be re-created, e.g. if you are using standard incompressible fluid element of Kratos, you have to pass "Fluid2D" as an argument.
3. Condition const& rReferenceBoundaryCondition:
is the name of condition that is to be re-created. If you are not using any condition (such as Face Pressure, or any of your custom-defined condition), just pass "Condition2D". Otherwise the name of you condition (e.g. "FacePressure2D")
4. NodeEraseProcess& node_erase:
the function re-generate permits one to control the mesh quality. For this reason some nodes that cause mesh distortion (e.g. get too close to one another) can to be removed from the mesh. The "Trigen Mesh refine" identifies those nodes automatically, but to remove those nodes from the database, one needs to use the NodeEraseProcess, which is a standard process of Kratos. Therefore, you should create this process beforehand, and then pass it to the ReGenerateMesh function as an argument. See the example section to undesrtand exactly how to implement this.
5. bool rem_nodes:
is the switch used to allow derefinement, i.e. nodes removal. If you shall like to remove the nodes that come too close (see the h_factor explanation for an exact meaning of "too close"), you have to pass "True" here. If you do not want any node to be removed, pass "False". By default the switch is activated as "True".
6. bool add_nodes :
is a smiliar switch, now permitting refinement, i.e. adding nodes. Pass "True" or "False" depending on weather you would like to refine your mesh or not.
7. double my_alpha = 1.4:
is a alpha shape parameter, that defines a geomertrical criterion for recognition of the free surfaces. Values between 1.4 and 1.7 are usually used. For details see "Alpha Shape Method"
8. double h_factor = 0.5:
if you decide to control the mesh quality, you shouldnt allow nodes to approach each other "too much". This "too much" is defined as the fraction of the original length of the elemental edge e. If you set the h_factor to 0.5 all the nodes that find themselves closer to the given one than 0.5*e will be removed
Brief description of the algorithm
Trigen PFEM refine performs the re-meshing in 3 major steps.
Step 1 consists in identification and erasing (in case the bool erase_nodes is activated as "True") of the nodes that come too close to each other. A binary tree search is performed in order to find all the nodes that lie within h_factor*e radius around the given node. e is approximately equal to the initial inter-nodal distance of the undeformed mesh, that is stored before the beginning of the simulation. The nodes that lie within this radius will be removed.
Step 2 consists in the free surface identification and removal of the distorted elements. This is done using the alpha-shape techniques. The radius of the circles containing three nodes of an element are compared with an average element size. The final domain boundary and the preliminary mesh of the domain interior is created at this step.
Step 3 consists in the refinement of the mesh generated at Step 2 in case bool add_nodes is activated as "True". If wants to perform the "user-defined" adaptive refinement, one has to pass the desired element size to the mesher. This is done by prescribing the corresponding NODAL_H (that is roughly the average of the height lengths of the elements in the patch around the given node). Then the mesher will intent to create an element whose volume will be as close as possible to the equidistant triangle of the volume 0.5*(1.5*prescribed_h*1.5*prescribed_h), where prescibed_h=0.333(NODAL_H1+NODAL_H2+NODAL_H3), NODAL_Hi being the NODAL_H of the elements' nodes.
How to use it
We shall now illustrate which steps need to be taken in order to use the Meshing Interface. This is done in the python script.
First of all you need to import Meshing application and the interface:
applications_interface.Import_MeshingApplication = True
from KratosMeshingApplication import *
Then you must create a mesher (we give an example of using the TriGenPFEMModeler):
Mesher = TriGenPFEMModeler()
To execute the ReMesh function, you need to
Mesher.ReGenerateMesh("TestElement2D", "Condition2D", model_part, node_erase_process, remove_nodes, add_nodes, alpha_shape, h_factor)
TestElement2D is the name of your element, you wish to regenerate (can be e.g. Fluid2D)
Condition2D is the name of condition in case it exists
node_erase_process: is a process that permits you to erase nodes in case FLAG remove_nodes is set to TRUE in this case you need to create it prior to passing to this function by doing "node_erase_process = NodeEraseProcess(model_part)
add_nodes is a flag that should be set to TRUE in case you want to perform adaptive refinement
alpha_shape is a parameter (usually around 1.5) for free-surface detection
h_factor is the coefficient that serves as a criterion for identification of excessively close nodes (ratio between permissible internodal distance A and the element size: h_factor=A/h) "
Common Usage Errors
A very common error is not setting the IS_FLUID flag to 1.0 for the nodes of the model part that you pass for re-meshing. Then all the elements will be deleted. Do (in python):
for node in model_part_for_remeshing.Nodes: node.SetSolutionStepValue(IS_FLUID, 0, 1.0)
to make sure it is set.
Another error results, when the element name you pass to ReGenerateMesh does not coincide with the one inside the original *.elem file of your model.
This error can even beat ERROR n1 as to the frequency that it occurs. When you are using the TrigenPFEMRefine in conjunction with some application (say IncompressibleFluidApplication), do not forget to add the MeshingApplication to the application interface.
applications_interface.Import_MeshingApplication = True
If you are assigning manually the desired element size (by prescribing NODAL_H), make sure that the bool add_nodes switch of the ReGenerateMesh function is set to True. Otherwise alpha shape will delete the elements and now new ones will be created.
- Description of the Validation Example
Validation example is provided inside: adaptive_mesher2d.gid folder. This example performs refinement and de-refinement in a sequence of circular areas withing a rectangular domain. The analytical solution of the DISPLACEMENT_X is initially prescribed to all existing nodes at the beginning. Then the mesher adds and removes nodes, and the difference between the prescribed value and the value that results when the new node is inserted (and therefore interpolation is performed) is checked. Once this difference exceeds the prescribed tolerance (10E-8), an error is thrown if you are running an example in the benchamarking mode.
Next we show the snapshots of refinement:
- File to be looked at to get an example of use
Trigen Mesh Suite
Trigen Mesh suite is equivalent to the Trigen Pfem Refine with the unique difference that it does not allow adaptive refinement/derefinement.
Tetgen Mesh Suite
Tetgen Mesh suite is equivalent to the Tetgen Pfem Refine with the unique difference that it does not allow adaptive refinement/derefinement.
WORK IN PROGRESS
2D Constrained Delaunay and 3D Constrained Delaunay
Non-Overlappling Mesh Projection Methods
The operation of passing information between non matching meshes needs three different ingredients: - A searching algorithm to identify the neighbor nodes of a given element - A method to identify if a neighbor point is or not contained into a given element; - A projection function that pass the information from an element to a point contained in each element;
The searching algorithm is the key point of the efficiency of the method , in fact, as we will see in the following sections, it is the most time consuming part.
The searching algorithm
The problem of searching the neighbors of a given point is a very complicated and time consuming procedure. Nowadays many different data structures are use to try to optimize the problem. In Kratos many data structures are implemented to perform the neighbor search, among the different choices there are the kd-tree data structure and the bins data structure. In computer science, a kd-tree (short for k-dimensional tree) is a space- partitioning data structure for organizing points in a k-dimensional space. The kd-tree is based on a recursive subdivision of space into disjoint hyper- rectangular regions called cells. Each node of the tree is associated with such region, called a box, and is associated with a set of data points that lie within this box. The root node of the tree is associated with a bounding box that contains all the data points. Consider an arbitrary node in the tree. As long as the number of data points associated with this node is greater than a small quantity, called the bucket size, the box is split into two boxes by an axis-orthogonal hyperplane that intersects this box. A representation of how the kd-tree works can be seen in figg. 4.1 and 4.2. There are a number of different splitting rules, which determine how this hyperplane is selected . Another simple but effective data structure for storing and finding ob- jects in two and three dimensional spaces is Bins. It divides the domain into a regular nx × ny × nz sub-domains and holds an array of buckets storing its elements. Fig. 4.3 shows a two dimensional domain divided by bins and fig. 4.4 shows the structure of bins. This structure provides an extremely fast spatial searching when entities are more or less uniformly distributed over the domain. The good performance for well distributed entities and simplicity makes bins one of the most popular data structure in different finite element applications.
Strong mapping: interpolation through the nodes of different meshes
One possibility for the transmission of data between non matching meshes is the direct interpolation between the nodal values of the origin and the destination mesh. Given an origin model part and a destination model part, for each node of the destination mesh, we find out which element of the origin mesh it is contained in. Once the element is identified, a simple finite element interpolation with the shape functions of the element of the origin model part is performed. Finally the value of the variable in the new point is given. Once we found that only node i of fig. 4.5 is inside element ABC, we can calculate the value of variable q simply through FE interpolation
qi = NA (xi )qA + NB (xi )qB + NC (xi )qC ;
This procedure is very simple and give quite good results until the dimension of the two meshes we hare interpolating between, is not too different, otherwise other approaches have to be found to solve the problem.
How to use the projection utility
All the instructions to use the projection algorithm are included in How to use the projection between no matching meshes