slimPre.py 41.5 KB
Newer Older
1
2
3
4
5
6
7
"""Second-generation Louvain-la-Neuve Ice-ocean model - User interface

    Contact: jonathan.lambrechts@uclouvain.be
    webpage: www.climate.be/slim

    SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.

8
    slimPre contains the tools writing the input files necessary for the simulations (for simulations, see slim.py file), such as partitioned mesh, bathymetry, coriolis, wind, tides, ...
9
10
11
12
    
    To generate the input files, run the script my_prepro_script.py:
        - On a single core: rundgpy my_prepro_script.py
        - On multiple cores: mpirun -np N rundgpy my_prepro_script.py (replace N by the number of cores)
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

    Mesh:
        A Mesh object contains a mesh, loaded from a .msh file, and the mesh projection system

    Region:
        A Region object contains a Mesh object, and a list of physical tags. The coordinates of all the mesh nodes contained on the physical tags is loaded in variable coordinates

    Coordinate_system:
        A Coordinate_system object contains a Region object, and the projection system of some data needed by the user. The variable coordinates gives the node coordinates of the Region in the data reference system.
        Vector field data can be rotate using the "rotate" and "rotate3D" functions

    Time:
        A Time object contains an array of times at which the data will be given. Time can be built using either a vector of instants relative to some initial time, or an initial and a final time and a time step. Some periodic window can also be built.

    write_file:
28
        The write_file function writes a netcdf file, using a Region object (giving the nodes for which there is data), a time object (giving the times of the data), and a list of data tuples.
29
30
        The first variable of a data tuple is the data name.
        The secod variable of a data tuple is the data array, with as many lines as the number of nodes, and as many columns as the number of time steps
31
32
33
    
    netcdf_to_msh:
        The function transforms a netcdf file to a .msh file.
34
35
36
37
38
39
40
41

    tpxo_tide:
        The function returns the data array of sea surface elevation, and the two horizontal transports (rotated in the coordinate system of the mesh, in m^2 s^-1) of the tides, using TPXO 7.2

    residual_flow:
        The function returns the data array of the two horizontal transports (in m^2 s^-1) at the nodes of some physical tags, such that the velocity is constant at all the nodes and that the total transport (in m^3 s^-1) is equal to the one prescribed by the user.

    smooth_bathymetry:
42
        The function smooths the given bathymetry, following a smoothing criterion, and writes it to a file. It can also write a bathymetry constant per element (P0), for a better stability. If the bathymetry is P0, it will be written in a idx file. Otherwise, it will be written in a netcdf file.
43
44
45
46
47
48
49
    
    interpolate_from_structured_grid
        The function interpolates from structured data 
    
    interpolate_from_structured_file
        The function interpolates structured data from a GMSH-file
    
50
51
52
53
    get_data_from_netcdf:
        The function returns the data array contained in the netcdf file given as an argument

    get_distance:
54
        The function returns an array with the minimum distance from each node of a given Region object to the physical tags given by the user. 
55
56

    fetch_ftp:
57
        This function downloads some data from the SLIM servers if available. The data is downloaded in the user global_data_directory. Data is only downloaded if it does not exist locally or is older than the server data. global_data_directory is by default ~/.cache/slim_data/
58
59
60
61
62

    set_global_data_directory:
        This functions gives the possibility to edit the global_data_directory

    make_directory:
63
        This function makes a directory if it does not exist
64
        
65
66
67
    partition_mesh:
        This function partitions the mesh (if run in multiprocessing) 
    
68
69
70
71
72
73
    partition_id:
        This function returns the number of the current partition as a string
    
    partition_nb:
        This function returns the number of partitions used as a string
    
74
    exit:
75
        This function closes the program
76
"""
Philippe Delandmeter's avatar
Philippe Delandmeter committed
77
import dgpy
78
from dgpy.scripts import slim_private
Philippe Delandmeter's avatar
Philippe Delandmeter committed
79
80
81
import numpy as np

class Mesh :
82
    """Create the slimPre Mesh"""
83
    def __init__(self, mesh_file_name, mesh_proj=None) :
84
85
86
        """keyword arguments:
        
        * mesh_file_name            
87
            path to the mesh file (.msh format)
88
        * mesh_proj           
89
            mesh projection system (Default: None). 
90
91
            
            This argument must not be given if the mesh is written in 3D (on a sphere), or if all the data are in the same projection system than the mesh.
92
        """
93
94
95
        isPart = slim_private._check_if_partionned(mesh_file_name)
        if (dgpy.Msg.GetCommSize() > 1) and not isPart:
            dgpy.Msg.Fatal("Calling slimPre.Mesh in multiprocessing requires providing a partitioned mesh (first use slimPre.partition_mesh).")
96
97
#        elif (dgpy.Msg.GetCommSize() == 1) and isPart:
#            dgpy.Msg.Fatal("Calling slimPre.Mesh on a single core requires providing a non-partitioned mesh.")
98
        self._groups = dgpy.dgGroupCollection(mesh_file_name)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
99
        self._mesh_proj = mesh_proj
100
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
101

Philippe Delandmeter's avatar
Philippe Delandmeter committed
102
class Region :
103
104
105
106
    """Create the slimPre Region

        region.coordinates is an array of shape (n,3), n being the number of nodes in the region, with the coordinates of all the region nodes
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
107
    def __init__(self, mesh, physical_tags='') :
108
109
110
        """keyword arguments:
        
        * mesh
111
            the slimPre Mesh object
112
        * physical_tags
