Ignore:
Timestamp:
06/07/17 10:50:54 (8 years ago)
Author:
Eric.Larour
Message:

CHG: merged branch back to trunk-jpl 21754.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/model.py

    r21097 r21759  
    11#module imports {{{
    2 import numpy
     2import numpy as np
    33import copy
    44import sys
    … …  
    4545from steadystate import steadystate
    4646from transient import transient
    47 from gia import gia
     47from giaivins import giaivins
    4848from autodiff import autodiff
    4949from inversion import inversion
    5050from outputdefinition import outputdefinition
    5151from qmu import qmu
     52from amr import amr
    5253from results import results
    5354from radaroverlay import radaroverlay
    … …  
    107108                self.levelset         = levelset()
    108109                self.calving          = calving()
    109                 self.gia              = gia()
     110                self.gia              = giaivins()
    110111
    111112                self.autodiff         = autodiff()
    112113                self.inversion        = inversion()
    113114                self.qmu              = qmu()
     115                self.amr                                         = amr()
    114116
    115117                self.results          = results()
    … …  
    150152                        'levelset',\
    151153                        'calving',\
    152                                                 'gia',\
     154                                        'gia',\
    153155                        'autodiff',\
    154156                        'inversion',\
    155157                        'qmu',\
     158                        'amr',\
    156159                        'outputdefinition',\
    157160                        'results',\
    … …  
    194197                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("inversion","[%s,%s]" % ("1x1",obj.inversion.__class__.__name__),"parameters for inverse methods"))
    195198                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("qmu","[%s,%s]" % ("1x1",obj.qmu.__class__.__name__),"dakota properties"))
     199                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("amr","[%s,%s]" % ("1x1",obj.amr.__class__.__name__),"adaptive mesh refinement properties"))
    196200                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("outputdefinition","[%s,%s]" % ("1x1",obj.outputdefinition.__class__.__name__),"output definition"))
    197201                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("results","[%s,%s]" % ("1x1",obj.results.__class__.__name__),"model results"))
    … …  
    231235                #get elements that are inside area
    232236                flag_elem=FlagElements(md1,area)
    233                 if not numpy.any(flag_elem):
     237                if not np.any(flag_elem):
    234238                        raise RuntimeError("extracted model is empty")
    235239
    236240                #kick out all elements with 3 dirichlets
    237                 spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    238                 spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
    239                 flag=numpy.ones(md1.mesh.numberofvertices)
     241                spc_elem=np.nonzero(np.logical_not(flag_elem))[0]
     242                spc_node=np.unique(md1.mesh.elements[spc_elem,:])-1
     243                flag=np.ones(md1.mesh.numberofvertices)
    240244                flag[spc_node]=0
    241                 pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]
     245                pos=np.nonzero(np.logical_not(np.sum(flag[md1.mesh.elements-1],axis=1)))[0]
    242246                flag_elem[pos]=0
    243247
    244248                #extracted elements and nodes lists
    245                 pos_elem=numpy.nonzero(flag_elem)[0]
    246                 pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
     249                pos_elem=np.nonzero(flag_elem)[0]
     250                pos_node=np.unique(md1.mesh.elements[pos_elem,:])-1
    247251
    248252                #keep track of some fields
    249253                numberofvertices1=md1.mesh.numberofvertices
    250254                numberofelements1=md1.mesh.numberofelements
    251                 numberofvertices2=numpy.size(pos_node)
    252                 numberofelements2=numpy.size(pos_elem)
    253                 flag_node=numpy.zeros(numberofvertices1)
     255                numberofvertices2=np.size(pos_node)
     256                numberofelements2=np.size(pos_elem)
     257                flag_node=np.zeros(numberofvertices1)
    254258                flag_node[pos_node]=1
    255259
    256260                #Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
    257                 Pelem=numpy.zeros(numberofelements1,int)
    258                 Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
    259                 Pnode=numpy.zeros(numberofvertices1,int)
    260                 Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
     261                Pelem=np.zeros(numberofelements1,int)
     262                Pelem[pos_elem]=np.arange(1,numberofelements2+1)
     263                Pnode=np.zeros(numberofvertices1,int)
     264                Pnode[pos_node]=np.arange(1,numberofvertices2+1)
    261265
    262266                #renumber the elements (some node won't exist anymore)
    … …  
    283287                        #get field
    284288                        field=getattr(md1,fieldi)
    285                         fieldsize=numpy.shape(field)
     289                        fieldsize=np.shape(field)
    286290                        if hasattr(field,'__dict__') and not m.ismember(fieldi,['results'])[0]:    #recursive call
    287291                                object_fields=vars(field)
    … …  
    289293                                        #get field
    290294                                        field=getattr(getattr(md1,fieldi),fieldj)
    291                                         fieldsize=numpy.shape(field)
     295                                        fieldsize=np.shape(field)
    292296                                        if len(fieldsize):
    293297                                                #size = number of nodes * n
    … …  
    295299                                                        setattr(getattr(md2,fieldi),fieldj,field[pos_node])
    296300                                                elif fieldsize[0]==numberofvertices1+1:
    297                                                         setattr(getattr(md2,fieldi),fieldj,numpy.vstack((field[pos_node],field[-1,:])))
     301                                                        setattr(getattr(md2,fieldi),fieldj,np.vstack((field[pos_node],field[-1,:])))
    298302                                                #size = number of elements * n
    299303                                                elif fieldsize[0]==numberofelements1:
    … …  
    305309                                                setattr(md2,fieldi,field[pos_node])
    306310                                        elif fieldsize[0]==numberofvertices1+1:
    307                                                 setattr(md2,fieldi,numpy.hstack((field[pos_node],field[-1,:])))
     311                                                setattr(md2,fieldi,np.hstack((field[pos_node],field[-1,:])))
    308312                                        #size = number of elements * n
    309313                                        elif fieldsize[0]==numberofelements1:
    … …  
    320324                if md1.mesh.__class__.__name__=='mesh3dprisms':
    321325                        md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
    322                         pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]
    323                         md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1]
     326                        pos=np.where(~np.isnan(md2.mesh.uppervertex))[0]
     327                        md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1]
    324328
    325329                        md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node]
    326                         pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0]
    327                         md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1]
     330                        pos=np.where(~np.isnan(md2.mesh.lowervertex))[0]
     331                        md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1]
    328332
    329333                        md2.mesh.upperelements=md1.mesh.upperelements[pos_elem]
    330                         pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0]
    331                         md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1]
     334                        pos=np.where(~np.isnan(md2.mesh.upperelements))[0]
     335                        md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1]
    332336
    333337                        md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem]
    334                         pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0]
    335                         md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1]
     338                        pos=np.where(~np.isnan(md2.mesh.lowerelements))[0]
     339                        md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1]
    336340
    337341                #Initial 2d mesh
    338342                if md1.mesh.__class__.__name__=='mesh3dprisms':
    339                         flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)]
    340                         pos_elem_2d=numpy.nonzero(flag_elem_2d)[0]
    341                         flag_node_2d=flag_node[numpy.arange(0,md1.mesh.numberofvertices2d)]
    342                         pos_node_2d=numpy.nonzero(flag_node_2d)[0]
    343 
    344                         md2.mesh.numberofelements2d=numpy.size(pos_elem_2d)
    345                         md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
     343                        flag_elem_2d=flag_elem[np.arange(0,md1.mesh.numberofelements2d)]
     344                        pos_elem_2d=np.nonzero(flag_elem_2d)[0]
     345                        flag_node_2d=flag_node[np.arange(0,md1.mesh.numberofvertices2d)]
     346                        pos_node_2d=np.nonzero(flag_node_2d)[0]
     347
     348                        md2.mesh.numberofelements2d=np.size(pos_elem_2d)
     349                        md2.mesh.numberofvertices2d=np.size(pos_node_2d)
    346350                        md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
    347351                        md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1]
    … …  
    354358                #Edges
    355359                if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
    356                         if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some numpy.nans...
     360                        if np.ndim(md2.mesh.edges)>1 and np.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some np.nans...
    357361                                #renumber first two columns
    358                                 pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
     362                                pos=np.nonzero(md2.mesh.edges[:,3]!=-1)[0]
    359363                                md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
    360364                                md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
    … …  
    362366                                md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
    363367                                #remove edges when the 2 vertices are not in the domain.
    364                                 md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
     368                                md2.mesh.edges=md2.mesh.edges[np.nonzero(np.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
    365369                                #Replace all zeros by -1 in the last two columns
    366                                 pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0]
     370                                pos=np.nonzero(md2.mesh.edges[:,2]==0)[0]
    367371                                md2.mesh.edges[pos,2]=-1
    368                                 pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0]
     372                                pos=np.nonzero(md2.mesh.edges[:,3]==0)[0]
    369373                                md2.mesh.edges[pos,3]=-1
    370374                                #Invert -1 on the third column with last column (Also invert first two columns!!)
    371                                 pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0]
     375                                pos=np.nonzero(md2.mesh.edges[:,2]==-1)[0]
    372376                                md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3]
    373377                                md2.mesh.edges[pos,3]=-1
    … …  
    376380                                md2.mesh.edges[pos,0]=values
    377381                                #Finally remove edges that do not belong to any element
    378                                 pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
    379                                 md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0)
     382                                pos=np.nonzero(np.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
     383                                md2.mesh.edges=np.delete(md2.mesh.edges,pos,axis=0)
    380384
    381385                #Penalties
    382                 if numpy.any(numpy.logical_not(numpy.isnan(md2.stressbalance.vertex_pairing))):
    383                         for i in xrange(numpy.size(md1.stressbalance.vertex_pairing,axis=0)):
     386                if np.any(np.logical_not(np.isnan(md2.stressbalance.vertex_pairing))):
     387                        for i in xrange(np.size(md1.stressbalance.vertex_pairing,axis=0)):
    384388                                md2.stressbalance.vertex_pairing[i,:]=Pnode[md1.stressbalance.vertex_pairing[i,:]]
    385                         md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[numpy.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:]
    386                 if numpy.any(numpy.logical_not(numpy.isnan(md2.masstransport.vertex_pairing))):
    387                         for i in xrange(numpy.size(md1.masstransport.vertex_pairing,axis=0)):
     389                        md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[np.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:]
     390                if np.any(np.logical_not(np.isnan(md2.masstransport.vertex_pairing))):
     391                        for i in xrange(np.size(md1.masstransport.vertex_pairing,axis=0)):
    388392                                md2.masstransport.vertex_pairing[i,:]=Pnode[md1.masstransport.vertex_pairing[i,:]]
    389                         md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[numpy.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:]
     393                        md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[np.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:]
    390394
    391395                #recreate segments
    … …  
    394398                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)[0]
    395399                        md2.mesh.segments=contourenvelope(md2)
    396                         md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
     400                        md2.mesh.vertexonboundary=np.zeros(numberofvertices2,bool)
    397401                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
    398402                else:
    … …  
    401405                        md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)[0]
    402406                        segments=contourenvelope(md2)
    403                         md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
     407                        md2.mesh.vertexonboundary=np.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
    404408                        md2.mesh.vertexonboundary[segments[:,0:2]-1]=True
    405                         md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
     409                        md2.mesh.vertexonboundary=np.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
    406410                        #Then do it for 3d as usual
    407411                        md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)[0]
    … …  
    410414                #Boundary conditions: Dirichlets on new boundary
    411415                #Catch the elements that have not been extracted
    412                 orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    413                 orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
     416                orphans_elem=np.nonzero(np.logical_not(flag_elem))[0]
     417                orphans_node=np.unique(md1.mesh.elements[orphans_elem,:])-1
    414418                #Figure out which node are on the boundary between md2 and md1
    415                 nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
     419                nodestoflag1=np.intersect1d(orphans_node,pos_node)
    416420                nodestoflag2=Pnode[nodestoflag1].astype(int)-1
    417                 if numpy.size(md1.stressbalance.spcvx)>1 and numpy.size(md1.stressbalance.spcvy)>2 and numpy.size(md1.stressbalance.spcvz)>2:
    418                         if numpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1:
     421                if np.size(md1.stressbalance.spcvx)>1 and np.size(md1.stressbalance.spcvy)>2 and np.size(md1.stressbalance.spcvz)>2:
     422                        if np.size(md1.inversion.vx_obs)>1 and np.size(md1.inversion.vy_obs)>1:
    419423                                md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2]
    420424                                md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
    421425                        else:
    422                                 md2.stressbalance.spcvx[nodestoflag2]=numpy.nan
    423                                 md2.stressbalance.spcvy[nodestoflag2]=numpy.nan
     426                                md2.stressbalance.spcvx[nodestoflag2]=np.nan
     427                                md2.stressbalance.spcvy[nodestoflag2]=np.nan
    424428                                print "\n!! extract warning: spc values should be checked !!\n\n"
    425429                        #put 0 for vz
    426430                        md2.stressbalance.spcvz[nodestoflag2]=0
    427                 if numpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))):
    428                         md2.thermal.spctemperature[nodestoflag2,0]=1
     431                if np.any(np.logical_not(np.isnan(md1.thermal.spctemperature))):
     432                        md2.thermal.spctemperature[nodestoflag2]=1
    429433
    430434                #Results fields
    … …  
    441445                                                        #get subfields
    442446                                                        for solutionsubfield,subfield in fieldi.__dict__.iteritems():
    443                                                                 if   numpy.size(subfield)==numberofvertices1:
     447                                                                if   np.size(subfield)==numberofvertices1:
    444448                                                                        setattr(fieldr,solutionsubfield,subfield[pos_node])
    445                                                                 elif numpy.size(subfield)==numberofelements1:
     449                                                                elif np.size(subfield)==numberofelements1:
    446450                                                                        setattr(fieldr,solutionsubfield,subfield[pos_elem])
    447451                                                                else:
    … …  
    455459                                                #get subfields
    456460                                                for solutionsubfield,subfield in field.__dict__.iteritems():
    457                                                         if   numpy.size(subfield)==numberofvertices1:
     461                                                        if   np.size(subfield)==numberofvertices1:
    458462                                                                setattr(fieldr,solutionsubfield,subfield[pos_node])
    459                                                         elif numpy.size(subfield)==numberofelements1:
     463                                                        elif np.size(subfield)==numberofelements1:
    460464                                                                setattr(fieldr,solutionsubfield,subfield[pos_elem])
    461465                                                        else:
    … …  
    510514                                raise TypeError("extrusionexponent must be >=0")
    511515                        numlayers=args[0]
    512                         extrusionlist=(numpy.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1]
     516                        extrusionlist=(np.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1]
    513517
    514518                elif len(args)==3:    #two polynomial laws
    … …  
    520524                                raise TypeError("lower and upper extrusionexponents must be >=0")
    521525
    522                         lowerextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2.
    523                         upperextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2.
    524                         extrusionlist=numpy.unique(numpy.concatenate((lowerextrusionlist,1.-upperextrusionlist)))
     526                        lowerextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2.
     527                        upperextrusionlist=(np.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2.
     528                        extrusionlist=np.unique(np.concatenate((lowerextrusionlist,1.-upperextrusionlist)))
    525529
    526530                if numlayers<2:
    … …  
    550554                md.mesh.extractedelements           = mesh2d.extractedelements
    551555
    552                 x3d=numpy.empty((0))
    553                 y3d=numpy.empty((0))
    554                 z3d=numpy.empty((0))    #the lower node is on the bed
     556                x3d=np.empty((0))
     557                y3d=np.empty((0))
     558                z3d=np.empty((0))    #the lower node is on the bed
    555559                thickness3d=md.geometry.thickness    #thickness and bed for these nodes
    556560                bed3d=md.geometry.base
    … …  
    558562                #Create the new layers
    559563                for i in xrange(numlayers):
    560                         x3d=numpy.concatenate((x3d,md.mesh.x))
    561                         y3d=numpy.concatenate((y3d,md.mesh.y))
     564                        x3d=np.concatenate((x3d,md.mesh.x))
     565                        y3d=np.concatenate((y3d,md.mesh.y))
    562566                        #nodes are distributed between bed and surface accordingly to the given exponent
    563                         z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))
    564                 number_nodes3d=numpy.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
     567                        z3d=np.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))
     568                number_nodes3d=np.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
    565569
    566570                #Extrude elements
    567                 elements3d=numpy.empty((0,6),int)
     571                elements3d=np.empty((0,6),int)
    568572                for i in xrange(numlayers-1):
    569                         elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices))))    #Create the elements of the 3d mesh for the non extruded part
    570                 number_el3d=numpy.size(elements3d,axis=0)    #number of 3d nodes for the non extruded part of the mesh
     573                        elements3d=np.vstack((elements3d,np.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices))))    #Create the elements of the 3d mesh for the non extruded part
     574                number_el3d=np.size(elements3d,axis=0)    #number of 3d nodes for the non extruded part of the mesh
    571575
    572576                #Keep a trace of lower and upper nodes
    573                 lowervertex=-1*numpy.ones(number_nodes3d,int)
    574                 uppervertex=-1*numpy.ones(number_nodes3d,int)
    575                 lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
    576                 uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
     577                lowervertex=np.nan*np.ones(number_nodes3d,int)
     578                uppervertex=np.nan*np.ones(number_nodes3d,int)
     579                lowervertex[md.mesh.numberofvertices:]=np.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
     580                uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=np.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
    577581                md.mesh.lowervertex=lowervertex
    578582                md.mesh.uppervertex=uppervertex
    579583
    580584                #same for lower and upper elements
    581                 lowerelements=-1*numpy.ones(number_el3d,int)
    582                 upperelements=-1*numpy.ones(number_el3d,int)
    583                 lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
    584                 upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
     585                lowerelements=np.nan*np.ones(number_el3d,int)
     586                upperelements=np.nan*np.ones(number_el3d,int)
     587                lowerelements[md.mesh.numberofelements:]=np.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
     588                upperelements[:(numlayers-2)*md.mesh.numberofelements]=np.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
    585589                md.mesh.lowerelements=lowerelements
    586590                md.mesh.upperelements=upperelements
    … …  
    605609
    606610                #bedinfo and surface info
    607                 md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
    608                 md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
     611                md.mesh.vertexonbase=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
     612                md.mesh.vertexonsurface=project3d(md,'vector',np.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
    609613                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
    610614
    … …  
    630634
    631635                #connectivity
    632                 md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
    633                 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
    634                 if not numpy.isnan(md.mesh.elementconnectivity).all():
     636                md.mesh.elementconnectivity=np.tile(md.mesh.elementconnectivity,(numlayers-1,1))
     637                md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
     638                if not np.isnan(md.mesh.elementconnectivity).all():
    635639                        for i in xrange(1,numlayers-1):
    636640                                md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \
    637641                                                =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
    638                                 md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
     642                                md.mesh.elementconnectivity[np.nonzero(md.mesh.elementconnectivity<0)]=0
    639643
    640644                md.materials.extrude(md)
    … …  
    674678
    675679                #observations
    676                 if not numpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)
    677                 if not numpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)
    678                 if not numpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)
    679                 if not numpy.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers)
    680                 if isinstance(md.inversion.min_parameters,numpy.ndarray):
     680                if not np.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)
     681                if not np.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)
     682                if not np.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)
     683                if not np.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers)
     684                if isinstance(md.inversion.min_parameters,np.ndarray):
    681685                    if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers)
    682                 if isinstance(md.inversion.max_parameters,numpy.ndarray):
     686                if isinstance(md.inversion.max_parameters,np.ndarray):
    683687                    if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers)
    684                 if not numpy.isnan(md.smb.mass_balance).all():
     688                if not np.isnan(md.smb.mass_balance).all():
    685689                        md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers)
    686690               
    687691                #results
    688                 if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)
    689                 if not numpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy)
    690                 if not numpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz)
    691                 if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)
    692                 if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
    693                 if not numpy.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1)
    694                 if not numpy.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1)
    695                 if not numpy.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1)
    696                 if not numpy.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1)
    697 
    698                 #gia
    699                 if not numpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)
    700                 if not numpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)
     692                if not np.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)
     693                if not np.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy)
     694                if not np.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz)
     695                if not np.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)
     696                if not np.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
     697                if not np.isnan(md.initialization.pressure).all(): md.initialization.pressure=project2d(md,md.initialization.pressure,1)
     698                if not np.isnan(md.initialization.sediment_head).all(): md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1)
     699                if not np.isnan(md.initialization.epl_head).all(): md.initialization.epl_head=project2d(md,md.initialization.epl_head,1)
     700                if not np.isnan(md.initialization.epl_thickness).all(): md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1)
     701
     702                #giaivins
     703                if not np.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)
     704                if not np.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)
    701705
    702706                #elementstype
    703                 if not numpy.isnan(md.flowequation.element_equation).all():
     707                if not np.isnan(md.flowequation.element_equation).all():
    704708                        md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1)
    705709                        md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1)
    … …  
    725729                md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers)
    726730                md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers)
    727                 if not numpy.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)
     731                if not np.isnan(md.damage.spcdamage).all(): md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)
    728732                md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1)
    729733
    … …  
    752756                md.geometry.thickness=project2d(md,md.geometry.thickness,1)
    753757                md.geometry.base=project2d(md,md.geometry.base,1)
    754                 if isinstance(md.geometry.bed,numpy.ndarray):
     758                if isinstance(md.geometry.bed,np.ndarray):
    755759                    md.geometry.bed=project2d(md,md.geometry.bed,1)
    756760                md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1)
    757761                md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1)
    758 
    759                 #lat long
    760                 if isinstance(md.mesh.lat,numpy.ndarray):
    761                     if md.mesh.lat.size==md.mesh.numberofvertices:  md.mesh.lat=project2d(md,md.mesh.lat,1)
    762                 if isinstance(md.mesh.long,numpy.ndarray):
    763                     if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1)
    764762
    765763                #Initialize with the 2d mesh
    … …  
    770768                mesh.numberofelements=md.mesh.numberofelements2d
    771769                mesh.elements=md.mesh.elements2d
    772                 if not numpy.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)
    773                 if not numpy.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)
     770                if not np.isnan(md.mesh.vertexonboundary).all(): mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)
     771                if not np.isnan(md.mesh.elementconnectivity).all(): mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)
     772                if isinstance(md.mesh.lat,np.ndarray):
     773                        if md.mesh.lat.size==md.mesh.numberofvertices:  mesh.lat=project2d(md,md.mesh.lat,1)
     774                if isinstance(md.mesh.long,np.ndarray):
     775                        if md.mesh.long.size==md.mesh.numberofvertices: mesh.long=project2d(md,md.mesh.long,1)
     776                mesh.epsg=md.mesh.epsg
    774777                md.mesh=mesh
    775778
Note: See TracChangeset for help on using the changeset viewer.