Changeset 20223 for issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py

Timestamp:

02/19/16 16:26:56 (9 years ago)

Author:

dlcheng

Message:

ADD: mesh3dsurface python version. CHG: gmshplanet python fscanf translations and verification

File:

• issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (modified) (6 diffs)

issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py

@@ -4 +4 @@
 from numpy import *
 from pairoptions import *
+from mesh3dsurface import *
+import subprocess
 
 def gmshplanet(*varargin):
@@ -30 +32 @@
 
         #initialize mesh: 
-        #mesh=mesh3dsurface()
-        mesh=mesh2d()
+        mesh=mesh3dsurface()
+        
 
         #create .geo file:  {{{
@@ -44 +46 @@
         fid.write('resolution=%g;\n'%resolution)
         fid.write('radius=%g;\n'%radius)
-        fid.write('Point(1) = [0.0,0.0,0.0,resolution];\n')
-        fid.write('Point(2) = [radius,0.0,0.0,resolution];\n')
-        fid.write('Point(3) = [0,radius,0.0,resolution];\n')
-        fid.write('Circle(1) = [2,1,3];\n')
-        fid.write('Point(4) = [-radius,0,0.0,resolution];\n')
-        fid.write('Point(5) = [0,-radius,0.0,resolution];\n')
-        fid.write('Circle(2) = [3,1,4];\n')
-        fid.write('Circle(3) = [4,1,5];\n')
-        fid.write('Circle(4) = [5,1,2];\n')
-        fid.write('Point(6) = [0,0,-radius,resolution];\n')
-        fid.write('Point(7) = [0,0,radius,resolution];\n')
-        fid.write('Circle(5) = [3,1,6];\n')
-        fid.write('Circle(6) = [6,1,5];\n')
-        fid.write('Circle(7) = [5,1,7];\n')
-        fid.write('Circle(8) = [7,1,3];\n')
-        fid.write('Circle(9) = [2,1,7];\n')
-        fid.write('Circle(10) = [7,1,4];\n')
-        fid.write('Circle(11) = [4,1,6];\n')
-        fid.write('Circle(12) = [6,1,2];\n')
-        fid.write('Line Loop(13) = [2,8,-10];\n')
-        fid.write('Ruled Surface(14) = [13];\n')
-        fid.write('Line Loop(15) = [10,3,7];\n')
-        fid.write('Ruled Surface(16) = [15];\n')
-        fid.write('Line Loop(17) = [-8,-9,1];\n')
-        fid.write('Ruled Surface(18) = [17];\n')
-        fid.write('Line Loop(19) = [-11,-2,5];\n')
-        fid.write('Ruled Surface(20) = [19];\n')
-        fid.write('Line Loop(21) = [-5,-12,-1];\n')
-        fid.write('Ruled Surface(22) = [21];\n')
-        fid.write('Line Loop(23) = [-3,11,6];\n')
-        fid.write('Ruled Surface(24) = [23];\n')
-        fid.write('Line Loop(25) = [-7,4,9];\n')
-        fid.write('Ruled Surface(26) = [25];\n')
-        fid.write('Line Loop(27) = [-4,12,-6];\n')
-        fid.write('Ruled Surface(28) = [27];\n')
-        fid.write('Surface Loop(29) = [28,26,16,14,20,24,22,18];\n')
-        fid.write('Volume(30) = [29];\n')
-        fid.write('Physical Surface(1) = [28,26,16,14,20,24,22,18];\n')
+        fid.write('Point(1) = {0.0,0.0,0.0,resolution};\n')
+        fid.write('Point(2) = {radius,0.0,0.0,resolution};\n')
+        fid.write('Point(3) = {0,radius,0.0,resolution};\n')
+        fid.write('Circle(1) = {2,1,3};\n')
+        fid.write('Point(4) = {-radius,0,0.0,resolution};\n')
+        fid.write('Point(5) = {0,-radius,0.0,resolution};\n')
+        fid.write('Circle(2) = {3,1,4};\n')
+        fid.write('Circle(3) = {4,1,5};\n')
+        fid.write('Circle(4) = {5,1,2};\n')
+        fid.write('Point(6) = {0,0,-radius,resolution};\n')
+        fid.write('Point(7) = {0,0,radius,resolution};\n')
+        fid.write('Circle(5) = {3,1,6};\n')
+        fid.write('Circle(6) = {6,1,5};\n')
+        fid.write('Circle(7) = {5,1,7};\n')
+        fid.write('Circle(8) = {7,1,3};\n')
+        fid.write('Circle(9) = {2,1,7};\n')
+        fid.write('Circle(10) = {7,1,4};\n')
+        fid.write('Circle(11) = {4,1,6};\n')
+        fid.write('Circle(12) = {6,1,2};\n')
+        fid.write('Line Loop(13) = {2,8,-10};\n')
+        fid.write('Ruled Surface(14) = {13};\n')
+        fid.write('Line Loop(15) = {10,3,7};\n')
+        fid.write('Ruled Surface(16) = {15};\n')
+        fid.write('Line Loop(17) = {-8,-9,1};\n')
+        fid.write('Ruled Surface(18) = {17};\n')
+        fid.write('Line Loop(19) = {-11,-2,5};\n')
+        fid.write('Ruled Surface(20) = {19};\n')
+        fid.write('Line Loop(21) = {-5,-12,-1};\n')
+        fid.write('Ruled Surface(22) = {21};\n')
+        fid.write('Line Loop(23) = {-3,11,6};\n')
+        fid.write('Ruled Surface(24) = {23};\n')
+        fid.write('Line Loop(25) = {-7,4,9};\n')
+        fid.write('Ruled Surface(26) = {25};\n')
+        fid.write('Line Loop(27) = {-4,12,-6};\n')
+        fid.write('Ruled Surface(28) = {27};\n')
+        fid.write('Surface Loop(29) = {28,26,16,14,20,24,22,18};\n')
+        fid.write('Volume(30) = {29};\n')
+        fid.write('Physical Surface(1) = {28,26,16,14,20,24,22,18};\n')
         fid.write('Physical Volume(2) = 30;\n')
         fid.close()
