Converting Nodal/Element Error into Metric (Scalar or Tensor2D/3D) for MMGProcess

[ssaxena93-tum]

I am writing a process within Convection-Diffusion Application to calculate error (L2 norm) in temperature gradient based on C1 and C0 continuous definitions of temperature gradient. The genesis of this question is to understand how to convert my error quantity, which may either be computed on nodes or gauss points, to a METRIC_SCALAR or METRIC_TENSOR_2D/3D quantity to interface with MMGProcess.py for the purposes of remeshing with MMG Library.

Below is the reply from Marc Nunez,

---

Frédéric Alauzet - Metric-Based Anisotropic Mesh Adaptation.pdf

Hi Saransh!

I can see you already have some knowledge on what a metric is and how MMG works, but in case you need more information about this, I have attached a book by F.Alauzet and P.Frey, where they introduce the concept of metric and some error-driven refinement strategies (I will refer to this later). Of course, it is not strictly necessary to go through the book to use the library but it can help to understand the working principle.

In short (and as you may know), what the MMG library needs to refine is an initial mesh and a metric field (defined at the nodes). With this, MMG will generate a new mesh in which each of the elements has a unitary volume if expressed in the metric provided. In Kratos we have already implemented (work by V.Mataix and R. Rossi) the interface between Kratos and MMG, so we only have to worry about computing the metrics (that it is done at Kratos level). Therefore, one can simply do (at python level)

mmg_process = KratosMultiphysics.MeshingApplication.MmgProcess2D(main_model_part, mmg_parameters)
mmg_process.Execute()

Where the main model part needs to already contain the metrics. How do we compute the metrics? One can define the metrics in any way, using the strategy that better suits your problem. Inside the MeshingApplication, you can find several strategies in the custom processes. I will describe a bit more two of these strategies later, but as a very first test, you can even define your metrics at python level!

For example, using the error that you mention, you could define the metrics as a function of the initial mesh size. We could say that if a node has an error above a certain tolerance, we will reduce its size by a factor of 10:

for node in main_model_part.Nodes:
     error=  self.GetError() # Here you would retrieve the error from your problem
     nodal_h = node.GetValue(KratosMultiphysics.NODAL_H) # Here you would need to use the  ind_nodal_h_process first
     if error > tolerance_error:
          nodal_h = nodal_h /10
     metric_scalar = 1 / (nodal_h**2) # The metric is related to the nodal_h by 1/h^2, more info on the book attached
     node.SetValue(KratosMultiphysics.MeshingApplication.METRIC_SCALAR, metric_scalar) # Setting th e METRIC_SCALAR variable

With this, you could directly call MMG to get the refined mesh. Of course, you can come up with any strategy to define your metric. Also, here I am using the METRIC_SCALAR, but both METRIC_TENSOR_2D and METRIC_TENSOR_3D could be used, that allow to define anisotrpy in your mesh.

What strategies are available in Kratos?

• metrics_error_process. I have not used this process, but it seems to be a similar case of what you want to do. It seems it expects the user to have defined an ELEMENT_ERROR in the elements, and it takes care of projecting it to the nodes. For further questions on using this process, I would suggest to open an issue in the Kratos repository tagging Vicente Mataix (@loumalouomega on github)
• metrics_hessian_process: This process is the main process that we normally use and that is described in the book attached. It uses the HESSIAN of a field (in your case it would be the temperature I assume), and uses it to estimate the interpolation error in the mesh, in order to refine more at the parts of the domain where this error is higher. Although you will be implementing your own metrics strategy, I suggest testing this one out, since it should also target the parts of the domain where the interpolation error is higher (maybe you can combine both?)

Here I have a cool image of the mesh that results from using the HESSIAN process on the velocity field in a fluid dynamics problem.

If these strategies are not suitable for your case, you can define your own metrics_***_process. As an extra example, here we were using a custom metric process that uses the divergence on the velocity field to refine.

To sum up: you have complete freedom to define the metrics as you like. The Kratos-MMG interface expects you to have defined a METRIC_SCALAR or a METRIC_TENSOR before using the MMG process. To do so, you can create a cpp process where you define the metric (for instance as a function on the nodal_h 1/h^2) or use one of the available ones.

Some of these strategies are included in the Examples repository: https://github.com/KratosMultiphysics/Examples/tree/master/mmg_remeshing_examples

There is also a python script to use the hessian process direclty from a .json:
https://github.com/KratosMultiphysics/Kratos/blob/master/applications/MeshingApplication/python_scripts/mmg_process.py

I hope this was helpful! For more specific questions you can contact me again or open an issue in the Kratos repository!

