Changeset 14017
- Timestamp:
- 11/27/12 14:32:53 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/parameterization/setflowequation.py
r14015 r14017 50 50 #Flag the elements that have not been flagged as filltype 51 51 if strcmpi(filltype,'hutter'): 52 hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]= 152 hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True 53 53 elif strcmpi(filltype,'macayeal'): 54 macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]= 154 macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True 55 55 elif strcmpi(filltype,'pattyn'): 56 pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]= 156 pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True 57 57 58 58 #check that each element has at least one flag … … 63 63 if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1): 64 64 print "setflowequation warning message: some elements have several types, higher order type is used for them" 65 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]= 066 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]= 067 macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]= 065 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False 66 hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False 67 macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False 68 68 69 69 #Check that no pattyn or stokes for 2d mesh … … 77 77 78 78 #Initialize node fields 79 nodeonhutter=numpy.zeros(md.mesh.numberofvertices )80 nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]= 181 nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices )82 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]= 183 nodeonl1l2=numpy.zeros(md.mesh.numberofvertices )84 nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]= 185 nodeonpattyn=numpy.zeros(md.mesh.numberofvertices )86 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]= 187 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )88 noneflag=numpy.zeros(md.mesh.numberofelements )79 nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool) 80 nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True 81 nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool) 82 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 83 nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool) 84 nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True 85 nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool) 86 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 87 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 88 noneflag=numpy.zeros(md.mesh.numberofelements,bool) 89 89 90 90 #First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) … … 97 97 # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront 98 98 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 99 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]= 0100 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]= 199 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False 100 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 101 101 102 102 #Then complete with NoneApproximation or the other model used if there is no stokes 103 103 if any(stokesflag): 104 104 if any(pattynflag): #fill with pattyn 105 pattynflag[numpy.logical_not(stokesflag)]= 1106 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]= 1105 pattynflag[numpy.logical_not(stokesflag)]=True 106 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 107 107 elif any(macayealflag): #fill with macayeal 108 macayealflag[numpy.logical_not(stokesflag)]= 1109 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]= 1108 macayealflag[numpy.logical_not(stokesflag)]=True 109 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 110 110 else: #fill with none 111 noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]= 1111 noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True 112 112 113 113 #Now take care of the coupling between MacAyeal and Pattyn 114 114 md.diagnostic.vertex_pairing=numpy.array([]) 115 nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices )116 nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices )117 nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices )118 macayealpattynflag=numpy.zeros(md.mesh.numberofelements )119 macayealstokesflag=numpy.zeros(md.mesh.numberofelements )120 pattynstokesflag=numpy.zeros(md.mesh.numberofelements )115 nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool) 116 nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool) 117 nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool) 118 macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool) 119 macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool) 120 pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool) 121 121 if strcmpi(coupling_method,'penalties'): 122 122 #Create the border nodes between Pattyn and MacAyeal and extrude them … … 135 135 if any(macayealflag) and any(pattynflag): #coupling macayeal pattyn 136 136 #Find node at the border 137 nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]= 1137 nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True 138 138 #Macayeal elements in contact with this layer become MacAyealPattyn elements 139 139 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0]) 140 140 commonelements=numpy.sum(matrixelements,axis=1)!=0 141 commonelements[numpy.nonzero(pattynflag)]= 0#only one layer: the elements previously in macayeal142 macayealflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealpattynelements143 macayealpattynflag[numpy.nonzero(commonelements)]= 1144 nodeonmacayeal[:]= 0145 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]= 1141 commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal 142 macayealflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements 143 macayealpattynflag[numpy.nonzero(commonelements)]=True 144 nodeonmacayeal[:]=False 145 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 146 146 147 147 #rule out elements that don't touch the 2 boundaries … … 151 151 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) 152 152 pos1=numpy.nonzero(elist==1)[0] 153 macayealflag[pos[pos1]]= 1154 macayealpattynflag[pos[pos1]]= 0153 macayealflag[pos[pos1]]=True 154 macayealpattynflag[pos[pos1]]=False 155 155 pos2=numpy.nonzero(elist==-1)[0] 156 pattynflag[pos[pos2]]= 1157 macayealpattynflag[pos[pos2]]= 0156 pattynflag[pos[pos2]]=True 157 macayealpattynflag[pos[pos2]]=False 158 158 159 159 #Recompute nodes associated to these elements 160 nodeonmacayeal[:]= 0161 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]= 1162 nodeonpattyn[:]= 0163 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]= 1164 nodeonmacayealpattyn[:]= 0165 nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]= 1160 nodeonmacayeal[:]=False 161 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 162 nodeonpattyn[:]=False 163 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 164 nodeonmacayealpattyn[:]=False 165 nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True 166 166 167 167 elif any(pattynflag) and any(stokesflag): #coupling pattyn stokes 168 168 #Find node at the border 169 nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]= 1169 nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True 170 170 #Stokes elements in contact with this layer become PattynStokes elements 171 171 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0]) 172 172 commonelements=numpy.sum(matrixelements,axis=1)!=0 173 commonelements[numpy.nonzero(pattynflag)]= 0#only one layer: the elements previously in macayeal174 stokesflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealpattynelements175 pattynstokesflag[numpy.nonzero(commonelements)]= 1176 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )177 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]= 1173 commonelements[numpy.nonzero(pattynflag)]=False #only one layer: the elements previously in macayeal 174 stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealpattynelements 175 pattynstokesflag[numpy.nonzero(commonelements)]=True 176 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 177 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 178 178 179 179 #rule out elements that don't touch the 2 boundaries … … 183 183 elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool) 184 184 pos1=numpy.nonzero(elist==1)[0] 185 stokesflag[pos[pos1]]= 1186 pattynstokesflag[pos[pos1]]= 0185 stokesflag[pos[pos1]]=True 186 pattynstokesflag[pos[pos1]]=False 187 187 pos2=numpy.nonzero(elist==-1)[0] 188 pattynflag[pos[pos2]]= 1189 pattynstokesflag[pos[pos2]]= 0188 pattynflag[pos[pos2]]=True 189 pattynstokesflag[pos[pos2]]=False 190 190 191 191 #Recompute nodes associated to these elements 192 nodeonstokes[:]= 0193 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]= 1194 nodeonpattyn[:]= 0195 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]= 1196 nodeonpattynstokes[:]= 0197 nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]= 1192 nodeonstokes[:]=False 193 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 194 nodeonpattyn[:]=False 195 nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True 196 nodeonpattynstokes[:]=False 197 nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True 198 198 199 199 elif any(stokesflag) and any(macayealflag): 200 200 #Find node at the border 201 nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]= 1201 nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True 202 202 #Stokes elements in contact with this layer become MacAyealStokes elements 203 203 matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0]) 204 204 commonelements=numpy.sum(matrixelements,axis=1)!=0 205 commonelements[numpy.nonzero(macayealflag)]= 0#only one layer: the elements previously in macayeal206 stokesflag[numpy.nonzero(commonelements)]= 0#these elements are now macayealmacayealelements207 macayealstokesflag[numpy.nonzero(commonelements)]= 1208 nodeonstokes=numpy.zeros(md.mesh.numberofvertices )209 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]= 1205 commonelements[numpy.nonzero(macayealflag)]=False #only one layer: the elements previously in macayeal 206 stokesflag[numpy.nonzero(commonelements)]=False #these elements are now macayealmacayealelements 207 macayealstokesflag[numpy.nonzero(commonelements)]=True 208 nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool) 209 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 210 210 211 211 #rule out elements that don't touch the 2 boundaries … … 215 215 elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) 216 216 pos1=numpy.nonzero(elist==1)[0] 217 macayealflag[pos[pos1]]= 1218 macayealstokesflag[pos[pos1]]= 0217 macayealflag[pos[pos1]]=True 218 macayealstokesflag[pos[pos1]]=False 219 219 pos2=numpy.nonzero(elist==-1)[0] 220 stokesflag[pos[pos2]]= 1221 macayealstokesflag[pos[pos2]]= 0220 stokesflag[pos[pos2]]=True 221 macayealstokesflag[pos[pos2]]=False 222 222 223 223 #Recompute nodes associated to these elements 224 nodeonmacayeal[:]= 0225 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]= 1226 nodeonstokes[:]= 0227 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]= 1228 nodeonmacayealstokes[:]= 0229 nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]= 1224 nodeonmacayeal[:]=False 225 nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True 226 nodeonstokes[:]=False 227 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True 228 nodeonmacayealstokes[:]=False 229 nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True 230 230 231 231 elif any(stokesflag) and any(hutterflag):
Note:
See TracChangeset
for help on using the changeset viewer.