113
114
115
            the list of all the tags on which put some data (optional)
            
            If physical_tags are not set, all the mesh nodes are gathered
116
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
117
        self._mesh = mesh
118
119
        if slim_private._is_string(physical_tags):
            physical_tags = [physical_tags]
Philippe Delandmeter's avatar
Philippe Delandmeter committed
120
121
122
123
124
        self._physical_tags = physical_tags
        self._nodes = []
        self._groups = []
        self._node_map = {}
        self._node_ids = []
125
        self._empty = False
Philippe Delandmeter's avatar
Philippe Delandmeter committed
126
127
        num_ = 0
        self._type = "interface"
128
        for iGroup in range(mesh._groups.getNbElementGroups()+mesh._groups.getNbGhostGroups()) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
129
            group = mesh._groups.getElementGroup(iGroup)
130
            if ((self._physical_tags != ['']) and (not group.getPhysicalTag() in self._physical_tags)) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
                continue
            self._type = "volume"
            self._groups.append(group)
            nbNodes = group.getNbNodes()
            ig = len(self._groups) - 1
            for iElement in range(group.getNbElements()) :
                element = group.getElement(iElement)
                for iNode in range(nbNodes) :
                    nodeId = mesh._groups.mesh().gmshVertexId(element.vertex(iNode))
                    if not nodeId in self._node_map :
                        self._nodes.append((ig, nbNodes*iElement+iNode))
                        num_ += 1
                        self._node_map[nodeId] = num_
                        self._node_ids.append(nodeId)
        if self._type == "volume" :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
146
            self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
147
148
149
150
            return

        for iFaceGroup in range(mesh._groups.getNbFaceGroups()) :
            faceGroup = mesh._groups.getFaceGroup(iFaceGroup)
151
            if (not faceGroup.physicalTag() in self._physical_tags) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
                continue
            self._groups.append(faceGroup)
            nbNodes = faceGroup.getNbNodes()
            ig = len(self._groups) - 1
            for iFace in range(faceGroup.size()) :
                group = faceGroup.elementGroup(0)
                iElement = faceGroup.elementId(iFace, 0)
                element = group.getElement(iElement)
                cl = faceGroup.closure(iFace, 0)
                for iNode in range(nbNodes) :
                    nodeId = mesh._groups.mesh().gmshVertexId(element.vertex(cl[iNode]))
                    if not nodeId in self._node_map :
                        self._nodes.append((ig, nbNodes*iFace+iNode))
                        num_ += 1
                        self._node_map[nodeId] = num_
                        self._node_ids.append(nodeId)
        self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
169
170
        if len(self._nodes) == 0:
            self._empty = True
Philippe Delandmeter's avatar
Philippe Delandmeter committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

    def _evaluateFunctor(self, f, size) :
        fc = dgpy.functorCache(dgpy.functorCache.NODE_GROUP_MODE, self._mesh._groups)
        val = np.empty((len(self._nodes), size))
        valg = []
        for g in self._groups :
            if self._type == "volume":
                fc.setGroup(g)
            else:
                fc.setInterfaceGroup(g)
            valg.append(np.copy(fc.get(f)))
        for i,(iG, iN) in enumerate(self._nodes) :
            val[i,:] = valg[iG][iN]
        return val


class Coordinate_system :
188
189
190
191
    """Create the slimPre Coordinate_system

        Coordinate_system.coordinates is an array of shape (n,3), n being the number of nodes in the region, with the coordinates of all the region nodes in the data projection system.
        Coordinate_system.rotate(u,v) returns the rotated vector array (u,v)
192
        Coordinate_system.rotate3D ... 
193
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
194
    def __init__(self, region, data_proj) :
195
196
197
        """keyword arguments:
        
        * region
198
            the slimPre Region object
199
        * data_proj
200
201
            data projection system.

202
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
203
        self._region = region
204
205
        if region._empty:
            return
206
        slim_private._print('Projection system for tags:')
207
        if region._physical_tags == '':
208
            slim_private._print('    the entire mesh')
209
210
        else:
            if(slim_private._is_string(region._physical_tags)):
211
                slim_private._print('    '+region._physical_tags)
212
213
            else:
                for tag in region._physical_tags:
214
                    slim_private._print('    '+tag)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
215
        mesh_proj = region._mesh._mesh_proj
Philippe Delandmeter's avatar
Philippe Delandmeter committed
216
        if mesh_proj is not None :
217
218
219
            slim_private._print('  Projection used:')
            slim_private._print('    mesh proj: %s' % mesh_proj)
            slim_private._print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
220
221
222
            _mesh_proj = dgpy.pj_init_plus(mesh_proj)
            _data_proj = dgpy.pj_init_plus(data_proj)
        elif mesh_proj is None :
223
224
225
            slim_private._print('  Projection used:')
            slim_private._print('    mesh proj: no projection. Coordinates are 3D')
            slim_private._print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
226
227
228
229
            mesh_proj2 = "+proj=latlong +ellps=WGS84"
            _mesh_proj = dgpy.pj_init_plus(mesh_proj2)
            _data_proj = dgpy.pj_init_plus(data_proj)

230
        xyz = region.coordinates      
Philippe Delandmeter's avatar
Philippe Delandmeter committed
231
        if mesh_proj :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
