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;
bug2cylmag.geo
Description: Binary data
_______________________________________________ gmsh mailing list [email protected] http://onelab.info/mailman/listinfo/gmsh
