[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
juha.ritala at aalto.fi
Thu Feb 20 15:18:02 CET 2014
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[]};