232
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, xyz)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
233
234
235
236
237
238
239
240
241
242
243
244
            epsVect = np.zeros(xyz.shape) 
            epsVect[:,0] = np.hypot(xyz[1,0]-xyz[0,0],xyz[1,1]-xyz[0,1]) / 1000.
            xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
            xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
            if dgpy.pj_is_latlong(_data_proj) :
                xyz1[:,0] = np.where( xyz1[:,0]<xyz0[:,0], xyz1[:,0]+2*np.pi, xyz1[:,0])
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
        else :
            lonlat = np.zeros(xyz.shape)
            lonlat[:,0] = np.arctan2(xyz[:,1], xyz[:,0])
            lonlat[:,1] = np.arcsin(xyz[:,2] / np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2]) )
            lonlat[:,2] = np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2]) 
Philippe Delandmeter's avatar
Philippe Delandmeter committed
245
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, lonlat)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
246
247
248
249
            epsVect = np.zeros(xyz.shape) 
            epsVect[:,0] = 0.01
            xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
            xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
250
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
251
252
253
254
255
256
257

        if _mesh_proj:
            dgpy.pj_free(_mesh_proj)
        if _data_proj:
            dgpy.pj_free(_data_proj)

    def rotate(self, u, v):
258
        """Return the rotated vector data array (u,v). Returned vector will be oriented in the mesh projection system
259
260
261
262
263
264
        
        keyword arguments:
        
        * u
            first component of the vector (in the data coordinate system)
        * v
265
            second component of the vector (in the data coordinate system)
266
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
267
268
269
        u2 =  u * np.cos(self._angle) + v * np.sin(self._angle)
        v2 = -u * np.sin(self._angle) + v * np.cos(self._angle)
        return (u2, v2)
270
        
271
    def rotate3D(self, F_lon, F_colat, F_r):
272
273
274
275
276
277
278
279
280
281
282
        """Return ...
        
        keyword arguments:
        
        * F_lon
            this var ...
        * F_colat
            this var ... 
        * F_r
            this var ...
        """
283
284
        lon = self.coordinates[0]
        colat = 0.5*np.pi - self.coordinates[1]
285
286
287
288
        Fx = np.sin(colat)*np.cos(lon)*F_r + np.cos(colat)*np.cos(lon)*F_colat - np.sin(lon)*F_lon
        Fy = np.sin(colat)*np.sin(lon)*F_r + np.cos(colat)*np.sin(lon)*F_colat + np.cos(colat)*F_lon
        Fz = np.cos(colat)*F_r             - np.sin(colat)*F_colat
        return (Fx,Fy,Fz)    
Philippe Delandmeter's avatar
Philippe Delandmeter committed
289
290
291


class Time :
292
293
    """Create the slimPre Time"""

294
    def __init__(self, time_vector=None, initial_time='1970-01-01 00:00:00', final_time=None, time_step=None, periodic=False) :
295
296
297
        """keyword arguments:

        * time_vector
298
            vector with all the time steps. Time steps are given in seconds from initial_time. time_vector cannot be given if final_time and/or time_step are set
299
        * initial_time
300
            initial_time given in the format 'YYYY-MM-DD hh:mm:ss'. If not set, initial_time is '1970-01-01 00:00:00'
301
        * final_time
302
            final_time given in the format 'YYYY-MM-DD hh:mm:ss'. final_time cannot be given if the time_vector is set. If final_time is set, time_step must be given.
303
        * time_step
304
305
306
307
308
            time_step in seconds. time_step cannot be given if the time_vector is set. If time_step is set, final_time must be given.
        * periodic
            flag to define a periodic variable (Default: False)
            
            If periodic, the first time should be 0 and the last time should be the period (e.g. time_vector should be from 0 to T), and the first value is repeated.
309
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
310
311
        self._periodic = periodic
        fmt = '%Y-%m-%d %H:%M:%S'
312
313
314
        if not final_time and not time_step :
            self._time = np.array(time_vector)
            fmt = '%Y-%m-%d %H:%M:%S'
315
316
            date_ref  = slim_private.datetime.datetime.strptime('1970-01-01 00:00:00', fmt)
            date_user = slim_private.datetime.datetime.strptime(initial_time, fmt)
317
318
319
            date_shift = (date_user-date_ref)
            time_shift = date_shift.days * 86400 + date_shift.seconds
            self._time += time_shift
320
        elif not time_vector and initial_time and final_time and time_step:
321
322
323
324
            date  = slim_private.datetime.datetime.strptime(initial_time, fmt)
            initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
            date  = slim_private.datetime.datetime.strptime(final_time, fmt)
            final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
325
            self._time = np.arange(initial_time, final_time+time_step*0.9999, time_step)
326
327
        else:
            dgpy.Msg.Fatal("Wrong combination of inputs to create a Time object")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
328
329


