source: issm/oecreview/Archive/16554-17801/ISSM-17208-17209.diff

../trunk-jpl/test/NightlyRun/test336.py

@@ -7 +7 @@
 from solve import *
 from MatlabFuncs import *
 
-md=triangle(model(),'../Exp/Square.exp',150000.)
+md=triangle(model(),'../Exp/Square.exp',100000.)
 md=setmask(md,'','')
 md=parameterize(md,'../Par/SquareSheetConstrained.py') 
 md=setflowequation(md,'SSA','all')
@@ -19 +19 @@
 md.transient.isthermal=False
 md.transient.isgroundingline=False
 md.transient.isgia=False
-md.transient.islevelset=True
+md.transient.islevelset=False
 
 # init levelset function
 ymin=numpy.min(md.mesh.y)
 ymax=numpy.max(md.mesh.y)
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-xmed=(xmax+xmin)/2
-ymed=(ymax+ymin)/2
-distx=numpy.absolute(md.mesh.x.reshape(-1,1)-xmed)
-disty=numpy.absolute(md.mesh.y.reshape(-1,1)-ymed)
-md.mask.ice_levelset=numpy.maximum(distx,disty)-1.e5
+xmed=(xmax+xmin)/2.
+ymed=(ymax+ymin)/2.
+md.mask.ice_levelset=numpy.sqrt(numpy.power(md.mesh.x-xmed,2.)+numpy.power(md.mesh.y-ymed,2.)) - (xmax-xmin)/3.
 
+# set spcs
+mask=1.*numpy.ones((md.mesh.numberofvertices,1))
+nrverts=md.mesh.elements.shape[1]
+for i in range(0,md.mesh.numberofelements):
+        elt=numpy.copy(md.mesh.elements[i])
+        elt-=1
+        isiceinelement=False
+        for iv in range(0,nrverts):
+                if(md.mask.ice_levelset[elt[iv]]<=0.):
+                        isiceinelement=True
+        if(isiceinelement):
+                for iv in range(0,nrverts):
+                        mask[elt[iv]]=2.
+
+v=0.
+for i in range(0,md.mesh.numberofvertices):
+        if(mask[i]==1.):
+                md.stressbalance.spcvx[i]=v
+                md.stressbalance.spcvy[i]=v
+                md.stressbalance.spcvz[i]=v
+
 md=solve(md,TransientSolutionEnum())
 
 #Fields and tolerances to track changes
 field_names     =['Vx','Vy','Vel','Pressure','MaskIceLevelset']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]#,1e-13,1e-13]
 field_values=[\
         md.results.TransientSolution[0].Vx,\
         md.results.TransientSolution[0].Vy,\
         md.results.TransientSolution[0].Vel,\
         md.results.TransientSolution[0].Pressure,\
-        md.results.TransientSolution[0].MaskIceLevelset,\
+#       md.results.TransientSolution[0].MaskIceLevelset,\
         md.results.TransientSolution[1].Vx,\
         md.results.TransientSolution[1].Vy,\
         md.results.TransientSolution[1].Vel,\
         md.results.TransientSolution[1].Pressure,\
-        md.results.TransientSolution[1].MaskIceLevelset,\
+#       md.results.TransientSolution[1].MaskIceLevelset,\
         ]