source: issm/trunk-jpl/src/m/model/mesh/triangle.py@ 12112

from numpy import *
import TriMesh as tm
import NodeConnectivity as nc
import ElementConnectivity as ec

def triangle(md, domainname, resolution,riftname=''):
        #TRIANGLE - create model mesh using the triangle package
        #
        #   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
        #   where md is a @model object, domainname is the name of an Argus domain outline file, 
        #   and resolution is a characteristic length for the mesh (same unit as the domain outline
        #   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
        #
        #   Usage:
        #      md=triangle(md,domainname,resolution)
        #   or md=triangle(md,domainname, resolution, riftname)
        #
        #   Examples:
        #      md=triangle(md,'DomainOutline.exp',1000);
        #      md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp');


        #Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would 
        #be made of 1000*1000 area squares). 

        #Check that mesh was not already run, and warn user: 
        if md.mesh.numberofelements != 0.:
                choice = input('This model already has a mesh. Are you sure you want to go ahead? (y/n)')
                if choice != 'y':
                        print 'no meshing done ... exiting'
                        return []
                
        area = resolution**2.

        #Mesh using TriMesh
        [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=tm.TriMesh(domainname,riftname,area)

        #Fill in rest of fields:
        md.mesh.numberofelements = size(md.mesh.elements)
        md.mesh.numberofvertices = size(md.mesh.x)
        md.mesh.z = zeros(md.mesh.numberofvertices)
        md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices)
        md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)] = 1.
        md.mesh.vertexonbed = ones(md.mesh.numberofvertices)
        md.mesh.vertexonsurface = ones(md.mesh.numberofvertices)
        md.mesh.elementonbed = ones(md.mesh.numberofelements)
        md.mesh.elementonsurface = ones(md.mesh.numberofelements)

        #Now, build the connectivity tables for this mesh.
        [md.mesh.vertexconnectivity]= nc.NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
        [md.mesh.elementconnectivity] = ec.ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)

        #type of model
        md.mesh.dimension = 2.
        return md