Philippe Delandmeter's avatar
Philippe Delandmeter committed
330
def write_file(file_name, region=None, time=None, data=[]):
331
    """Create a SLIM netcdf file, standard input for SLIM simulations.
332
333
334
335
336
337

    keyword arguments:

    * file_name
        name of the netcdf file to write
    * region
338
        slimPre Region object (optional).
339
    * time
340
        slimPre Time object (optional).
341
    * data
342
343
        data is a list of data tuples. Each tuple follows the form 
        
344
        ("data_name", data_array)
345
        
346
347
        where data_array is an array of n lines and m columns, n being the number of nodes in the region and m being the number of time steps in the time.
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
348
    oF = dgpy.slimNetCDFWrap(file_name)
349
350
351
    space_dep = region is not None
    time_dep = time is not None
    if space_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
352
        nbNodes = len(region._node_ids)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
353
354
        nodeDim = oF.create_dimension('node', nbNodes)
        nodeVar = oF.create_variable_int('node', [nodeDim])
355
    if time_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
356
357
358
359
360
        nTimes = np.shape(time._time)[0]
        timeArray = np.array([time._time])
        if not time._periodic:
            timeUnits = 'seconds since 1970-01-01 00:00:00'
        else:
361
            timeUnits = 'seconds since the beginning of the period'
362
363
        timeDim = oF.create_dimension('time', nTimes)
        timeVar = oF.create_variable('time', [timeDim], timeUnits, time._periodic)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
364
365
366
    nFields = len(data)
    valVar = [0]*nFields
    for iField in range(nFields) :
367
        if space_dep and time_dep:
368
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim, nodeDim], None, False)
369
370
371
        elif space_dep and (not time_dep):
            valVar[iField] = oF.create_variable(data[iField][0], [nodeDim], None, False)
        elif (not space_dep) and time_dep:
372
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim], None, False)
373
374
        else:
            valVar[iField] = oF.create_variable(data[iField][0], [], None, False)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
375
    oF.enddef()
376
377
378
379
    if space_dep :
        oF.set_variable_int(nodeVar, region._node_ids)        
    if time_dep :
        oF.set_variable(timeVar, timeArray)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
380
381
382
383
384
385
386
387
388
389
390
391
392
393
    for iField in range(nFields) :
        if not isinstance(data[iField][1], np.ndarray) :
            try :
                values = np.array([float(data[iField][1])])
            except :
                values = np.array(data[iField][1])
        else :
            values = data[iField][1]
        if len(values.shape) == 1:
            values = values[None,:]

        oF.set_variable(valVar[iField], values)


394
def netcdf_to_msh(mesh_file_name, nc_file_name, variable_name, export_file_name, time = 0):
395
    """Export a field in the netcdf format to msh format.
Philippe Delandmeter's avatar
Philippe Delandmeter committed
396
397

    keyword arguments:
398
399
400
401
402
403
404
405
406
407
408
    
    * mesh_file_name
        path to the mesh file (.msh format)
    * nc_file_name    
        path to the netcdf file (.nc format)
    * variable_name 
        name of the variable in the netcdf file
    * export_file_name 
        path to the field name (.msh format)
    * time  
        time at which the data will be exported in msh format (default: 0)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
409
410
    """
    groups = dgpy.dgGroupCollection(mesh_file_name)
411
    f = slim_private._load_function((nc_file_name,variable_name),groups)
412
    dgpy.function.getTime().set(time)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
413
414
    groups.exportFunctionMsh(f, export_file_name)

415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
def read_temporal_serie(file_name, time=None, output_netcdf=None):
    """Return a temporal serie from a text data file.
    If time is given, the function returns a vector with the data interpolated
    at the instants of time object.
    If time is not given, the function returns a tuple (time, data), where
    time is a Time object with the instants of the data file
    data is a vector with the values of the data file

    keyword arguments:
    
    * file_name
        path to the file.
        Data must be given in the following format:
            1991-06-10 00:00:00, 20.7
            1991-06-11 00:00:00, 40.2
            1991-06-15 00:00:00, 13.8
        Different date must be strictly increasing but the time step
        does not have to be constant
    * time
        time is a slimPre Time object.
        if time is given, data will be interpolated to the time instants.
    * output_netcdf
        output_netcdf = ('netcdf_file_name', 'variable_name')
        if output_netcdf is given, the data are written in a netcdf file


    """
    f = open(file_name, 'r')
    fmt = '%Y-%m-%d %H:%M:%S'
    if time:
        refTime = time._time
        valVec = np.zeros(refTime.shape)
        l = f.readline()
        (t, val) = l.split(',')
        date = slim_private.datetime.datetime.strptime(t, fmt)
        t0 = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
        val0 = float(val)
        l = f.readline()
        (t, val) = l.split(',')
        date = slim_private.datetime.datetime.strptime(t, fmt)
        t1 = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
        val1 = float(val)
        if refTime[0] < t0:
            dgpy.Msg.Fatal('Minimum time requested is before initial time from data file')
        for i in range(refTime.shape[0]):
            while refTime[i] > t1:
                t1 = t0
                val1 = val0
                l = f.readline()
                if len(l) < 10:
                    dgpy.Msg.Fatal('Requested time is after final time from data file')
                (t, val) = l.split(',')
                date = slim_private.datetime.datetime.strptime(t, fmt)
                t1 = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
                val1 = float(val)
            xsi = (refTime[i] - t0) / (t1 - t0)
            valVec[i] = (1-xsi) * val0 + xsi * val1
        if output_netcdf:
            write_file(output_netcdf[0], region=None, time=time, data=[(output_netcdf[1], valVec)])
        return valVec
    else:
        num_lines = sum(1 for line in open(file_name))
477
478
        valVec = np.zeros((num_lines))
        timeVec = np.zeros((num_lines))
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
        l = f.readline()
        it = 0
        while l != '' and it < 1e7:
            (t, val) = l.split(',')
            date = slim_private.datetime.datetime.strptime(t, fmt)
            timeVec[it] = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
            valVec[it] = float(val)
            it = it+1
            l = f.readline()
        if it > 1e7-10:
            dgpy.Msg.Fatal('temporal serie data file is too long')
        time = Time(time_vector=timeVec)
        if output_netcdf:
            write_file(output_netcdf[0], region=None, time=time, data=[(output_netcdf[1], valVec)])
        return (time, valVec)

    


Philippe Delandmeter's avatar
Philippe Delandmeter committed
498
def tpxo_tide(region, time, data_file_name=None): 
499
500
501
502
    """Give TPXO tides data
    
    This function will download the files tpxo and generate the data at the requested time steps.
    It will return the tpxo tides as variables (ssh, ssu, ssv)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
