Changeset 13520
- Timestamp:
- 10/03/12 10:38:28 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/mesh/triangle.py
r13009 r13520 1 from numpy import * 2 import TriMesh as tm 3 import NodeConnectivity as nc 4 import ElementConnectivity as ec 1 import numpy 2 from TriMesh import * 3 from NodeConnectivity import * 4 from ElementConnectivity import * 5 from MatlabFuncs import * 5 6 6 def triangle(md, domainname, resolution,riftname=''): 7 #TRIANGLE - create model mesh using the triangle package 8 # 9 # This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution 10 # where md is a @model object, domainname is the name of an Argus domain outline file, 11 # and resolution is a characteristic length for the mesh (same unit as the domain outline 12 # unit). Riftname is an optional argument (Argus domain outline) describing rifts. 13 # 14 # Usage: 15 # md=triangle(md,domainname,resolution) 16 # or md=triangle(md,domainname, resolution, riftname) 17 # 18 # Examples: 19 # md=triangle(md,'DomainOutline.exp',1000); 20 # md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp'); 7 def triangle(md,domainname,*args): 8 """ 9 TRIANGLE - create model mesh using the triangle package 21 10 11 This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution 12 where md is a @model object, domainname is the name of an Argus domain outline file, 13 and resolution is a characteristic length for the mesh (same unit as the domain outline 14 unit). Riftname is an optional argument (Argus domain outline) describing rifts. 15 16 Usage: 17 md=triangle(md,domainname,resolution) 18 or md=triangle(md,domainname, resolution, riftname) 19 20 Examples: 21 md=triangle(md,'DomainOutline.exp',1000); 22 md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp'); 23 """ 22 24 23 25 #Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m resolution node would 24 26 #be made of 1000*1000 area squares). 25 27 28 if len(args)==1: 29 resolution=args[0] 30 riftname='' 31 if len(args)==2: 32 riftname=args[0] 33 resolution=args[1] 34 26 35 #Check that mesh was not already run, and warn user: 27 if md.mesh.numberofelements != 0.:28 choice = input('This model already has a mesh. Are you sure you want to go ahead? (y/n)')29 if choice != 'y':36 if md.mesh.numberofelements: 37 choice = raw_input('This model already has a mesh. Are you sure you want to go ahead? (y/n)') 38 if not strcmp(choice,'y'): 30 39 print 'no meshing done ... exiting' 31 return []32 33 area = resolution**2 .40 return None 41 42 area = resolution**2 34 43 35 44 #Mesh using TriMesh 36 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=tm.TriMesh(domainname,riftname,area) 37 45 [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area) 38 46 39 47 #Fill in rest of fields: 40 md.mesh.numberofelements = len(md.mesh.elements)41 md.mesh.numberofvertices = len(md.mesh.x)42 md.mesh.z = zeros(md.mesh.numberofvertices)43 md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices)48 md.mesh.numberofelements = numpy.size(md.mesh.elements,axis=0) 49 md.mesh.numberofvertices = numpy.size(md.mesh.x) 50 md.mesh.z = numpy.zeros(md.mesh.numberofvertices) 51 md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices) 44 52 md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1. 45 md.mesh.vertexonbed = ones(md.mesh.numberofvertices)46 md.mesh.vertexonsurface = ones(md.mesh.numberofvertices)47 md.mesh.elementonbed = ones(md.mesh.numberofelements)48 md.mesh.elementonsurface = ones(md.mesh.numberofelements)53 md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices) 54 md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices) 55 md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements) 56 md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements) 49 57 50 58 #Now, build the connectivity tables for this mesh. 51 [md.mesh.vertexconnectivity]= nc.NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)52 [md.mesh.elementconnectivity] = ec.ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)53 59 [md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) 60 [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) 61 54 62 #type of model 55 63 md.mesh.dimension = 2. 64 56 65 return md 66
Note:
See TracChangeset
for help on using the changeset viewer.