I am interested in running simulations on a perforated torus, however I am having trouble with the mesh generation. Currently, I create boundaries corresponding to a square and then another set of boundaries coming from circular obstacles via a poisson point process and use buildmesh to generate the resulting mesh. I then use periodic finite element spaces to identify the correct sides of the square.
However the meshing will not work when the obstacle intersect the boundaries of the square. Is there a clean way of doing this in FreeFem? Has anyone else already dealt with this problem?
I include some sample code with specifically chosen values for the obstacle positions however the goal is to have code robust against any given input, provided that none of the obstacles intersect each other and have centers in the square \[0,1\]\\times\[0,1\].
real r = 0.05; // Obstacle radius
real eps = 0.01; // epsilon for mesh
// Set obstacle configuration
int NumPart = 3;
real[int] xC(NumPart),yC(NumPart);
// (Does not work) wish to be robust against this
xC[0] = 0.5; yC[0]=0.0;
xC[1] = 0.0; yC[1]=0.05;
xC[2] = 0.5; yC[2]=0.5;
// Example points which will mesh
// xC[0] = 0.25; yC[0]=0.25;
// xC[1] = 0.75; yC[1]=0.75;
// xC[2] = 0.5; yC[2]=0.5;
// Borders definition for obstacles
border obst(t=0, 2*pi; i)
{
x = r*cos(t) + xC[i];
y = r*sin(t) + yC[i];
label = 4;
}
// Discretisation level for obstacles
int[int] nMesh(NumPart);
for(int i = 0; i < NumPart; i++){
nMesh[i] = -2*pi*r/eps ;
}
// Create borders for square
border a(t=0,1){x=t;y=0;label=0;};
border b(t=0,1){x=1;y=t;label=1;};
border c(t=0,1){x=t;y=1;label=2;};
border d(t=0,1){x=0;y=t;label=3;};
int nBorder = 1/eps;
// Build mesh
mesh Th = buildmesh(a(nBorder) + b(nBorder) + c(-nBorder) + d(-nBorder) + obst(nMesh),fixedborder=1 );
plot(Th,cmm = "Mesh");
fespace VhPeriodic(Th,P1,periodic=[[0,x],[2,x],[1,y],[3,y]]);
I think there is no “automatic” procedure to remove the outer part of each disc. What you can do is compute the intersection of each disc with the square by hand (in terms of the parameters), and then define the borders of the obstacles by only accounting for the part which is inside the square.
Thank you for responding. I had feared that might be the only way to do it, but it will be very messy to code so I hoped that there would be a more ‘‘elegant’’ solution.
Are there any ‘‘hidden’’ meshing tricks in FreeFem I could exploit? One idea I had was to create a larger square (an additional period in each direction) and by hand extend my holes periodically by a 1/2 period so that a hole towards the left will also appear on the right. Then intersect this mesh with a unit square and identify the relevant boundaries. Are there any efficient tools to do this: i.e. to intersect two meshes and then remesh the boundaries so that I can identify them for the periodic fespace?
My other idea was to see if it is possible embed the torus in R^3, and remove the holes there. However I am not sure how to do this.
If anyone has any examples doing something like this before, any help would be appreciated.
Up to my knowledge there is no tool to compute the intersection of two meshes. There is obviously a difficulty which is that in general the cells do not match correctly at the frontier.
Essentially there was no way around considering every possible case for the intersection of a circle with the boundary, then constructing the corresponding arcs and updating the borders for the exterior of the square. I have considered most cases, but as my application is for stochastic homogenisation, I have neglected some cases which occur with probability zero.
I have found that it works quite well for my application (low density perforations and with a uniform distance in between them). In cases where a set of points do not work, periodically shifting them by a small value seems to help.
