Mesh3 from MeshS with obstacles

General Discussion
1

Dear all,

I have a surface mesh (obstacle.mesh).

Next, I create the surface mesh of the bounding box such that the obstacle.mesh is fully including in it.

Then when I want to merge the two surface mesh as

meshS MESHTOT = SURFbox + obstacle;

I get

current line = 25
Assertion fail : (hmin > Norme2(Pn - Px) / 1e9)
line :2112, in file msh3.cpp
Assertion fail : (hmin > Norme2(Pn - Px) / 1e9)
line :2112, in file msh3.cpp
err code 6 , mpirank 0

Does someone know what is the problem with my mesh ?

I put the mesh and the code to reproduce the error in attachment.

Thank you in advance,

Best regards,

Loïc,

MWE_Loic.edp (902 Bytes)
obstacle.mesh (2.3 MB)

2

I think that your mesh has some topological problem, since FF says “nbnomanifold=158 Warning no manifold obj nb:80 adj 158 of dim =2”

3

Thanks @fb77. Indeed, I have been able to repair my surface mesh by removing small artifacts in it.

However, now I have a small problem. I want to create a volume mesh of the space between the box and the obstacle.

To do this kind of mesh, previously I was using TetGen by defining the holes and the domain.

However in this case the mesh obstacle are disconnected (in several parts) and I do not know a priori the exact localization of the obstacles.

Is there any way to do it (with Tetgen or others) without defining explicitly the holes ?

I have put the repaired mesh in attachment and the code to load the surface mesh.

Thank you in advance for your help.

Best regards,

Loïc,

image
image694×636 50.2 KB

MWE_Loic.edp (912 Bytes)
obstacle_repaired.mesh (2.2 MB)

4

Normally you should do something like

real[int] domain = [9.,9.,1.,0,1.e-3];
real[int] myholes = [10.,10.,2.,9.5,10.5,2.5];
mesh3 Th=tetg(MESHTOT,switch="paAYY",nbofregions=1,regionlist=domain,holelist=myholes);

the “domain” gives the 3 coordinates of a point in your desired domain, the label you want for this region and the max volume of the tetrahedra.
The “holelist” gives a list of 3 coordinates of points in the holes (hence 3 * number of holes numbers). Here I have put 2 holes but I took the coordinates randomly so they are wrong.
You should know at least one point in each hole in order to do this correctly and exclude the holes from your mesh.

5

Yes, exactly, @fb77. I usually do this, for example, when I create spheres or cylinder obstacles, i.e., when I know the exact positions of the obstacles.

But now, since the obstacle geometry is given, I was wondering if there’s a way to automatically fill the space in the box without manually specifying one point per hole. That doesn’t seem feasible in this case, so I was hoping you might know of a simple way to (semi-automatically) extract one point per hole, since I have many such domains to process. I don’t know actually how to get easily one point per holes in case of several obstacles whose I don’t know the locations.

Thank you for help,

Loïc,

6

I think that tetgen is not able to do that automatically. I would try, knowing the meshS, to make a loop on its triangles and generates 2 points close to each triangle, one on one side and one on the other side (or maybe with orientation one can select always the same side). And then try to separate the connected components… Indeed it means making the job that tetgen does not, it is a pain. Unless you can do it with another software.

7

I have been able to overcome the difficulty by important each part of the mesh separately as

load “medit”
include “cube.idp”
load “tetgen”
include “MeshSurface.idp”

/************************* SURFACE MESH OF THE COIL **********************/
meshS ThS0 = readmeshS(“cut1_0.mesh”);
meshS ThS1 = readmeshS(“cut1_1.mesh”);
meshS ThS2 = readmeshS(“cut1_2.mesh”);
meshS ThS4 = readmeshS(“cut1_4.mesh”);
meshS ThS5 = readmeshS(“cut1_5.mesh”);

/************** SURFACE MESH OF THE BOUNDING BOX **************************/

real eta = 0.2;
real[int,int] BB=[[11-eta,13+eta],[9-eta,11+eta],[1-eta,3+eta]]; // bounding bax
int[int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right
int[int] N = [20, 20, 20];
meshS SURFbox = SurfaceHex(N, BB, L, 1); //discretization, size, label, orientation

/************************* COMBINE SURFACE MESH ***********************************/
meshS MESHTOT = SURFbox + ThS0 + ThS1 + ThS2 + ThS4 + ThS5 ;

/************************* CREATE VOLUME MESH ***********************************/

real x0 = ThS0[0][0].x;
real y0 = ThS0[0][0].y;
real z0 = ThS0[0][0].z;

real x1 = ThS1[0][0].x;
real y1 = ThS1[0][0].y;
real z1 = ThS1[0][0].z;

real x2 = ThS2[0][0].x;
real y2 = ThS2[0][0].y;
real z2 = ThS2[0][0].z;

real x4 = ThS4[0][0].x;
real y4 = ThS4[0][0].y;
real z4 = ThS4[0][0].z;

real x5 = ThS5[0][0].x;
real y5 = ThS5[0][0].y;
real z5 = ThS5[0][0].z;

real xb = 11-eta;
real yb = 9-eta;
real zb = 1-eta;

real[int] domain = [xb,yb,zb, 53, 1.e-4];
real[int] myholes = [x0,y0,z0, x1,y1, z1, x2,y2, z2, x4,y4, z4, x5,y5, z5];
mesh3 Th = tetg(MESHTOT, switch=“paAAQYY”, holelist=myholes, nbofregions=1, regionlist=domain);
medit(“mesh”, Th);
savemesh(Th,“RVE1.mesh”);

Everything works well.

But now, my last problem:
I
want to apply periodic boundary conditions on the border of the box as

fespace Vh(Th, [P2,P2,P2,P1], periodic=[[1, x, z], [3, x, z], [2, y, z], [4, y, z], [5, x, y], [6, x, y]]);

But I get this error:

Assertion fail : (hmn>1.0e-20)
line :215, in file lgmesh3.cpp

Normally, when we construct the box using SurfaceHex(N, BB, L, 1), the opposite faces of the box are generated with matching meshes, which allows to support periodic boundary conditions. Therefore, I’m not sure what is causing the issue here.

Could you have an idea how to fix this problem ?

I put all the materials to reproduce the problem in attachment.

Thank you in advance,

Best regards,

Loïc,
cut1_0.mesh (40.8 KB)
cut1_1.mesh (2.7 MB)
cut1_2.mesh (80.0 KB)
cut1_4.mesh (68.2 KB)
cut1_5.mesh (63.1 KB)
MWE_PER.edp (1.9 KB)

8

Great that you succeeded in finding points in the holes!
Your labels are not correctly defined, it should be
int[int,int] L=[[4,2],[1,3],[5,6]]; // the label of the 6 face left,right, front, back, down, up

9

Thanks! It works.

Loïc,