[Gmsh] Extremely small ratio of the smallest characteristic length to the size of the model leads to an error (self intersecting surface mesh)

[Juha Ritala]

Hi,

I have used the geometry module of Gmsh to model one quarter of a 
cylinder symmetric system representing the cantilever and tip of an 
atomic force microscope (AFM) on top of a sample (see the attached .geo 
file). My goal is to compute the electric field between the AFM probe 
and the back plate when a bias voltage is applied between those. The 
distance between the tip and the sample (h) can be as small as 1 nm, 
whereas the radius and height of the container (Lc) is 20 mm to justify 
setting the electrostatic potential to zero at the boundary.

If I set the distance between the AFM tip and the sample to 1 nm, the 
characteristic length of the mesh around the tip (tip_ls) must obviously 
be smaller than that. Using a characteristic length of exactly 1 nm 
usually results in a successful mesh generation, but if I set the 
characteristic length to be smaller than 1 nm, the 3D meshing process 
gives me a self intersecting surface mesh error. I suspect this is 
caused by the extremely small ratio of the characteristic length at the 
tip to the size of the model. If I set the characteristic length at the 
tip to 0.2 nm, this ratio is 10^-8.


I have no prior knowledge of mesh generation algorithms, but it seems to 
me that this is either a tolerance or a numerical precision problem. I 
have made the following attempts to solve the problem:

The most notable attempt I made was to change the tolerance of the 
Tetgen algorithm, which is the first part of the 3D meshing, and where 
the self intersecting surface mesh error arises. Changing the tolerance 
from the default value of 10^-8 to 10^-9 seems to have a clear positive 
impact on the success of meshing. Without this change, the meshing fails 
always if the characteristic length at the tip is set to 0.2 nm. With 
this change, I have managed to generate a good 3D mesh, although the 
meshing still randomly fails. I have not yet tried to set the tolerance 
even smaller, as changing the tolerance of the Tetgen requires 
recompilation of the source.

I have tried to tune the options that are on the first lines of my .geo 
file. The option Mesh.RandomFactor seems to have a clear impact on the 
meshing process. It seems that a smaller value is better, but the 2D 
meshing fails if the value is too small. According to the documentation, 
too small value of Mesh.RandomFactor leads to problems with numerical 
precision.


The questions I have are:
Has someone else had this kind of a problem? Am I on the right with my 
solution attempts? Is there something else that can be done? Where is 
the value of Mesh.RandomFactor actually used?


Best regards,
Juha

-- 
Juha Ritala
Department of Applied Physics
Aalto University, Finland
juha.ritala at aalto.fi

-------------- next part --------------
// Options
Geometry.Tolerance = 1e-11;
//Mesh.CharacteristicLengthMin = 1e-9;
Mesh.CharacteristicLengthMax = 2e-3;
//Mesh.CharacteristicLengthFactor = 1;
Mesh.RandomFactor = 1e-11;
Mesh.ToleranceEdgeLength = 1e-10;
Mesh.LcIntegrationPrecision = 1e-10;

// Parameters
H = 15e-6;
R = 20e-9;
Rd = 35e-6;
td = 0.5e-6;
ts = 1e-3;
h = 1e-9;
Lc = 1e6*R;
alpha = 15*Pi/180;
beta = Pi/2-alpha;

// Characteristic lengths
tip_ls = 1e-9;
cone_ls = 2e-7;
cantilever_ls = 1e-6;
far_ls = 2e-3;
sample_bottom_ls = 1e-3;


// Geometry of the probe
p = newp;
Point(p) = {0, 0, h, tip_ls};
Point(p+1) = {0, R*Sin(beta), h+R-R*Cos(beta), tip_ls};
Point(p+2) = {0, R*Sin(beta)+H*Tan(alpha), h+R-R*Cos(beta)+H, cone_ls};
Point(p+3) = {0, Rd, h+R-R*Cos(beta)+H, cantilever_ls};
Point(p+4) = {0, Rd, h+R-R*Cos(beta)+H+td, cantilever_ls};
Point(p+5) = {0, 0, h+R-R*Cos(beta)+H+td, cantilever_ls};
Point(p+6) = {0, 0, h+R, tip_ls};
p_probe_top = p+5;
p_probe_bottom = p;

l = newl;
Circle(l) = {p,p+6,p+1};
Line(l+1) = {p+1,p+2};
Line(l+2) = {p+2,p+3};
Line(l+3) = {p+3,p+4};
Line(l+4) = {p+4,p+5};
outer_line_probe[] = {l:l+4};

surf_probe[] = {};
For j In {0:4}
  tmp[] = Extrude {{0,0,1},{0,0,0},-Pi/2} {Line{l+j};};
  surf_probe[] += tmp[1];
EndFor


// The air cylinder
p = newp;
Point(p) = {0, 0, 0, tip_ls};
Point(p+1) = {0, Lc, 0, far_ls};
Point(p+2) = {0, Lc, Lc, far_ls};
Point(p+3) = {0, 0, Lc, far_ls};
p_sample_center = p;
p_sample_edge = p+1;

l = newl;
Line(l) = {p_probe_top,p+3};
Line(l+1) = {p+3,p+2};
Line(l+2) = {p+2,p+1};
Line(l+3) = {p+1,p};
Line(l+4) = {p,p_probe_bottom};
line_sample_top = l+3;

ll = newll;
Line Loop(ll) = {l:l+4, outer_line_probe[]};
s = news;
Plane Surface(s) = {ll};

tmp[] = Extrude {{0,0,1},{0,0,0},-Pi/2} {Surface{s};};
vol_air = tmp[1];
pb_surfs[] = {s,tmp[0]};
surf_sample_top = tmp[4];


// The sample
p = newp;
Point(p) = {0, 0, -ts, sample_bottom_ls};
Point(p+1) = {0, Lc, -ts, far_ls};
p_bottom_center = p;

l = newl;
Line(l) = {p_sample_center,p};
Line(l+1) = {p,p+1};
Line(l+2) = {p+1,p_sample_edge};

ll = newll;
Line Loop(ll) = {l:l+2, line_sample_top};
s = news;
Plane Surface(s) = {ll};

tmp[] = Extrude {{0,0,1},{0,0,0},-Pi/2} {Surface{s};};
vol_sample = tmp[1];
pb_surfs[] += {s,tmp[0]};


// Back-plate and infinite boundary
surf_inf[] = CombinedBoundary{ Volume{vol_air,vol_sample}; };
surf_inf[] -= {surf_probe[],pb_surfs[]};


// Mesh size fields
exp1=1.1;
a = (far_ls-tip_ls)/(Lc/2)^exp1;
b = tip_ls;
Printf("a: %g", a);
Printf("b: %g", b);
Printf("h: %g", h);

Field[1] = Box;
Field[1].VIn = a;
Field[1].VOut = 3*a;
Field[1].XMin = -1e-9;
Field[1].XMax = Lc+1e-9;
Field[1].YMin = -1e-9;
Field[1].YMax = Lc+1e-9;
Field[1].ZMin = 0;
Field[1].ZMax = Lc+1e-9;

Field[2] = MathEval;
Field[2].F = Sprintf("F1*Sqrt(x^2+y^2+(z-%g)^2)^%g + %g", h, exp1, b);
Background Field = 2;


// Define the physical objects
PROBE = 1000;
SAMPLE = 1001;
AIR = 1002;
GROUND = 1003;
Physical Surface(PROBE) = {surf_probe[]};
Physical Volume(SAMPLE) = {vol_sample};
Physical Volume(AIR) = {vol_air};
Physical Surface(GROUND) = {surf_inf[]};