503
504
    
    keyword arguments:
505
    
506
    * region
507
        slimPre Region object on which put the data
508
    * time
509
        Time object of the instant at which put the data
510
    * data_file_name
511
        name of the netcdf data file where tides are exported (optional)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
512
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
513
    lonlat = Coordinate_system(region, data_proj='+proj=latlong +ellps=WGS84')
514
    h = u = v = np.empty((1,1))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
515
516
    h_file_name = 'h_tpxo7.2.nc'
    u_file_name = 'u_tpxo7.2.nc'
517
518
    h_file_name = slim_private.dgftp.get( ("/slim_data/tides/%s" % h_file_name) )
    u_file_name = slim_private.dgftp.get( ("/slim_data/tides/%s" % u_file_name) )
Philippe Delandmeter's avatar
Philippe Delandmeter committed
519
520
521
    tide_ssh = dgpy.slimFunctionTpxo(h_file_name,'ha','hp','lon_z','lat_z')
    tide_u = dgpy.slimFunctionTpxo(u_file_name,'Ua','up','lon_u','lat_u')
    tide_v = dgpy.slimFunctionTpxo(u_file_name,'Va','vp','lon_v','lat_v')
522
523
524
525
    if not region._empty:
        lonlat_coordinates = lonlat.coordinates * 180 / np.pi
        n = np.size(time._time)
        m = len(lonlat_coordinates)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
526

527
528
529
530
531
532
533
        h = np.empty((n,m))
        u = np.empty((n,m))
        v = np.empty((n,m))
        for i in range(n):
            h[i,:] = tide_ssh.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
            u[i,:] = tide_u.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
            v[i,:] = tide_v.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
Philippe Delandmeter's avatar
Philippe Delandmeter committed
534

535
        (u, v) = lonlat.rotate(u,v)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
536
537

    if data_file_name is not None:
Philippe Delandmeter's avatar
Philippe Delandmeter committed
538
        write_file(data_file_name, region=region, time=time, data=[('h', h), ('u', u), ('v', v)])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
539
540
541
542


    return (h, u, v)

Philippe Delandmeter's avatar
Philippe Delandmeter committed
543
def residual_flow(region, flow, bathymetry, data_file_name=None): 
544
545
546
    """Set a residual flow
    
    It will return the velocity (u, v) such that speeds at all the nodes are constant and the velocity is normal to the boundary. The total transport (in m^3 s^-1) is equal to the given flow.
547
548
549
    The physical tags of the region object must be boundary physical tags.
    
    keyword arguments:
550
    
551
    * region
552
        slimPre Region object on which put the data
553
554
555
556
557
    * flow
        value of the flow crossing the region (must be a boundary region), in m^3 s^-1
    * bathymetry
        bathymetry in netcdf format or msh format. It is necessary to compute the transport
    * data_file_name
558
        name of the netcdf data file where the residual flow is exported (optional)
559
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
560

Philippe Delandmeter's avatar
Philippe Delandmeter committed
561
    groups = region._mesh._groups
Philippe Delandmeter's avatar
Philippe Delandmeter committed
562
563
564
    integNormals = dgpy.dgFunctionIntegratorInterface(groups, dgpy.function.getNormals())
    nx = 0
    ny = 0
Philippe Delandmeter's avatar
Philippe Delandmeter committed
565
    for tag in region._physical_tags:
566
567
568
569
570
        #normals = integNormals.computeVector(tag)
        normals = dgpy.fullMatrixDouble(1,3,True)
        integNormals.compute(tag, normals)
        nx += normals.get(0,0)
        ny += normals.get(0,1)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
571
572
573
574
575
576
577
578
579
580
581
    norm = np.sqrt(nx*nx+ny*ny)
    nx = nx/norm
    ny = ny/norm

    bathFunc = slim_private._load_function(bathymetry, groups)
    def orientedSurfNumpy(cmap, val, h, n) : 
        val[:] = h[:] * (n[:,0]*nx + n[:,1]*ny)
    f = dgpy.functionNumpy(1, orientedSurfNumpy, [bathFunc, dgpy.function.getNormals()])

    integInterface = dgpy.dgFunctionIntegratorInterface(groups, f)
    section = 0
582
    sectionFM = dgpy.fullMatrixDouble(1,1,True)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
583
    for tag in region._physical_tags:
584
585
586
        #section += integInterface.compute(tag)
        integInterface.compute(tag, sectionFM)
        section += sectionFM.get(0,0)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
587
588
589
590
591

    def transportNumpy(cmap, val, h) : # eta is neglected here
        val[:,0] = -h[:] * flow / section * nx
        val[:,1] = -h[:] * flow / section * ny
    f = dgpy.functionNumpy(2, transportNumpy, [bathFunc])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
592
    vals = region._evaluateFunctor(f, 2)
593
594
    if region._empty:
        vals = np.empty((1,2))
595
    if data_file_name:
596
        write_file(data_file_name, region, None, [("uflow", vals[:,0]), ("vflow", vals[:,1])])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
597
598
    return vals

599
def smooth_bathymetry(mesh, bathymetry, output_file_name, coefficient=0.5, transform_p0=True):
600
    """Smooth the bathymetry in order to locally reduce the noise in the momentum
