Changeset 26744 for issm/trunk/src/m/plot/processmesh.py

Timestamp:

12/22/21 10:39:44 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 26742

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/plot/processmesh.py (modified) (2 diffs)

issm/trunk/src/m/plot/processmesh.py

@@ -41 +41 @@
     if hasattr(md.mesh, 'elements2d'):
         elements2d = md.mesh.elements2d
+        numofvertices2d = md.mesh.numberofvertices2d
+        numofelements2d = md.mesh.numberofelements2d
+    else:
+        numofvertices2d = np.nan
+        numofelements2d = np.nan
 
     if options.exist('amr'):
         step = options.getfieldvalue('amr')
-        x = md.results.TransientSolution[step].MeshX
-        y = md.results.TransientSolution[step].MeshY
-        elements = md.results.TransientSolution[step].MeshElements
+        nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshX))
+        x = md.results.TransientSolution[step].MeshX[nonan]
+        y = md.results.TransientSolution[step].MeshY[nonan]
+        nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshElements))
+        elements = md.results.TransientSolution[step].MeshElements[nonan] - 1
+        eldim = np.shape(md.results.TransientSolution[step].MeshElements)[1]
+        elements = np.reshape(elements, ((int(len(elements) / eldim), eldim)))
+
     else:
-        elements = md.mesh.elements
+        elements = md.mesh.elements - 1
         if options.getfieldvalue('coord', 'xy') != 'latlon':
             x = md.mesh.x
@@ -65 +75 @@
         z = getattr(md, z)
 
+    force2D = numofelements2d in np.shape(data) or numofvertices2d in np.shape(data)
     #is it a 2D plot?
-    if md.mesh.dimension() == 2 or options.getfieldvalue('layer', 0) >= 1:
+    if md.mesh.dimension() == 2 or options.getfieldvalue('layer', 0) >= 1 or force2D:
         is2d = 1
     else:
         is2d = 0
 
-    elements = md.mesh.elements - 1
-
     #layer projection?
-    if options.getfieldvalue('layer', 0) >= 1:
+    if options.getfieldvalue('layer', 0) >= 1 or force2D:
         if options.getfieldvalue('coord', 'xy') == 'latlon':
             raise Exception('processmesh error message: cannot work with 3D meshes for now')