Changeset 13710

Timestamp:

10/17/12 10:28:51 (12 years ago)

Author:

jschierm

Message:

FIX: Fixes to run python rifts (plus a matlab fix as well).

Location:

issm/trunk-jpl/src/m

Files:

• classes/rifts.m (modified) (1 diff)
• classes/rifts.py (modified) (2 diffs)
• mesh/rifts/meshprocessrifts.py (modified) (4 diffs)

issm/trunk-jpl/src/m/classes/rifts.m

@@ -38 +38 @@
                                         md = checkmessage(md,['model should be processed for rifts (run meshprocessrifts)!']);
                                 end
-                                md = checkfield(md,'rifts.riftstruct.fill','values',[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
+                                for i=1:numrifts,
+                                        md = checkfield(md,sprintf('rifts.riftstruct(%d).fill',i),'values',[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
+                                end
                         else
                                 if ~isnans(obj.riftstruct),

issm/trunk-jpl/src/m/classes/rifts.py

@@ -53 +53 @@
                                 #We have segments with rift markers, but no rift structure!
                                 md.checkmessage("model should be processed for rifts (run meshprocessrifts)!")
-                        md = checkfield(md,'rifts.riftstruct.fill','values',[WaterEnum(),AirEnum(),IceEnum(),MelangeEnum()])
+                        for i,rift in enumerate(self.riftstruct):
+                                md = checkfield(md,"rifts.riftstruct[%d]['fill']" % i,'values',[WaterEnum(),AirEnum(),IceEnum(),MelangeEnum()])
                 else:
                         if numpy.any(numpy.logical_not(isnans(self.riftstruct))):
-                                md.checkmessage("riftstruct shoud be NaN since numrifts is 0!")
+                                md.checkmessage("riftstruct should be NaN since numrifts is 0!")
 
                 return md
@@ -70 +71 @@
 
                 numpairs=0
-                for i in xrange(0,numrifts):
-                        numpairs+=numpy.size(self.riftstruct[i].penaltypairs,0)
+                for rift in self.riftstruct:
+                        numpairs+=numpy.size(rift['penaltypairs'],axis=0)
 
                 # 2 for nodes + 2 for elements+ 2 for  normals + 1 for length + 1 for fill + 1 for friction + 1 for fraction + 1 for fractionincrement + 1 for state.
                 data=numpy.zeros((numpairs,12))
                 count=0
-                for i in xrange(0,numrifts):
-                        numpairsforthisrift=numpy.size(self.riftstruct[i]['penaltypairs'],0)
-                        data[count:count+numpairsforthisrift-1,0:6]=self.riftstruct[i]['penaltypairs']
-                        data[count:count+numpairsforthisrift-1,7]=self.riftstruct[i]['fill']
-                        data[count:count+numpairsforthisrift-1,8]=self.riftstruct[i]['friction']
-                        data[count:count+numpairsforthisrift-1,9]=self.riftstruct[i]['fraction']
-                        data[count:count+numpairsforthisrift-1,10]=self.riftstruct[i]['fractionincrement']
-                        data[count:count+numpairsforthisrift-1,11]=self.riftstruct[i]['state']
+                for rift in self.riftstruct:
+                        numpairsforthisrift=numpy.size(rift['penaltypairs'],0)
+                        data[count:count+numpairsforthisrift,0:7]=rift['penaltypairs']
+                        data[count:count+numpairsforthisrift,7]=rift['fill']
+                        data[count:count+numpairsforthisrift,8]=rift['friction']
+                        data[count:count+numpairsforthisrift,9]=rift['fraction']
+                        data[count:count+numpairsforthisrift,10]=rift['fractionincrement']
+                        data[count:count+numpairsforthisrift,11]=rift['state'].reshape(-1)
                         count+=numpairsforthisrift
 

issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py

@@ -1 +1 @@
 import numpy
-#from TriMeshProcessRifts import *
+from TriMeshProcessRifts import *
 from ContourToMesh import *
 from meshprocessoutsiderifts import *
@@ -23 +23 @@
         #Call MEX file
         [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers)
-        if not isinstance(md.rifts.riftstruct,'list') or not md.rifts.riftstruct:
+        md.mesh.x=md.mesh.x.reshape(-1)
+        md.mesh.y=md.mesh.y.reshape(-1)
+        if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
                 raise RuntimeError("TriMeshProcessRifts did not find any rift")
 
@@ -40 +42 @@
         #get coordinates of rift tips
         for rift in md.rifts.riftstruct:
-                rift.tip1coordinates=numpy.hstack((md.mesh.x[rift.tips[0].astype(int)-1].reshape(-1,1),md.mesh.y[rift.tips[0].astype(int)-1].reshape(-1,1)))
-                rift.tip2coordinates=numpy.hstack((md.mesh.x[rift.tips[1].astype(int)-1].reshape(-1,1),md.mesh.y[rift.tips[1].astype(int)-1].reshape(-1,1)))
+                rift['tip1coordinates']=numpy.hstack((md.mesh.x[rift['tips'][0,0].astype(int)-1].reshape(-1,1),md.mesh.y[rift['tips'][0,0].astype(int)-1].reshape(-1,1)))
+                rift['tip2coordinates']=numpy.hstack((md.mesh.x[rift['tips'][0,1].astype(int)-1].reshape(-1,1),md.mesh.y[rift['tips'][0,1].astype(int)-1].reshape(-1,1)))
 
         #In case we have rifts that open up the domain outline, we need to open them: 
-        flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
+        [flags,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
         found=0
         for rift in md.rifts.riftstruct:
-                if flags[rift.tips[0].astype(int)-1]==0:
+                if flags[rift['tips'][0,0].astype(int)-1]==0:
                         found=1
                         break
-                if flags[rift.tips[1].astype(int)-1]==0:
+                if flags[rift['tips'][0,1].astype(int)-1]==0:
                         found=1
                         break
@@ -59 +61 @@
         aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y)
         pos=numpy.nonzero(aires<0)[0]
-        md.mesh.elements[pos,:]=numpy.array([[md.mesh.elements[pos,1],md.mesh.elements[pos,0],md.mesh.elements[pos,2]]])
+        md.mesh.elements[pos,:]=numpy.hstack((md.mesh.elements[pos,1].reshape(-1,1),md.mesh.elements[pos,0].reshape(-1,1),md.mesh.elements[pos,2].reshape(-1,1)))
 
         return md