Dear Gmsh team, Recently I encountered an issue when meshing some geometries with ellipse arc, and after some time working on this issue, I think currently gmsh cannot deal with ellipses that start and end points are symmetric with respect to major axis (or minor axis). All line numbers and reference to code are based on 3.0.6 source code. example.geo is a small example that can reproduce this issue on 3.0.6 windows version, and GmshTest.py is a python script based on gmsh ellipse handling code that can print variable values during processing (I do not have gmsh build environment on my machine..). Just a note that I was using the problem I encountered as test case (instead of example.geo) in GmshTest.py.
Basically what I found is that in Geo/Geo.cpp, near line 364, if for (x1,x3) and (y1,y3), one pair has the same value and the other has the same absolute value but opposite sign, sys matrix will have two identical rows and sys2x2() on line 371 will fail. If I understand the code correctly, the code here is primarily trying to calculate the semi-major axis length f1 and semi-minor axis length f2. In case that two points are symmetric wrt axis after rotation, I think we have to use the major axis point coordinate to get f1, and only solve f2 using the coordinate of an end point. I can work around this issue in my work, though it will be nice if Gmsh is able to handle this case.. Thanks, Shengjie
from math import *
def norm3(vec):
return sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2])
def norme(vec):
d=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2])
if d==0.0:
return vec
else:
return [vec[0]/d,vec[1]/d,vec[2]/d]
def prodve(a,b):
return [a[1] * b[2] - a[2] * b[1],-a[0] * b[2] + a[2] * b[0],a[0] *
b[1] - a[1] * b[0]]
def direction(a,b):
return [b[0]-a[0],b[1]-a[1],b[2]-a[2]]
def projette(vec,mat):
X = vec[0] * mat[0][0] + vec[1] * mat[0][1] + vec[2] * mat[0][2]
Y = vec[0] * mat[1][0] + vec[1] * mat[1][1] + vec[2] * mat[1][2]
Z = vec[0] * mat[2][0] + vec[1] * mat[2][1] + vec[2] * mat[2][2]
return [X,Y,Z]
def gmsh_SQU(x):
return (x*x)
def sys2x2(mat,b):
norm = gmsh_SQU(mat[0][0]) + gmsh_SQU(mat[1][1]) + gmsh_SQU(mat[0][1])
+ gmsh_SQU(mat[1][0])
det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]
res=[0.0,0.0]
if norm == 0.0 or (fabs(det) / norm) < 1.e-12:
print('sys2x2:
norm='+str(norm)+',det='+str(det)+',det/norm='+str(fabs(det)/norm))
return res
ud = 1.0 / det
res[0] = b[0] * mat[1][1] - mat[0][1] * b[1]
res[1] = mat[0][0] * b[1] - mat[1][0] * b[0]
res[0]*=ud
res[1]*=ud
return res
def myasin(a):
if a <= -1.0:
return -pi / 2.0
elif a >= 1.0:
return pi / 2.0
else:
return asin(a)
def myatan2(y,x):
if y == 0.0 and x == 0.0:
return 0.0;
return atan2(y, x);
def angle_02pi(A):
DP = 2 * pi
res=A
while res > DP or res < 0.0:
if res > 0:
res -= DP
else:
res += DP
return res
def list2str(list):
return '['+','.join(str(e) for e in list)+']'
def mat2str(mat):
str='[\n'
for list in mat:
str+=list2str(list)
str+='\n'
str+=']\n'
return str
def trace(firstPoint,center,lastPoint,majorAxisPoint):
fmt_str='{0:3}: {1:5}= {2}'
v=[firstPoint,center,lastPoint,majorAxisPoint]
dir12=direction(v[1],v[0])
dir32=direction(v[1],v[2])
dir42=direction(v[1],v[3])
print(fmt_str.format(269,'dir12',list2str(dir12)))
print(fmt_str.format(270,'dir32',list2str(dir32)))
print(fmt_str.format(272,'dir42',list2str(dir42)))
v0=dir12
v2=dir32
v3=dir42
dir12=norme(dir12)
dir32=norme(dir32)
print(fmt_str.format(288,'dir12',list2str(dir12)))
print(fmt_str.format(289,'dir32',list2str(dir32)))
n=prodve(dir12,dir32)
print(fmt_str.format(291,'n',list2str(n)))
isValid=True
if norm3(n) < 1.e-15:
isValid = False
print(fmt_str.format(294,'isValid',isValid))
else:
n=norme(n)
print(fmt_str.format(297,'n',list2str(n)))
if (fabs(n[0]) < 1.e-5 and fabs(n[1]) < 1.e-5 and fabs(n[2]) <
1.e-5):
isValid = False
print(fmt_str.format(299,'isValid',isValid))
if not isValid:
print('n = norme(c->Circle.n) unable to simulate; set to
[0,0,1]')
n=[0.0,0.0,1.0]
m=prodve(n,dir12)
print(fmt_str.format(310,'m',list2str(m)))
m=norme(m)
print(fmt_str.format(311,'m',list2str(m)))
mat=[dir12,m,n]
print(fmt_str.format(321,'mat',mat2str(mat)))
if n[0]==0.0 and n[1]==0.0 and False:
mat=[[1,0,0],[0,1,0],[0,0,1]]
print(fmt_str.format(333,'mat',mat2str(mat)))
v0=projette(v0,mat)
v2=projette(v2,mat)
v3=projette(v3,mat)
print(fmt_str.format(336,'v0',list2str(v0)))
print(fmt_str.format(337,'v2',list2str(v2)))
print(fmt_str.format(339,'v3',list2str(v3)))
R=sqrt(v0[0] * v0[0] + v0[1] * v0[1])
R2=sqrt(v2[0] * v2[0] + v2[1] * v2[1])
print(fmt_str.format(340,'R',R))
print(fmt_str.format(341,'R2',R2))
A4 = myatan2(v3[1],v3[0]);
print(fmt_str.format(358,'A4',A4))
A4 = angle_02pi(A4);
print(fmt_str.format(359,'A4',A4))
x1 = v0[0] * cos(A4) + v0[1] * sin(A4);
y1 = -v0[0] * sin(A4) + v0[1] * cos(A4);
x3 = v2[0] * cos(A4) + v2[1] * sin(A4);
y3 = -v2[0] * sin(A4) + v2[1] * cos(A4);
print(fmt_str.format(360,'x1',x1))
print(fmt_str.format(361,'y1',y1))
print(fmt_str.format(362,'x3',x3))
print(fmt_str.format(363,'y3',y3))
sys=[[x1*x1,y1*y1],[x3*x3,y3*y3]]
print(fmt_str.format(368,'sys',mat2str(sys)))
sol=sys2x2(sys,[1.0,1.0])
print(fmt_str.format(371,'sol',list2str(sol)))
if sol[0] <= 0 or sol[1] <= 0:
print('373: Ellipse is wrong')
A1=0.0
A3=0.0
f1=R
f2=R
print(fmt_str.format(374,'A1,A3',0.0))
print(fmt_str.format(375,'f1,f2',R))
else:
f1 = sqrt(1.0 / sol[0])
f2 = sqrt(1.0 / sol[1])
print(fmt_str.format(378,'f1',f1))
print(fmt_str.format(379,'f2',f2))
if x1 < 0:
A1 = -myasin(y1 / f2) + A4 + pi
print(fmt_str.format(384,'A1',A1))
else:
A1 = myasin(y1 / f2) + A4
print(fmt_str.format(386,'A1',A1))
if(x3 < 0):
A3 = -myasin(y3 / f2) + A4 + pi
print(fmt_str.format(388,'A3',A3))
else:
A3 = myasin(y3 / f2) + A4
print(fmt_str.format(390,'A3',A3))
A1 = angle_02pi(A1)
A3 = angle_02pi(A3)
print(fmt_str.format(399,'A1',A1))
print(fmt_str.format(400,'A3',A3))
if A1>=A3:
A3+=2*pi
print(fmt_str.format(402,'A3',A3))
if __name__=="__main__":
rightIntersection=[10.,5.,0.]
leftIntersection=[-9.99552,5.,0.]
center=[0.,5.,0.]
pmaj=[0.,20.,0.]
nmaj=[0.,-10.,0.]
p8=[-6.54654,-6.33893,0.]
p9=[6.54654,-6.33893,0.]
p4=leftIntersection
p3=rightIntersection
p5=center
p6=pmaj
p7=nmaj
trace(p8,p5,p9,p7)
example.geo
Description: Binary data
This is the output of "python GmshTest.py" 269: dir12= [-6.54654,-11.338930000000001,0.0] 270: dir32= [6.54654,-11.338930000000001,0.0] 272: dir42= [0.0,-15.0,0.0] 288: dir12= [-0.5000003272060961,-0.8660252148718285,0.0] 289: dir32= [0.5000003272060961,-0.8660252148718285,0.0] 291: n = [0.0,0.0,0.8660257816092879] 297: n = [0.0,0.0,1.0] 310: m = [0.8660252148718285,-0.5000003272060961,0.0] 311: m = [0.8660252148718286,-0.5000003272060962,0.0] 321: mat = [ [-0.5000003272060961,-0.8660252148718285,0.0] [0.8660252148718286,-0.5000003272060962,0.0] [0.0,0.0,1.0] ] 336: v0 = [13.093071431734419,0.0,0.0] 337: v2 = [6.546527147598827,11.338937420334043,0.0] 339: v3 = [12.990378223077427,7.5000049080914435,0.0] 340: R = 13.093071431734419 341: R2 = 13.09307143173442 358: A4 = 0.5235991534233956 359: A4 = 0.5235991534233956 360: x1 = 11.33893 361: y1 = -6.54654 362: x3 = 11.338930000000003 363: y3 = 6.546540000000002 368: sys = [ [128.5713335449,42.8571859716] [128.57133354490006,42.857185971600025] ] sys2x2: norm=36734.652397836486,det=0.0,det/norm=0.0 371: sol = [0.0,0.0] 373: Ellipse is wrong 374: A1,A3= 0.0 375: f1,f2= 13.093071431734419 399: A1 = 0.0 400: A3 = 0.0 402: A3 = 6.283185307179586
_______________________________________________ gmsh mailing list [email protected] http://onelab.info/mailman/listinfo/gmsh
