[Gmsh] meshing trouble

Wed Feb 15 06:00:14 CET 2017

Hi, 

I have fought with this issue for a long time, but still can't resolve it, any help would be appreciated. I have solid objects in the air (with spherical out surface) that needed to be excluded. Unfortunately the solid objects (a frame and two magnets) touch each other. So I made a surface loop to include all the solids to define the air body. Gmsh can't mesh the model, complaining about intersecting surfaces. If you delete the frame (see code comment on how to delete), everything else being the same, Gmsh can mesh it without any problems. 

Very strange!

Regards,
Yingjie

=====================Code starts from here (also in attachment)========================

// define geometry-specific parameters
mm = 1.e-3;
gap = 1*mm; //Magnet-Magnet gap [m], gap > 0

DefineConstant[
  cub = {50*mm, Name "Parameters/1Magnet bottom size [m]"}
  hite = {10*mm, Name "Parameters/2Magnet hieght [m]"}
  fdep = {60*mm, Name "Parameters/3Frame depth [m]"}
  inf = {100*mm, Name "Parameters/4Air box distance [m]"}
  lc1 = {TotalMemory <= 2048 ? 30*mm : 10*mm, Name "Parameters/5Mesh size on magnets [m]"}
  lc2 = {TotalMemory <= 2048 ? 90*mm : 20*mm, Name "Parameters/6Mesh size at infinity [m]"}
];

// change global Gmsh options
Mesh.Optimize = 1; // optimize quality of tetrahedra
Mesh.VolumeEdges = 0; // hide volume edges
Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)

//parameters used by this macro
// cx, cy, cz: center of sphere
// rad: radius of sphere
// slc: mesh size

Macro MakeCircleZ //facing z direction
  cp = newp; Point(cp) = {cx, cy, cz, slc}; //center
  dd[] = {rad, 0, -rad, 0};
  For i In {0:3}
    spoint[i] = newp;
Point(spoint[i]) = {cx+dd[i], cy+dd[3-i], cz, slc};
If (i>0)
  quarc[i-1] = newreg; //quarter arcs
  Circle(quarc[i-1]) = {spoint[i-1], cp, spoint[i]};
EndIf
  EndFor
  //closing arc
  quarc[3] = newreg; Circle(quarc[3]) = {spoint[3], cp, spoint[0]};
  zloop = newreg; Line Loop(zloop) = {quarc[]};
  zcircle = newreg; Plane Surface(zcircle) = {zloop};
Return

cx = 0; cy=0; cz = -gap-hite;
rad = cub; slc = lc1;
Call MakeCircleZ;
mag1[] = Extrude {0, 0, hite} { Surface{zcircle}; };
jfl1 = zloop; //joint face loop 1

cx = 0; cy=0; cz = gap+hite;
rad = cub; slc = lc1;
Call MakeCircleZ;
mag2[] = Extrude {0, 0, -hite} { Surface{zcircle}; };
jfl2 = zloop; //joint face 2

//Delete {Volume{mag1[1]};}
//Delete {Volume{mag2[1]};}
Physical Volume("Magnets") = {mag1[1], mag2[1]};

//create steel frame around the magnet
p1 = newp; Point(p1) = {-2*cub, -fdep, -hite-gap, lc1};
p2 = newp; Point(p2) = { 2*cub, -fdep, -hite-gap, lc1};
p3 = newp; Point(p3) = { 2*cub, -fdep,  hite+gap, lc1};
p4 = newp; Point(p4) = {-2*cub, -fdep,  hite+gap, lc1};
l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};

q1 = newp; Point(q1) = {-2*cub, fdep, -hite-gap, lc1};
q2 = newp; Point(q2) = { 2*cub, fdep, -hite-gap, lc1};
q3 = newp; Point(q3) = { 2*cub, fdep,  hite+gap, lc1};
q4 = newp; Point(q4) = {-2*cub, fdep,  hite+gap, lc1};
e1 = newl; Line(e1) = {q1,q2}; e2 = newl; Line(e2) = {q2,q3};
e3 = newl; Line(e3) = {q3,q4}; e4 = newl; Line(e4) = {q4,q1};

d1 = newl; Line(d1) = {p1,q1}; d2 = newl; Line(d2) = {p2,q2};
d3 = newl; Line(d3) = {p3,q3}; d4 = newl; Line(d4) = {p4,q4};

ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};

lex1 = newll; Line Loop(lex1) = {l1, d2, -e1, -d1};
fex1 = news; Plane Surface(fex1) = {lex1}; 
hof1 = news; Plane Surface(hof1) = {lex1, jfl1}; //holed face 1

lex2 = newll; Line Loop(lex2) = {l2, d3, -e2, -d2};
fex2 = news; Plane Surface(fex2) = {lex2}; 

lex3 = newll; Line Loop(lex3) = {l3, d4, -e3, -d3};
fex3 = news; Plane Surface(fex3) = {lex3}; 
hof2 = news; Plane Surface(hof2) = {lex3, jfl2}; //holed face 2

lex4 = newll; Line Loop(lex4) = {l4, d1, -e4, -d4};
fex4 = news; Plane Surface(fex4) = {lex4}; 

ll1d = newll; Line Loop(ll1d) = {e1,e2,e3,e4};

