Changeset 24275

Timestamp:

10/25/19 02:44:52 (5 years ago)

Author:

bdef

Message:

NEW: adding some user required options to eportvtk

File:

• issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (modified) (5 diffs)

issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py

@@ -4 +4 @@
 
 
-def exportVTK(filename, md, *args, enveloppe=False):
+def exportVTK(filename, md, *args, enveloppe=False, **kwargs):
     '''
     vtk export
@@ -20 +20 @@
     enveloppe is an option keeping only the enveloppe of the md (it is False by default)
 
+    Options:
+        - clipping : allows to reduce your domain (cliping=[Xmin, Xmax, Ymin, Ymax])
+        - coarsetime : output one timestep every X (coarsetime=X, with X an integer)
+        - singletime : output only timestep X (singletime=X, with X an integer or -1 for last)
+
     TODO: - make time easily accessible
 
     Basile de Fleurian:
     '''
+
+    for key in kwargs.keys():
+        if key not in ['clipping', 'coarsetime', 'singletime']:
+            raise BadOption('Provided option "{}" is not supported possibilities are : {}'.format(key, ['cliping', 'coarsetime']))
+
+    if 'coarsetime' in kwargs.keys() and 'singletime' in kwargs.keys():
+        raise BadOption("You can't specify both 'coarsetime' and 'singletime'")
 
     # File checking and creation {{{
@@ -113 +125 @@
                                         3 : md.geometry.bed
                                         4 : 0\n''')
-            if mesh_alti == 1:
+            if mesh_alti == '1':
                 points = np.column_stack((md.mesh.x, md.mesh.y, md.geometry.surface))
-            elif mesh_alti == 2:
+            elif mesh_alti == '2':
                 points = np.column_stack((md.mesh.x, md.mesh.y, md.geometry.base))
-            elif mesh_alti == 3:
+            elif mesh_alti == '3':
                 points = np.column_stack((md.mesh.x, md.mesh.y, md.geometry.bed))
-            elif mesh_alti == 4:
+            elif mesh_alti == '4':
                 points = np.column_stack((md.mesh.x, md.mesh.y, 0. * md.mesh.x))
             else:
                 points = np.column_stack((md.mesh.x, md.mesh.y, md.geometry.surface))
         elif dim == 3:
-            mesh_alti = 1
+            mesh_alti = '1'
             points = np.column_stack((md.mesh.x, md.mesh.y, md.mesh.z))
         else:
             raise BadDimension('exportVTK does not support dimension {}'.format(dim))
+
+    if 'clipping' in kwargs.keys():
+        # first get the boundaries and check them
+        [Xmin, Xmax, Ymin, Ymax] = kwargs['clipping']
+        if Xmin > Xmax:
+            raise ClipError('Xmax ({}) should be larger than Xmin ({})'.format(Xmax, Xmin))
+        if Ymin > Ymax:
+            raise ClipError('Ymax ({}) should be larger than Ymin ({})'.format(Ymax, Ymin))
+        if Xmin > np.nanmax(points[:, 0]) or Xmax < np.nanmin(points[:, 0]):
+            raise ClipError('Your X boundaries [{}, {}] are outside of the model domain [{},{}]'.format(Xmin, Xmax, np.nanmin(points[:, 0]), np.nanmax(points[:, 0])))
+        if Ymin > np.nanmax(points[:, 1]) or Ymax < np.nanmin(points[:, 1]):
+            raise ClipError('Your Y boundaries [{}, {}] are outside of the model domain [{},{}]'.format(Ymin, Ymax, np.nanmin(points[:, 1]), np.nanmax(points[:, 1])))
+
+        #boundaries should be fine lets do stuff
+        InX = np.where(np.logical_and(points[:, 0] >= Xmin, points[:, 0] <= Xmax))
+        InY = np.where(np.logical_and(points[:, 1] >= Ymin, points[:, 1] <= Ymax))
+
+        Isinside = np.zeros(np.shape(points)[0], dtype=bool)
+        clip_convert_index = np.nan * np.ones(np.shape(points)[0])
+
+        Inclipping = np.intersect1d(InX, InY)
+        Isinside[Inclipping] = True
+        points = points[Inclipping, :]
+        num_of_points = np.shape(points)[0]
+
+        clipconnect = np.asarray([], dtype=int)
+        for elt in connect:
+            if set(elt).issubset(Inclipping):
+                clipconnect = np.append(clipconnect, elt, axis=0)
+
+        num_of_elt = int(np.size(clipconnect) / 3)
+        connect = clipconnect.reshape(num_of_elt, 3)
+
+        clip_convert_index = np.asarray([[i, np.where(Inclipping == i)[0][0]] for i, val in enumerate(clip_convert_index) if any(Inclipping == i)])
+        enveloppe_index = enveloppe_index[clip_convert_index[:, 0]]
+        for elt in range(0, num_of_elt):
+            try:
+                connect[elt, 0] = clip_convert_index[np.where(clip_convert_index == connect[elt, 0])[0], 1][0]
+            except IndexError:
+                connect[elt, 0] = -1
+            try:
+                connect[elt, 1] = clip_convert_index[np.where(clip_convert_index == connect[elt, 1])[0], 1][0]
+            except IndexError:
+                connect[elt, 1] = -1
+            try:
+                connect[elt, 2] = clip_convert_index[np.where(clip_convert_index == connect[elt, 2])[0], 1][0]
+            except IndexError:
+                connect[elt, 2] = -1
+
+        connect = connect[np.where(connect != -1)[0], :]
+        num_of_elt = np.shape(connect)[0]
+
     # }}}
     # write header and mesh {{{
     print('Now starting to write stuff')
-    for step in range(0, num_of_timesteps):
+
+    if 'coarsetime' in kwargs.keys():
+        steplist = range(0, num_of_timesteps, kwargs['coarsetime'])
+    elif 'singletime' in kwargs.keys():
+        steplist = [kwargs['singletime']]
+    else:
+        steplist = range(0, num_of_timesteps)
+
+    for step in steplist:
         print('Writing for step {}'.format(step))
         saved_cells = {}
@@ -142 +214 @@
         fid.write('POINTS {:d} float\n'.format(num_of_points))
     #updating z for mesh evolution
-        if moving_mesh and mesh_alti in [1, 2]:
+        if moving_mesh and mesh_alti in ['1', '2']:
             base = np.squeeze(res_struct.__dict__['TransientSolution'][step].__dict__['Base'][enveloppe_index])
             thick_change_ratio = (np.squeeze(res_struct.__dict__['TransientSolution'][step].__dict__['Thickness'][enveloppe_index]) / md.geometry.thickness[enveloppe_index])
@@ -344 +416 @@
 class BadDimension(Exception):
     """The required dimension is not supported yet."""
+
+
+class BadOption(Exception):
+    """The given option does not exist."""
+
+
+class ClipError(Exception):
+    """Error while trying to clip the domain."""