![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 20177
- Timestamp:
- 02/16/16 14:40:16 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py
r20173 r20177 30 30 31 31 #initialize mesh: 32 #mesh=mesh3dsurface() 32 33 mesh=mesh2d() 33 34 … … 35 36 fid=open('sphere.geo','w') 36 37 37 fid.write('Mesh.Algorithm = 1 \n')38 fid.write('Mesh.Algorithm = 1;\n') 38 39 if options.exist('refine'): 39 fid.write('Mesh.Algorithm = 7 \n')40 fid.write('Mesh.CharacteristicLengthFromPoints= 0 \n')41 fid.write('Mesh.SmoothRatio= 3 \n')42 fid.write('Mesh.RemeshAlgorithm= 1 \n')43 fid.write('resolution= #g\n',resolution)44 fid.write('radius= #g\n',radius)45 fid.write('Point(1) = [0.0,0.0,0.0,resolution] \n')46 fid.write('Point(2) = [radius,0.0,0.0,resolution] \n')47 fid.write('Point(3) = [0,radius,0.0,resolution] \n')48 fid.write('Circle(1) = [2,1,3] \n')49 fid.write('Point(4) = [-radius,0,0.0,resolution] \n')50 fid.write('Point(5) = [0,-radius,0.0,resolution] \n')51 fid.write('Circle(2) = [3,1,4] \n')52 fid.write('Circle(3) = [4,1,5] \n')53 fid.write('Circle(4) = [5,1,2] \n')54 fid.write('Point(6) = [0,0,-radius,resolution] \n')55 fid.write('Point(7) = [0,0,radius,resolution] \n')56 fid.write('Circle(5) = [3,1,6] \n')57 fid.write('Circle(6) = [6,1,5] \n')58 fid.write('Circle(7) = [5,1,7] \n')59 fid.write('Circle(8) = [7,1,3] \n')60 fid.write('Circle(9) = [2,1,7] \n')61 fid.write('Circle(10) = [7,1,4] \n')62 fid.write('Circle(11) = [4,1,6] \n')63 fid.write('Circle(12) = [6,1,2] \n')64 fid.write('Line Loop(13) = [2,8,-10] \n')65 fid.write('Ruled Surface(14) = [13] \n')66 fid.write('Line Loop(15) = [10,3,7] \n')67 fid.write('Ruled Surface(16) = [15] \n')68 fid.write('Line Loop(17) = [-8,-9,1] \n')69 fid.write('Ruled Surface(18) = [17] \n')70 fid.write('Line Loop(19) = [-11,-2,5] \n')71 fid.write('Ruled Surface(20) = [19] \n')72 fid.write('Line Loop(21) = [-5,-12,-1] \n')73 fid.write('Ruled Surface(22) = [21] \n')74 fid.write('Line Loop(23) = [-3,11,6] \n')75 fid.write('Ruled Surface(24) = [23] \n')76 fid.write('Line Loop(25) = [-7,4,9] \n')77 fid.write('Ruled Surface(26) = [25] \n')78 fid.write('Line Loop(27) = [-4,12,-6] \n')79 fid.write('Ruled Surface(28) = [27] \n')80 fid.write('Surface Loop(29) = [28,26,16,14,20,24,22,18] \n')81 fid.write('Volume(30) = [29] \n')82 fid.write('Physical Surface(1) = [28,26,16,14,20,24,22,18] \n')83 fid.write('Physical Volume(2) = 30 \n')40 fid.write('Mesh.Algorithm = 7;\n') 41 fid.write('Mesh.CharacteristicLengthFromPoints= 0;\n') 42 fid.write('Mesh.SmoothRatio= 3;\n') 43 fid.write('Mesh.RemeshAlgorithm= 1;\n') 44 fid.write('resolution=%g;\n'%resolution) 45 fid.write('radius=%g;\n'%radius) 46 fid.write('Point(1) = [0.0,0.0,0.0,resolution];\n') 47 fid.write('Point(2) = [radius,0.0,0.0,resolution];\n') 48 fid.write('Point(3) = [0,radius,0.0,resolution];\n') 49 fid.write('Circle(1) = [2,1,3];\n') 50 fid.write('Point(4) = [-radius,0,0.0,resolution];\n') 51 fid.write('Point(5) = [0,-radius,0.0,resolution];\n') 52 fid.write('Circle(2) = [3,1,4];\n') 53 fid.write('Circle(3) = [4,1,5];\n') 54 fid.write('Circle(4) = [5,1,2];\n') 55 fid.write('Point(6) = [0,0,-radius,resolution];\n') 56 fid.write('Point(7) = [0,0,radius,resolution];\n') 57 fid.write('Circle(5) = [3,1,6];\n') 58 fid.write('Circle(6) = [6,1,5];\n') 59 fid.write('Circle(7) = [5,1,7];\n') 60 fid.write('Circle(8) = [7,1,3];\n') 61 fid.write('Circle(9) = [2,1,7];\n') 62 fid.write('Circle(10) = [7,1,4];\n') 63 fid.write('Circle(11) = [4,1,6];\n') 64 fid.write('Circle(12) = [6,1,2];\n') 65 fid.write('Line Loop(13) = [2,8,-10];\n') 66 fid.write('Ruled Surface(14) = [13];\n') 67 fid.write('Line Loop(15) = [10,3,7];\n') 68 fid.write('Ruled Surface(16) = [15];\n') 69 fid.write('Line Loop(17) = [-8,-9,1];\n') 70 fid.write('Ruled Surface(18) = [17];\n') 71 fid.write('Line Loop(19) = [-11,-2,5];\n') 72 fid.write('Ruled Surface(20) = [19];\n') 73 fid.write('Line Loop(21) = [-5,-12,-1];\n') 74 fid.write('Ruled Surface(22) = [21];\n') 75 fid.write('Line Loop(23) = [-3,11,6];\n') 76 fid.write('Ruled Surface(24) = [23];\n') 77 fid.write('Line Loop(25) = [-7,4,9];\n') 78 fid.write('Ruled Surface(26) = [25];\n') 79 fid.write('Line Loop(27) = [-4,12,-6];\n') 80 fid.write('Ruled Surface(28) = [27];\n') 81 fid.write('Surface Loop(29) = [28,26,16,14,20,24,22,18];\n') 82 fid.write('Volume(30) = [29];\n') 83 fid.write('Physical Surface(1) = [28,26,16,14,20,24,22,18];\n') 84 fid.write('Physical Volume(2) = 30;\n') 84 85 fid.close() 85 86 #}}} … … 92 93 fid=open('sphere.pos','w') 93 94 94 fid.write('View "background mesh" [ \n')95 fid.write('View "background mesh" [;\n') 95 96 for i in range(meshini.numberofelements): 96 fid.write('ST( #g,#g,#g,#g,#g,#g,#g,#g,#g)[#g,#g,#g]\n',\97 fid.write('ST(%g,%g,%g,%g,%g,%g,%g,%g,%g)[%g,%g,%g];\n',\ 97 98 meshini.x(meshini.elements(i,0)), meshini.y(meshini.elements(i,0)), meshini.z(meshini.elements(i,0)),\ 98 99 meshini.x(meshini.elements(i,1)), meshini.y(meshini.elements(i,1)), meshini.z(meshini.elements(i,1)),\ … … 100 101 metric(meshini.elements(i,0)), metric(meshini.elements(i,1)), metric(meshini.elements(i,2))\ 101 102 ) 102 fid.write('] \n')103 fid.write('];\n') 103 104 104 105 fid.close() … … 107 108 #call gmsh 108 109 if options.exist('refine'): 109 eval( ['!gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos'])110 eval('!gmsh -tol 1e-8 -2 sphere.geo -bgm sphere.pos') 110 111 else: 111 112 #call gmsh 112 eval( ['!gmsh -tol 1e-8 -2 sphere.geo'])113 eval('!gmsh -tol 1e-8 -2 sphere.geo') 113 114 114 115 #import mesh: {{{ … … 116 117 117 118 #Get Mesh format 118 A=fscanf(fid,' #s',1)119 A=fscanf(fid,'%s',1) 119 120 if not strcmp(A,'$MeshFormat'): 120 121 raise RuntimeError(['Expecting $MeshFormat (', A, ')']) 121 122 122 A=fscanf(fid,' #f #i #i',[1, 3])123 A=fscanf(fid,' #s',1)123 A=fscanf(fid,'%f %i %i',[1, 3]) 124 A=fscanf(fid,'%s',1) 124 125 if not strcmp(A,'$EndMeshFormat'): 125 126 raise RuntimeError(['Expecting $EndMeshFormat (', A, ')']) 126 127 127 128 #Nodes 128 A=fscanf(fid,' #s',1)129 A=fscanf(fid,'%s',1) 129 130 if not strcmp(A,'$Nodes'): 130 131 raise RuntimeError(['Expecting $Nodes (', A, ')']) 131 132 132 mesh.numberofvertices=fscanf(fid,' #i',1)133 A=fscanf(fid,' #i #f #f #f',[4, mesh.numberofvertices])133 mesh.numberofvertices=fscanf(fid,'%i',1) 134 A=fscanf(fid,'%i %f %f %f',[4, mesh.numberofvertices]) 134 135 mesh.x = transpose(A[1,:]) 135 136 mesh.y = transpose(A[2,:]) 136 137 mesh.z = transpose(A[3,:]) 137 138 138 A=fscanf(fid,' #s',1)139 A=fscanf(fid,'%s',1) 139 140 if not strcmp(A,'$EndNodes'): 140 141 raise RuntimeError(['Expecting $EndNodes (', A, ')']) 141 142 142 143 #Elements 143 A=fscanf(fid,' #s',1)144 A=fscanf(fid,'%s',1) 144 145 if not strcmp(A,'$Elements'): 145 146 raise RuntimeError(['Expecting $Elements (', A, ')']) 146 mesh.numberofelements=fscanf(fid,' #i',1)147 A=fscanf(fid,' #i #i #i #i #i #i #i #i',[8, mesh.numberofelements])147 mesh.numberofelements=fscanf(fid,'%i',1) 148 A=fscanf(fid,'%i %i %i %i %i %i %i %i',[8, mesh.numberofelements]) 148 149 mesh.elements=transpose(A[6:8,:]) 149 A=fscanf(fid,' #s',1)150 A=fscanf(fid,'%s',1) 150 151 if not strcmp(A,'$EndElements'): 151 152 raise RuntimeError(['Expecting $EndElements (', A, ')']) … … 159 160 160 161 #erase files: 161 eval( ['!rm -rf sphere.geo sphere.msh sphere.pos'])162 eval('!rm -rf sphere.geo sphere.msh sphere.pos') 162 163 163 164
Note:
See TracChangeset
for help on using the changeset viewer.