hite2 = hite + gap + 1.5*cub;
p1 = newp; Point(p1) = {-3.5*cub, -fdep, -hite2, lc1};
p2 = newp; Point(p2) = { 3.5*cub, -fdep, -hite2, lc1};
p3 = newp; Point(p3) = { 3.5*cub, -fdep,  hite2, lc1};
p4 = newp; Point(p4) = {-3.5*cub, -fdep,  hite2, lc1};
l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};

q1 = newp; Point(q1) = {-3.5*cub, fdep, -hite2, lc1};
q2 = newp; Point(q2) = { 3.5*cub, fdep, -hite2, lc1};
q3 = newp; Point(q3) = { 3.5*cub, fdep,  hite2, lc1};
q4 = newp; Point(q4) = {-3.5*cub, fdep,  hite2, lc1};
e1 = newl; Line(e1) = {q1,q2}; e2 = newl; Line(e2) = {q2,q3};
e3 = newl; Line(e3) = {q3,q4}; e4 = newl; Line(e4) = {q4,q1};

d1 = newl; Line(d1) = {p1,q1}; d2 = newl; Line(d2) = {p2,q2};
d3 = newl; Line(d3) = {p3,q3}; d4 = newl; Line(d4) = {p4,q4};

ll2 = newll; Line Loop(ll2) = {l1,l2,l3,l4};

let1 = newll; Line Loop(let1) = {l1, d2, -e1, -d1};
fet1 = news; Plane Surface(fet1) = {let1}; 

let2 = newll; Line Loop(let2) = {l2, d3, -e2, -d2};
fet2 = news; Plane Surface(fet2) = {let2}; 

let3 = newll; Line Loop(let3) = {l3, d4, -e3, -d3};
fet3 = news; Plane Surface(fet3) = {let3}; 

let4 = newll; Line Loop(let4) = {l4, d1, -e4, -d4};
fet4 = news; Plane Surface(fet4) = {let4}; 

ll2d = newll; Line Loop(ll2d) = {e1,e2,e3,e4};

s1 = news; Plane Surface(s1) = {ll2, ll1};
s1d = news; Plane Surface(s1d) = {ll2d, ll1d};

//the frame, comment out the next five lines to delete it
slframe = newsl; 
Surface Loop(slframe)={ s1, fex1, fex2, fex3, fex4,
                        s1d, fet1, fet2, fet3, fet4};
frame = newv; Volume(frame) = {slframe};
Physical Volume("Frame") = {frame};

//now create the surface enclosing both mags and the frame
solids = newsl; 
Surface Loop(solids) = { s1d, fet1, fet2, fet3, fet4, s1,
    hof1, mag1[2], mag1[3], mag1[4], mag1[5], mag1[0], fex2,
    hof2, mag2[2], mag2[3], mag2[4], mag2[5], mag2[0], fex4};

//parameters used by this macro
// cx, cy, cz: center of sphere
// rad: radius of sphere
// slc: mesh size of sphere

Macro MakeSphereSurface
  cp = newp; Point(cp) = {cx, cy, cz, slc}; //center
  np = newp; Point(np) = {cx, cy, cz+rad, slc}; //north
  sp = newp; Point(sp) = {cx, cy, cz-rad, slc}; //south
  dd[] = {rad, 0, -rad, 0};
  For i In {0:3}
    ap = newp;  pe[i] = ap; //points on the equator
Point(ap) = {cx+dd[i], cy+dd[3-i], cz, slc};
na[i] = newreg; Circle(na[i]) = {ap, cp, np}; //north arc
sa[i] = newreg; Circle(sa[i]) = {ap, cp, sp}; //south arc
If (i>0)
  eq[i-1] = newreg; //equator
  Circle(eq[i-1]) = {pe[i-1], cp, ap};
EndIf
  EndFor
  //closing arc of the equator
  eq[3] = newreg; Circle(eq[3]) = {ap, cp, pe[0]};
  //now the non-plane surfaces: sloop[]
  For i In {0:3}
    lloop = newreg; Line Loop(lloop)={eq[i], -na[i], na[(i+1)%4]};
sloop[2*i] = newreg; Ruled Surface(sloop[2*i]) = {lloop};
    lloop = newreg; Line Loop(lloop)={eq[i], -sa[i], sa[(i+1)%4]};
    sloop[2*i+1] = newreg; Ruled Surface(sloop[2*i+1]) = {lloop};
  EndFor
  thesurf = newreg;
  Surface Loop(thesurf) = {sloop[]};
Return

// create air box around magnets
BoundingBox; // recompute model bounding box
cx = (General.MinX + General.MaxX) / 2;
cy = (General.MinY + General.MaxY) / 2;
cz = (General.MinZ + General.MaxZ) / 2;
lx = 2*inf + General.MaxX - General.MinX;
ly = 2*inf + General.MaxY - General.MinZ;
lz = 2*inf + General.MaxZ - General.MinZ;
rad = (lx + ly + lz)/3;
slc = lc2;

Call MakeSphereSurface ;

vair = newv; Volume(vair) = {thesurf, solids};
Physical Volume("Air") = {vair}; // air
Physical Surface("Infty") = sloop[];

// adjust value here for correct merge result
//Geometry.Tolerance = 1e-6; 
//Coherence Mesh;
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://onelab.info/pipermail/gmsh/attachments/20170215/5cf7306b/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bug2cylmag.geo
Type: application/octet-stream
Size: 7189 bytes
Desc: not available
URL: <http://onelab.info/pipermail/gmsh/attachments/20170215/5cf7306b/attachment-0001.geo>