Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 14016) +++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 14017) @@ -50,9 +50,9 @@ #Flag the elements that have not been flagged as filltype if strcmpi(filltype,'hutter'): - hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1 + hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True elif strcmpi(filltype,'macayeal'): - macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1 + macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True elif strcmpi(filltype,'pattyn'): - pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1 + pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True #check that each element has at least one flag @@ -63,7 +63,7 @@ if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1): print "setflowequation warning message: some elements have several types, higher order type is used for them" - hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=0 - hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0 - macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0 + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False + hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False + macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False #Check that no pattyn or stokes for 2d mesh @@ -77,14 +77,14 @@ #Initialize node fields - nodeonhutter=numpy.zeros(md.mesh.numberofvertices) - nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=1 - nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices) - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 - nodeonl1l2=numpy.zeros(md.mesh.numberofvertices) - nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=1 - nodeonpattyn=numpy.zeros(md.mesh.numberofvertices) - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - noneflag=numpy.zeros(md.mesh.numberofelements) + nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True + nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True + nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True + nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True + nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) + noneflag=numpy.zeros(md.mesh.numberofelements,bool) #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) @@ -97,26 +97,26 @@ # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront - stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 + stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True #Then complete with NoneApproximation or the other model used if there is no stokes if any(stokesflag): if any(pattynflag): #fill with pattyn - pattynflag[numpy.logical_not(stokesflag)]=1 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 + pattynflag[numpy.logical_not(stokesflag)]=True + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True elif any(macayealflag): #fill with macayeal - macayealflag[numpy.logical_not(stokesflag)]=1 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 + macayealflag[numpy.logical_not(stokesflag)]=True + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True else: #fill with none - noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1 + noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True #Now take care of the coupling between MacAyeal and Pattyn md.diagnostic.vertex_pairing=numpy.array([]) - nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices) - nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices) - macayealpattynflag=numpy.zeros(md.mesh.numberofelements) - macayealstokesflag=numpy.zeros(md.mesh.numberofelements) - pattynstokesflag=numpy.zeros(md.mesh.numberofelements) + nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool) + macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool) + macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool) + pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool) if strcmpi(coupling_method,'penalties'): #Create the border nodes between Pattyn and MacAyeal and extrude them @@ -135,13 +135,13 @@ if any(macayealflag) and any(pattynflag): #coupling macayeal pattyn #Find node at the border - nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1 + nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True #Macayeal elements in contact with this layer become MacAyealPattyn elements matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0]) commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal - macayealflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements - macayealpattynflag[numpy.nonzero(commonelements)]=1 - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 + commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal + macayealflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements + macayealpattynflag[numpy.nonzero(commonelements)]=True + nodeonmacayeal[:]=False + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True #rule out elements that don't touch the 2 boundaries @@ -151,29 +151,29 @@ elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] - macayealflag[pos[pos1]]=1 - macayealpattynflag[pos[pos1]]=0 + macayealflag[pos[pos1]]=True + macayealpattynflag[pos[pos1]]=False pos2=numpy.nonzero(elist==-1)[0] - pattynflag[pos[pos2]]=1 - macayealpattynflag[pos[pos2]]=0 + pattynflag[pos[pos2]]=True + macayealpattynflag[pos[pos2]]=False #Recompute nodes associated to these elements - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 - nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 - nodeonmacayealpattyn[:]=0 - nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=1 + nodeonmacayeal[:]=False + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True + nodeonpattyn[:]=False + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True + nodeonmacayealpattyn[:]=False + nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True elif any(pattynflag) and any(stokesflag): #coupling pattyn stokes #Find node at the border - nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1 + nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True #Stokes elements in contact with this layer become PattynStokes elements matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0]) commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(pattynflag)]=0 #only one layer: the elements previously in macayeal - stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealpattynelements - pattynstokesflag[numpy.nonzero(commonelements)]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 + commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements + pattynstokesflag[numpy.nonzero(commonelements)]=True + nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True #rule out elements that don't touch the 2 boundaries @@ -183,29 +183,29 @@ elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] - stokesflag[pos[pos1]]=1 - pattynstokesflag[pos[pos1]]=0 + stokesflag[pos[pos1]]=True + pattynstokesflag[pos[pos1]]=False pos2=numpy.nonzero(elist==-1)[0] - pattynflag[pos[pos2]]=1 - pattynstokesflag[pos[pos2]]=0 + pattynflag[pos[pos2]]=True + pattynstokesflag[pos[pos2]]=False #Recompute nodes associated to these elements - nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 - nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 - nodeonpattynstokes[:]=0 - nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=1 + nodeonstokes[:]=False + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True + nodeonpattyn[:]=False + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True + nodeonpattynstokes[:]=False + nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True elif any(stokesflag) and any(macayealflag): #Find node at the border - nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1 + nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True #Stokes elements in contact with this layer become MacAyealStokes elements matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0]) commonelements=numpy.sum(matrixelements,axis=1)!=0 - commonelements[numpy.nonzero(macayealflag)]=0 #only one layer: the elements previously in macayeal - stokesflag[numpy.nonzero(commonelements)]=0 #these elements are now macayealmacayealelements - macayealstokesflag[numpy.nonzero(commonelements)]=1 - nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 + commonelements[numpy.nonzero(macayealflag)]=False #only one layer: the elements previously in macayeal + stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealmacayealelements + macayealstokesflag[numpy.nonzero(commonelements)]=True + nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True #rule out elements that don't touch the 2 boundaries @@ -215,17 +215,17 @@ elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] - macayealflag[pos[pos1]]=1 - macayealstokesflag[pos[pos1]]=0 + macayealflag[pos[pos1]]=True + macayealstokesflag[pos[pos1]]=False pos2=numpy.nonzero(elist==-1)[0] - stokesflag[pos[pos2]]=1 - macayealstokesflag[pos[pos2]]=0 + stokesflag[pos[pos2]]=True + macayealstokesflag[pos[pos2]]=False #Recompute nodes associated to these elements - nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 - nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 - nodeonmacayealstokes[:]=0 - nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=1 + nodeonmacayeal[:]=False + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True + nodeonstokes[:]=False + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True + nodeonmacayealstokes[:]=False + nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True elif any(stokesflag) and any(hutterflag):