source: issm/trunk/src/m/extrusion/project2d.py@ 17806

import numpy as npy

def project2d(md3d,value,layer):
        '''
        returns the value of a field for a given layer of the mesh
        
        project 'value' vector taken at layer 'layer' from extruded 2d-3d
        mesh onto 2d mesh used to do the extrusion. This routine is used to
        compare values between a 2d-3d mesh at a certain layer, and the
        equivalent value (if it exists), on the original 2d mesh.       

        Note: for consistency with usage in ISSM, layers are indexed starting
        from one.  This routine accounts for the fact that python indexes from
        zero.  In other words, when using this function consider the bottom
        layer as layer=1.  

   Usage:
      projection_value=project2d(md3d,value,layer)

   Example:
      vel2=project2d(md3d,md3d.vel,2)
        '''

        if md3d.mesh.domaintype().lower() != '3d':
                raise StandardError("model passed to project2d function should be 3D")

        if layer<1 or layer>md3d.mesh.numberoflayers:
                raise ValueError("layer must be between 0 and %i" % md3d.mesh.numberoflayers)
        
        # coerce to array in case float is passed
        if type(value)!=npy.ndarray:
                print 'coercing array'
                value=npy.array(value)

        vec2d=False
        if value.ndim==2: 
                value=value.reshape(-1,)
                vec2d=True
        
        if value.size==1:
                projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d]
        elif value.shape[0]==md3d.mesh.numberofvertices:
                #print 'indices: ', (layer-1)*md3d.mesh.numberofvertices2d, layer*md3d.mesh.numberofvertices2d
                projection_value=value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d]
        elif value.shape[0]==md3d.mesh.numberofvertices+1:
                projection_value=[value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d], value[-1]]
        else:
                projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d]

        if vec2d:
                projection_value=projection_value.reshape(-1,1)

        return projection_value