Philippe Delandmeter's avatar
Philippe Delandmeter committed
601
602
    
    keyword arguments:
603
604
605
    
    * mesh 
        slimPre mesh object
606
607
    * bathymetry
        bathymetry file
608
        
609
        Format:
610
611
612
613
            * tuple containing the raw bathymetry data file (.nc format) and the name of the variable bathymetry 
            * bathymetry mesh file (.msh or .idx) 
    * output_file_name
        name of the file name containing the smooth bathymetry. The format is netcdf (if transform_p0 = False) or idx (if transform_p0 = True).
614
615
616
    * coefficient
        coefficient for the smoothing algorithm (default: 0.5). This forbids the bathymetry from moving from more than 100*coefficient % on a single element
    * transform_p0
617
        flag to transform the bathymetry to a bathymetry constant per element. This gives a good stability and has negligible impact the results (default: True)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
618
619
620
    """
    groups = mesh._groups
    dof = dgpy.dgDofContainer(groups, 1)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
621
    slim_private._load(dof, bathymetry)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
622
623
624
625
626
627
628
629
630
631
    #load data
    m = groups.cgMesh()
    n = m.nVertices()
    data = np.ndarray((n),"d")
    for ig in range(groups.getNbElementGroups()) :
        g = groups.getElementGroup(ig)
        gdata = dof.getGroupProxy(ig)
        for ie in range(g.getNbElements()):
            vs = m.vertices(ig, ie)
            for i, v in enumerate(vs) :
632
633
634
                if v >= n:
                    n *= 2
                    data.resize((n))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
                data[v] = gdata(i,ie)
    #smoothing
    count = 0
    factor = coefficient
    while True:
        count = count+1
        emax = 0
        for ig in range(groups.getNbElementGroups()):
            g = groups.getElementGroup(ig)
            for ie in range(g.getNbElements()):
                vs = m.vertices(ig, ie)
                vs = np.array(vs)
                bathe = data[vs]
                dmax = bathe.max() * factor
                for i in range(len(vs)):
                    for j in range(i+1, len(vs)):
                        d = abs(bathe[i] - bathe[j])
                        if (d > dmax):
                            mean = (bathe[i] + bathe[j])/2
                            if (d - dmax > emax):
                                emax = d - dmax
                            if (bathe[i] > bathe[j]):
                                 bathe[i] = mean+dmax/2
                                 bathe[j] = mean-dmax/2
                            else:
                                 bathe[i] = mean-dmax/2
                                 bathe[j] = mean+dmax/2
                data[vs,] = bathe
663
664
        #TODO better print
        print(count,emax)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
665
666
        if emax < .2 :
            break
667
    #read data
Philippe Delandmeter's avatar
Philippe Delandmeter committed
668
669
    for ig in range(groups.getNbElementGroups()) :
        g = groups.getElementGroup(ig)
670
        gdata = dof.getGroupProxy(ig)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
671
672
673
        for ie in range(g.getNbElements()):
            vs = m.vertices(ig, ie)
            for i, v in enumerate(vs) :
674
                gdata.set(i,ie, data[v])
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
    
    if (not transform_p0):
        region = Region(mesh)
        data = region._evaluateFunctor(dof.getFunction(), 1)
        if not slim_private._is_string(bathymetry):
            if bathymetry[0][-3:] == ".nc":
                write_file(output_file_name, region=region, time=None, data=[(bathymetry[1], data)])
        else:
            write_file(output_file_name, region=region, time=None, data=[('smoothed_bath', data)])
    else:
        nodalVolume = dgpy.dgNodalVolume(groups)
        converter = dgpy.dgPrismDofConverter(nodalVolume.get())
        bath_P0 = dgpy.dgDofContainer(groups,1)
        converter.convertToPZero2d(dof,bath_P0)
        if output_file_name[-3:] == ".nc":
690
            slim_private._print("The bathymetry is dg and saved in "+output_file_name[:-3]+"/"+slim_private.path.basename(output_file_name[:-3])+".idx")
691
692
693
694
            output_file_name = output_file_name[:-3]
        exporter = dgpy.dgIdxExporter(bath_P0, output_file_name)
        exporter.exportIdx()
        
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
def interpolate_from_structured_grid(x, y, ox, oy, dx, dy, data, fill_value=-9999.):
    """Interpolates the values from a structured grid

    keyword arguments:
    
    * x
        list of the x-coordinates of all the points to interpolate at
    * y
        list of the y-coordinates of all the points to interpolate at
    * ox
        x-origin of the structured data
    * oy
        y-origin of the structured data
    * dx
        x-step of the structured data
    * dy
        y-step of the structured data
    * data
        structured grid
    * fill_value
        default value given to a point outside the structured grid
    """    
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
    ny = data.shape[0]
    nx = data.shape[1]
    rx = x - ox
    ry = y - oy
    ix = (rx / dx).astype(int)
    iy = (ry / dy).astype(int)
    index=np.empty(x.shape)
    if (dx > 0):
        index = (x < ox) + (x > ox + (nx-1)*dx)
    else:
        index = (x > ox) + (x < ox + (nx-1)*dx)
    if (dy > 0):
        index+= (y < oy) + (y > oy + (ny-1)*dy)
    else:
        index+= (y > oy) + (y < oy + (ny-1)*dy)
    ix = np.minimum(np.maximum(ix,0),nx-1)
    iy = np.minimum(np.maximum(iy,0),ny-1)
    xsi = (rx - ix*dx) / dx
    eta = (ry - iy*dy) / dy
736
    val = np.where(index, fill_value, \
737
738
739
740
741
742
          (1-xsi)*(1-eta) * data[ iy      ,  ix      ] + \
          (  xsi)*(1-eta) * data[ iy      , (ix+1)%nx] + \
          (  xsi)*(  eta) * data[(iy+1)%ny, (ix+1)%nx] + \
          (1-xsi)*(  eta) * data[(iy+1)%ny,  ix      ])
    return val
        
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
def interpolate_from_structured_file(nodes, structured_file, fill_value=-9999, transpose_matrix=False):
    """Interpolates the values from a structured file of gmsh type:
    
    ox oy oz
    
    dx dy dz
    
    nx ny nz
    
    val11 val12 ...
    
    val21 val22 ...
    
    ...

