[Gmsh] Ellipse do not work when start and end points are symmetric wrt major axis
Shengjie Xu
xsj617603321 at gmail.com
Sat Jun 9 20:42:28 CEST 2018
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://onelab.info/pipermail/gmsh/attachments/20180609/b563dacf/attachment.html>
-------------- next part --------------
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)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: example.geo
Type: application/octet-stream
Size: 322 bytes
Desc: not available
URL: <http://onelab.info/pipermail/gmsh/attachments/20180609/b563dacf/attachment.geo>
-------------- next part --------------
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
More information about the gmsh
mailing list