source: issm/trunk-jpl/src/py3/partition/AreaAverageOntoPartition.py@ 23677

import numpy as np
import copy

def AreaAverageOntoPartition(md,vector,layer=None):
        '''AREAAVERAGEONTOPARTITION 
   compute partition values for a certain vector expressed on the vertices of the mesh.
   Use area weighted average.

   Usage:
      average=AreaAverageOntoPartition(md,vector)
      average=AreaAverageOntoPartition(md,vector,layer) #if in 3D, chose which layer is partitioned
'''
        #some checks
        if(md.mesh.dimension()==3):
                if layer == None:
                        raise RuntimeError('AreaAverageOntoPartition: layer should be provided onto which Area Averaging occurs')

                #save 3D model
                md3d = copy.deepcopy(md)

                md.mesh.elements = md.mesh.elements2d
                md.mesh.x = md.mesh.x2d
                md.mesh.y = md.mesh.y2d
                md.mesh.numberofvertices = md.mesh.numberofvertices2d
                md.mesh.numberofelements = md.mesh.numberofelements2d
                md.qmu.vertex_weight = []
                md.mesh.vertexconnectivity = []

                #run connectivity routine
                md = adjacency(md)

                #finally, project vector: 
                vector = project2d(md3d,vector,layer)
                md.qmu.partition = project2d(md3d,md3d.qmu.partition,layer)


        #ok, first check that part is Matlab indexed
        part = (md.qmu.partition).copy()
        part = part.flatten() + 1

        #some check: 
        if md.qmu.numberofpartitions != max(part):
                raise RuntimeError('AreaAverageOntoPartition error message: ''npart'' should be equal to max(md.qmu.partition)')


        #initialize output
        partvector = np.zeros((max(part)))

        #start weight average
        weightedvector = vector.flatten()*md.qmu.vertex_weight
        for i in range(max(part)):
                pos=np.where((part-1)==i)
                partvector[i]=sum(weightedvector[pos])/sum(md.qmu.vertex_weight[pos])


        #in 3D, restore 3D model:
        if(md.mesh.dimension()==3):
                md = copy.deepcopy(md3d)


        return partvector