758
759
    This module is only available in 2D !

760
761
762
763
764
765
766
767
768
769
770
    keyword arguments:
    
    * nodes
        coordinates of the nodes to interpolate at (from Region or Coordinate_system)
    * structured_file
        name of the file containing the structured data (binary files should end with '.bin' extension) 
    * fill_value
        default value given to a point outside the structured grid
    * transpose_matrix
        boolean assessing if the matrix should be transposed (Default: False) (True if nx is the number of lines and ny is the number of columns)
    """   
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
    if structured_file[-4:] == '.bin':
        import struct
        f = open(structured_file,'rb')
        O = struct.unpack("ddd", f.read(24))
        ox = O[0]
        oy = O[1]
        D = struct.unpack("ddd", f.read(24))
        dx = D[0]
        dy = D[1]
        N = struct.unpack("iii", f.read(12))
        nx = N[0]
        ny = N[1]
        if transpose_matrix:
            data = np.fromfile(f, "d", nx * ny).reshape((nx, ny))
        else:
            data = np.fromfile(f, "d", nx * ny).reshape((ny, nx))
    else:
        f = open(structured_file,'r')
        line = f.readline()
        origin = line.split(' ')
        ox = float(origin[0])
        oy = float(origin[1])
        line = f.readline()
        delta = line.split(' ')
        dx = float(delta[0])
        dy = float(delta[1])
        line = f.readline()
        number = line.split(' ')
        nx = int(number[0])
        ny = int(number[1])
        if transpose_matrix:
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((nx, ny))
        else:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
804
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((ny, nx))
805
806
807
808
    f.close()
    x = nodes[:,0]
    y = nodes[:,1]
    data = np.transpose(data)
809
    val = interpolate_from_structured_grid(x, y, ox, oy, dx, dy, data, fill_value=fill_value)    
810
811
    return val
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
812
def get_data_from_netcdf(nc_file_name, variable_name):
813
    """Return the array containing the values of the variable in the nectdf file (SLIM format)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
814
815

    keyword arguments:
816
    
817
818
819
820
    * nc_file_name
        path to the netcdf file (.nc format)
    * variable_name
        name of the variable in the netcdf file
Philippe Delandmeter's avatar
Philippe Delandmeter committed
821
822
    """
    #TODO use a wrapper instead of netCDF4
823
    slim_private._print("Loading netCDF4 for python, might not be installed on your computer...")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
824
    from netCDF4 import Dataset
825
    f = Dataset(nc_file_name,"r")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
826
827
828
    val = f.variables[variable_name][:]
    return val

829
def get_distance(region, physical_tags):
830
    """Return the minimum distance from each node on the region to the physical tags. The index of the 1-D array returned is the nodal index, same as Region.coordinates 
831
832
    
    keyword arguments:
833
    
834
    * region
835
        slimPre Region object of the nodes on which compute the distance
836
837
    * physical_tags
        vector containing the physical tags from which the minimum distance will be computed