Regards,

[loumalouomega]

Dear @ssaxena93-tum, there are different ways to convert the error into a metric. The first one is compute the size of the mesh as it is traditionally done and compute the metric from this mesh size. This one will provide a isotropic mesh, and it is the procedure done in metrics_error_process. Another way is to compute the spatial distribution of the error so you can obtain a tensor with defines an anisotropic mesh.

[ssaxena93-tum]

Hi @loumalouomega Thanks for the info. @marcnunezc and @loumalouomega I was curious as to what the difference between the way Metric Scalar is defined between the metrics_error_process and the code snippet marc provides above, where he is using NODAL_H, which as I understand is the smallest edge length of the elements that share that node. Above Metric Scalar is defined as 1/(NODAL_H)^2 while in the cpp file, it is just defined as h_min which is either average of the minimum element size of the lowest element size. Shouldn't it also be some form of 1/(NODAL_H)^2 which in this case would simply be (1/h_min^2) ?

[loumalouomega]

First, since recently you can define a single scalar metric, see

Kratos/applications/MeshingApplication/custom_processes/mmg_process.cpp

Line 276 in b9c79c8

if (it_node_begin->Has(r_tensor_variable)) {

and

Kratos/applications/MeshingApplication/custom_processes/mmg_process.cpp

Line 302 in b9c79c8

double& r_metric = it_node->GetValue(METRIC_SCALAR);

. It was added here https://github.com/KratosMultiphysics/Kratos/pull/5778/files With the METRIC_SCALAR the metric corresponds directly with the mesh size, aka NODAL_H

[ssaxena93-tum]

Ah ok, that makes things clear. I had another couple of questions regarding the MMGProcess itself with regards to selective remeshing in coupled problems. So I have 2 approaches to test in the remeshing process,

1. A generic Conjugate Heat Transfer problem where I have conforming meshes and essentially all the nodes distributed between solid and fluid domain. Let's say globally there are 6105 nodes and 5444 are in fluid domain and 661 in solid domain while each domain has its own elements written in respective mdpa files. I would like to remesh only the fluid domain. However, if I pass only the Fluid Domain its shows this error wrt reading the metric.
   

Here's the fluid mdpa as a text. The node IDs till 6105 (the global numbering convention of IDs) while there are only 5444 nodes here since its only the fluid domain.
test_cylinder_cooling_Re100_Pr2_Fluid.txt

How do I go about remeshing in this case? Also, would it preserve the interface in this case between the solid and fluid domain or would I have to map it?

2. In cases where I remesh a generic fluid domain using the ISOSURFACE Level-Set approach using the Distance to Skin of the Solid Part and cut the fluid domain where the distance is <0, essentially the region where solid is supposed to be, is there a way to remesh in a way that I get a mapped skin part between the solid and fluid?

[loumalouomega]

First, I assume the solid and fluid share an interface. Am I right? You should fix the interface with "fix_contour_model_parts"

Kratos/applications/MeshingApplication/python_scripts/mmg_process.py

Line 119 in 7da8973

"fix_contour_model_parts" : [],

With this you will not solve the problem, as when re-meshing the mesh is rebuilt (so the nodes are destroyed). This can be solved, but I need a minimal example in order to do the corresponding changes.

[ssaxena93-tum]

Yes, the solid and fluid share an interface with nodes with same Ids in fluid and solid. I've attached the sample problem below. It is one of the most minimal problems I could find.
Benchmark Problem.zip

So if I understand correctly, this fix contour model parts separates the interface but it has to be fixed after remeshing since the meshes after remeshing would be nonconforming?

[loumalouomega]

I need to know the branch you are using

[ssaxena93-tum]

https://github.com/KratosMultiphysics/Kratos/tree/geo-cht

[loumalouomega]

Thanks, I was unable to run the stuff, I will merge with the master BTW

[ssaxena93-tum]

Sure, let me know if you need anything else from me.

[loumalouomega]

There are plenty of conflicts, when was the last time you updated the branch?

[ssaxena93-tum]

October, I believe.

[loumalouomega]

Yes, that explains a lot. I recommend to update regularly, specially with big changes

[loumalouomega]

@ssaxena93-tum I have been trying to adapt the MMGProcess, but it is a little bit more complicated that I expected. Before doing a "chapuza" (botched job) I prefer to agree a proceeding to synchronize shared meshes between model parts inside a model. @KratosMultiphysics/technical-committee which can be the better way to remesh a modelpart, which shares an interface with another model part, and syncronize the shared nodes in the interface when remeshing?

[jcotela]

@loumalouomega thanks for your suggestions, they are of great help in clarifying our problem. Just for me to understand. Lets consider for now our simplest possible problem, we have a single modelpart (solid+fluid) which covers the entire domain for the thermal problem and submodelparts for each domain (+ for the interface and the usual boundary condition contours). If we remesh the entire problem in this context, are the submodelparts preserved (in the sense that they still exist and have entities, even if the actual nodes/elements changed)? Would we be able to use the fix_contour_model_parts feature to keep the interface fixed in this case? (the interface modelpart is well defined, but not an external boundary). That would already be a helpful first step in our case.

[RiccardoRossi]

to give my two cents: modelparts and submodelparts are mantained after the remeshing (albeit with new entities in all of them according to the remehing) the processes apply bcs in the InitializeSolutionStep so if you call again that function the remeshed mp has all the conditions applied correctly

…

On Mon, Dec 23, 2019, 8:18 PM jcotela ***@***.***> wrote: @loumalouomega <https://github.com/loumalouomega> thanks for your suggestions, they are of great help in clarifying our problem. Just for me to understand. Lets consider for now our simplest possible problem, we have a single modelpart (solid+fluid) which covers the entire domain for the thermal problem and submodelparts for each domain (+ for the interface and the usual boundary condition contours). If we remesh the entire problem in this context, are the submodelparts preserved (in the sense that they still exist and have entities, even if the actual nodes/elements changed)? Would we be able to use the fix_contour_model_parts feature to keep the interface fixed in this case? (the interface modelpart is well defined, but not an external boundary). That would already be a helpful first step in our case. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#6133?email_source=notifications&email_token=AB5PWEOVDQIR57VDISFTVRDQ2EFKFA5CNFSM4J5GN4YKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHRYJ4I#issuecomment-568558833>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB5PWEOQFZYQ4KTHRVDLTCTQ2EFKFANCNFSM4J5GN4YA> .

[loumalouomega]

Yes, my concern is about shared nodes. This may still exist in the not remeshed model part but as the base hypothesis of this is to share nodes this will lead to memory errors. We should think a way to communicate this kind of changes between model parts belonging to a same model. I could do it, but it will take time, and as it is something that may be used elsewhere I prefer to agree first the way to proceed that to invent a way that will be not compatible in a near future. El lun., 23 dic. 2019 20:28, Riccardo Rossi <[email protected]> escribió:

…

to give my two cents: modelparts and submodelparts are mantained after the remeshing (albeit with new entities in all of them according to the remehing) the processes apply bcs in the InitializeSolutionStep so if you call again that function the remeshed mp has all the conditions applied correctly On Mon, Dec 23, 2019, 8:18 PM jcotela ***@***.***> wrote: > @loumalouomega <https://github.com/loumalouomega> thanks for your > suggestions, they are of great help in clarifying our problem. Just for me > to understand. Lets consider for now our simplest possible problem, we have > a single modelpart (solid+fluid) which covers the entire domain for the > thermal problem and submodelparts for each domain (+ for the interface and > the usual boundary condition contours). If we remesh the entire problem in > this context, are the submodelparts preserved (in the sense that they still > exist and have entities, even if the actual nodes/elements changed)? Would > we be able to use the fix_contour_model_parts feature to keep the interface > fixed in this case? (the interface modelpart is well defined, but not an > external boundary). That would already be a helpful first step in our case. > > — > You are receiving this because you were assigned. > Reply to this email directly, view it on GitHub > < #6133?email_source=notifications&email_token=AB5PWEOVDQIR57VDISFTVRDQ2EFKFA5CNFSM4J5GN4YKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHRYJ4I#issuecomment-568558833 >, > or unsubscribe > < https://github.com/notifications/unsubscribe-auth/AB5PWEOQFZYQ4KTHRVDLTCTQ2EFKFANCNFSM4J5GN4YA > > . > — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#6133?email_source=notifications&email_token=AEYQZAG75DRRQCPE3CF6A2DQ2EGPTA5CNFSM4J5GN4YKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHRY5CY#issuecomment-568561291>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEYQZAAOM4XLGF7SXT3FSSDQ2EGPTANCNFSM4J5GN4YA> .

[ssaxena93-tum]

Also a small pointer, not all submodel parts are preserved. Only those that have elements/conditions are preserved.
@loumalouomega is correct about the shared interface. Since only one side of the shared interface is remeshed, the 2 meshes from the 2 domains are no longer conforming, leading to memory errors. The nodes on the shared interface/boundary have to be mapped post remeshing.

[marcnunezc]

Hi @loumalouomega Thanks for the info. @marcnunezc and @loumalouomega I was curious as to what the difference between the way Metric Scalar is defined between the metrics_error_process and the code snippet marc provides above, where he is using NODAL_H, which as I understand is the smallest edge length of the elements that share that node. Above Metric Scalar is defined as 1/(NODAL_H)^2 while in the cpp file, it is just defined as h_min which is either average of the minimum element size of the lowest element size. Shouldn't it also be some form of 1/(NODAL_H)^2 which in this case would simply be (1/h_min^2) ?

I would like to add more info about this question to clarify concepts.

The eigenvalues of the metric are related to the edge length (h) of the element by lambda_i = h_i^-2.

To generate an isotropic mesh, we can either: build a diagonal METRIC_TENSOR defining the eigen values equally with value h_i^-2, or with the METRIC_SCALAR (which is already h). This means that a given isotropic metric tensor in a given point M(x), will generate elements of the size corresponding to the eigen values of M(x).

Regardless of how you define it, you can choose the size "h" that you want to prescribe in the mesh.

In the metric_error_process, an h_min is used, for instance:

h_min = 0.01
metric_scalar = h_min

The output of this will give an element of edge length 0.01.

In the code snippet I wrote initially, I chose to do this as a function of the current edge length of the element. So for this case we have:

current_h = node.GetValue(KratosMultiphysics.NODAL_H)
new_h = current_h / 10
metric_scalar = new_h

The output of this will give an element with edge length ten times less than the initial one.

Note that we are doing the same thing, defining the metric as a function of an edge length, but we choose this new "h" using a different strategy.

EDIT: updated as the Metric Scalar was wrongly defined as 1/h^2 as pointed out by Vicente #6133 (comment)

[loumalouomega]

@marcnunezc METRIC_SCALAR is directly h, you don't need to compute 1/h^2

[marcnunezc]

@marcnunezc METRIC_SCALAR is directly h, you don't need to compute 1/h^2

Ok sorry! I assumed it was the eigen value. I will edit my comment so it doesnt contain wrong info.

In any case, the point of having the "freedom" to define the metric as you please still stands!

[ssaxena93-tum]

@loumalouomega can "no_surf_mesh" option be used to preserve the interface? (Atleast in 3D)

[loumalouomega]

@loumalouomega can "no_surf_mesh" option be used to preserve the interface? (Atleast in 3D)

Is that an advanced option in mmg?

[ssaxena93-tum]

Kratos/applications/MeshingApplication/python_scripts/mmg_process.py

Line 137 in c6acb6c

"no_surf_mesh" : false,

There are other advanced options such no move mesh, no insert mesh etc

[loumalouomega]

Kratos/applications/MeshingApplication/python_scripts/mmg_process.py

Line 137 in c6acb6c

"no_surf_mesh" : false,

Yes, it is an advanced parameter de mmg

[loumalouomega]

See https://github.com/MmgTools/mmg/blob/79c6e96bd616adec0af55beecff762452bff3201/src/mmg3d/API_functions_3d.c#L90

[ssaxena93-tum]

I'll try it with that. But I still have issue w.r.t how to structure my workflow to remesh.
I have 2 mdpa files 1- fluid and 2- solid. The fluid mdpa also has nodes for solid domain for some reason and this was throwing the error -

6105 = total # of nodes, 5444 = nodes in fluid domain. How do I fix this?

I tried passing just the fluid part, however, that too is throwing errors -

[loumalouomega]

It is the same example as before or a different one?

[ssaxena93-tum]

Same as before. Although, the problem would arise is any coupled simulation setup. It's how GiD defines the two domains.

[loumalouomega]

In that case, same problem as before. We need to agreed a way to communicate the shared interface, so i can be remeshed

[ssaxena93-tum]

Shouldn't no_surf_mesh preserve the shared interface as it is a surface in the fluid/solid domain (For 3D atleast)?

[loumalouomega]

Even in that case (even if you block the interface), the old nodes will be removed and replaced with the equivalent (new entities, with new ids, but same coordinates and Nodal values), so the simulation will crash. In order to remesh we need to define a communication mechanism. El dom., 5 ene. 2020 23:30, ssaxena93-tum <[email protected]> escribió:

…

Shouldn't no_surf_mesh preserve the shared interface as it is a surface in the fluid/solid domain (For 3D atleast)? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#6133?email_source=notifications&email_token=AEYQZAFG3ACNKWXJREM534DQ4JNQZA5CNFSM4J5GN4YKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIEBGTY#issuecomment-570954575>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEYQZADPAKFKJ6PAXEWOC7DQ4JNQZANCNFSM4J5GN4YA> .