Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24274) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24275) @@ -4,5 +4,5 @@ -def exportVTK(filename, md, *args, enveloppe=False): +def exportVTK(filename, md, *args, enveloppe=False, **kwargs): ''' vtk export @@ -20,8 +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,23 +125,83 @@ 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,5 +214,5 @@ 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,2 +416,10 @@ 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."""