838
839
    """
    xyz = region.coordinates
840
    dist = slim_private.distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
841
842
    return np.array(dist)

843
844
845
846
847
848
849
850
851
852
853
854
def extrude(mesh_file_name, bath_file_name, nb_layers=None, z_layers=None, layers_function=None, mesh_file_name_out='', factor_show=0, periodicity=None):


    if (nb_layers and z_layers) or (nb_layers and layers_function) or (z_layers and layers_function):
        dgpy.Msg.Fatal("Only one among nb_layers, z_layers and layers_function can be defined!")

    if periodicity:
        shiftOperation = periodicity[0]
        cutTags = periodicity[1]
        pasteTags = periodicity[2]
        mapFilename = periodicity[3]
        slim_private.periodicMap.periodicMap(mesh_file_name, '_tmp_mesh_file', 1, shiftOperation, cutTags, pasteTags, mapFilename)
855
        slim_private.shutil.move('_tmp_mesh_file', mesh_file_name)
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911

    groups = dgpy.dgGroupCollection(mesh_file_name)
    bathFunc = slim_private._load_function(bath_file_name, groups)
    bathDof = dgpy.dgDofContainer(groups, 1)
    bathDof.interpolate(bathFunc)
    
    globNodeId2Bath = {}
    dataProx = dgpy.fullMatrixDouble()
    for iGroup in range(groups.getNbElementGroups()) :
        group = groups.getElementGroup(iGroup)
        for iElement in range(group.getNbElements()) :
            el = group.getElement(iElement)
            bathDof.getElementProxy(iGroup,iElement, dataProx)
            for iPt in range(group.getNbNodes()) :
                vId = groups.mesh().gmshVertexId(el.vertex(iPt))
                globNodeId2Bath[vId] = dataProx(iPt,0)

    def layersForZ(h) :
        for i in range(len(z_layers)) :
            if z_layers[i] > h : 
                return i
        return len(z_layers)

    def getLayers(e, v) :
        h = globNodeId2Bath[v[3]]
        if nb_layers:
            return [z * h /nb_layers for z in  range(nb_layers + 1)]
        elif z_layers:
            return [z_layers[i] for i in range(layersForZ(h))]
        elif layers_function:
            return layers_function(v[0], v[1], h)
        else:
            dgpy.Msg.Fatal("Neither nb_layers nor z_layers nor layers_function was defined")

    def getLayersShow(e, v) :
        layers = getLayers(e, v)
        for i in range(len(layers)):
            layers[i] = factor_show * layers[i]
        return layers

    if mesh_file_name_out:
        mesh3d = mesh_file_name_out
        mesh3dShow = mesh_file_name_out[:-4] + '_show.msh'
    else:
        mesh3d = mesh_file_name[:-4] + '_3d.msh'
        mesh3dShow = mesh_file_name[:-4] + '_3d_show.msh'

    if periodicity:
        slim_private.extrudePeriodic.extrudePeriodic(slim_private.mesh.mesh(mesh_file_name), getLayers, mapFilename, mapFilename).write(mesh3d)
    else:
        slim_private.extrude.extrude(slim_private.mesh.mesh(mesh_file_name), getLayers).write(mesh3d)
    if factor_show > 0:
        slim_private.extrude.extrude(slim_private.mesh.mesh(mesh_file_name), getLayersShow).write(mesh3dShow)
    return


912
def fetch_ftp(file_name, data_dir=None):
913
914
    """This function downloads some data from the SLIM servers if available. The data is downloaded in the user global_data_directory. Data is only downloaded if it does not exist locally or is older than the server data. 

915
916
    global_data_directory is by default ~/.cache/slim_data/
    """
917
    file_name =  slim_private.dgftp.get(file_name)
918
919
    if data_dir: 
        if dgpy.Msg.GetCommRank() == 0:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
920
            slim_private.shutil.copyfile(file_name, data_dir+slim_private.path.basename(file_name))
921
            file_name = data_dir+slim_private.path.basename(file_name)
922
923
924
        dgpy.Msg.Barrier()
    return file_name

925
def set_global_data_directory(path):
926
    """Edit the global_data_directory"""
927
928
    slim_private.dgftp.global_data_local_directory = path

929
def make_directory(data_dir):
930
    """Create data_dir directory if it does not exist yet"""
931
932
    if not slim_private.path.isdir(data_dir) :
        slim_private.makedirs(data_dir)
933
    dgpy.Msg.Barrier()
934
    
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
def partition_mesh(mesh_file_name, partitioned_mesh_file_name):
    """Partition the mesh 
    
    keyword arguments:
    
    * mesh_file_name            
            path to the mesh file (.msh format)
    * partitioned_mesh_file_name
            name where the partitioned mesh is stored.
    """
    isPart = slim_private._check_if_partionned(mesh_file_name)
    if isPart:
        dgpy.Msg.Warning("The mesh "+mesh_file_name+" is already partitioned !")
        slim_private.shutil.copyfile(mesh_file_name, partitioned_mesh_file_name)
    elif dgpy.Msg.GetCommSize() == 1:
        dgpy.Msg.Info("No partition done on a single processor, the partition will be done later, if needed")
        slim_private.shutil.copyfile(mesh_file_name, partitioned_mesh_file_name)
    else:
        dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
        mesh_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'
        slim_private.shutil.copyfile(mesh_name, partitioned_mesh_file_name)          
        dgpy.Msg.Barrier()
        if dgpy.Msg.GetCommRank() == 0 :
            slim_private.remove(mesh_name)
    return partitioned_mesh_file_name

961
962
963
964
965
966
967
968
969
def partition_id():
    """Return the number of the current partition as a string
    """
    return str(dgpy.Msg.GetCommRank())

def partition_nb():
    """Return the number of partitions used as a string
    """
    return str(dgpy.Msg.GetCommSize())
970
971

def exit(i):
972
    """Exit properly the program"""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
973
    dgpy.Msg.Exit(i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997

def mesh_shp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection = "+latlong +ellps=WGS84",  boundary_names={}):
    """Generate a mesh file in gmsh .msh file format (http://gmsh.info/). The mesh is generated on a sphere of radius 1. The mesh size
    should be specified on this sphere (i.e. distances should be divided by the earth radius)

    arguments:

    * shp_inputs : list of string
        vector of path to files in the ESRI shape file format. The shape files can be composed of point, lines or polygons. The projection
        used in the shape files should be longitude-latitude (in degrees).
    * lon, lat : float
        position of one point inside the domain, the mesh region is the close surface containing this point.
    * mesh_size : (float x, float y, float z) -> float
        callback to get the desired element size, (on a sphere of radius 1) at one point on this sphere
    * mesh_size_min : float
        minimum value of the mesh size (approximatively the minimum value returned by mesh_size) on a sphere of radius 1
    * output : string
        name of the output file
    * projection : string
        description of the projection used in the output file in the libproj "+" format (the mesh itself is generated on the sphere)
    * boundary_name : dictionary of (int, string)
        to assign physicals names to the boundary elements of the mesh, a field (of type integer64) called "entity" should be present in the shape files. Then this parameter can be used to map the values of this field to physical names
    """
    dgpy.meshShp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection, boundary_names)