@@ -108 +110 @@
         #call gmsh
         if options.exist('refine'):
-                eval('!gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos')
+                subprocess.call('gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos',shell=True)
         else:
                 #call gmsh
-                eval('!gmsh -tol 1e-8 -2 sphere.geo')
+                subprocess.call('gmsh -tol 1e-8 -2 sphere.geo',shell=True)
 
         #import mesh:  {{{
@@ -117 +119 @@
 
         #Get Mesh format
-        A=fscanf(fid,'%s',1)
+        A=fid.readline().strip()
         if not strcmp(A,'$MeshFormat'):
                 raise RuntimeError(['Expecting $MeshFormat (', A, ')'])
 
-        A=fscanf(fid,'%f %i %i',[1, 3])
-        A=fscanf(fid,'%s',1)
+        A=fid.readline().split()
+        A=fid.readline().strip()
         if not strcmp(A,'$EndMeshFormat'):
                 raise RuntimeError(['Expecting $EndMeshFormat (', A, ')'])
 
         #Nodes
-        A=fscanf(fid,'%s',1)
+        A=fid.readline().strip()
         if not strcmp(A,'$Nodes'):
                 raise RuntimeError(['Expecting $Nodes (', A, ')'])
 
-        mesh.numberofvertices=fscanf(fid,'%i',1)
-        A=fscanf(fid,'%i %f %f %f',[4, mesh.numberofvertices])
-        mesh.x = transpose(A[1,:])
-        mesh.y = transpose(A[2,:])
-        mesh.z = transpose(A[3,:])
+        mesh.numberofvertices=int(fid.readline().strip())
+        mesh.x=empty(mesh.numberofvertices)
+        mesh.y=empty(mesh.numberofvertices)
+        mesh.z=empty(mesh.numberofvertices)
+        for i in range(mesh.numberofvertices):
+                A=fid.readline().split()
+                mesh.x[i]=float(A[0])
+                mesh.y[i]=float(A[1])
+                mesh.z[i]=float(A[2])
 
-        A=fscanf(fid,'%s',1)
+        A=fid.readline().strip()
         if not strcmp(A,'$EndNodes'):
                 raise RuntimeError(['Expecting $EndNodes (', A, ')'])
 
         #Elements
-        A=fscanf(fid,'%s',1)
+        A=fid.readline().strip()
         if not strcmp(A,'$Elements'):
                 raise RuntimeError(['Expecting $Elements (', A, ')'])
-        mesh.numberofelements=fscanf(fid,'%i',1)
-        A=fscanf(fid,'%i %i %i %i %i %i %i %i',[8, mesh.numberofelements])
-        mesh.elements=transpose(A[6:8,:]) 
-        A=fscanf(fid,'%s',1)
+        mesh.numberofelements=int(fid.readline().strip())
+        mesh.elements=empty([mesh.numberofelements,3])
+        for i in range(mesh.numberofelements):
+                A=fid.readline().split()
+                mesh.elements[i]=[float(A[5]),float(A[6]),float(A[7])]
+        A=fid.readline().strip()
         if not strcmp(A,'$EndElements'):
                 raise RuntimeError(['Expecting $EndElements (', A, ')'])
@@ -156 +164 @@
         #figure out other fields in mesh3dsurface: 
         mesh.r=sqrt(mesh.x**2+mesh.y**2+mesh.z**2)
-        mesh.lat = asin(mesh.z/mesh.r)/pi*180
-        mesh.long = atan2(mesh.y,mesh.x)/pi*180
+        mesh.lat = arcsin(mesh.z/mesh.r)/pi*180
+        mesh.long = arctan2(mesh.y,mesh.x)/pi*180
 
         #erase files: 
-        eval('!rm -rf sphere.geo sphere.msh sphere.pos')
-
+        subprocess.call('rm -rf sphere.geo sphere.msh sphere.pos',shell=True)
 
         #return mesh: