slimPre.py

"""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.

    slimPre contains the tools writing the input files necessary for the simulations (for simulations, see slim.py file), such as partitionned mesh, bathymetry, coriolis, wind, tides, ...
    
    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)

    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:
        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.
        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
    
    netcdf_to_msh:
        The function transforms a netcdf file to a .msh file.

    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:
        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.

    get_data_from_netcdf:
        The function returns the data array contained in the netcdf file given as an argument

    get_distance:
        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. 

    fetch_ftp:
        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/

    set_global_data_directory:
        This functions gives the possibility to edit the global_data_directory

    make_directory:
        The function makes a directory if it does not exist
    exit:
        The function closes the program
"""
import dgpy
import slim_private
import numpy as np

class Mesh :
    """Create the slimPre Mesh"""
    def __init__(self, mesh_file_name, mesh_proj=None) :
        """keyword arguments:
        
        * mesh_file_name            
            path to the mesh file (.msh format). The partitioned mesh will automatically be loaded in multiprocessing
        * mesh_proj           
            mesh projection system (Default: None). 
            This argument must not be given if the mesh is writen in 3D (on a sphere), or if all the data are in the same projection system than the mesh.
        """
        if dgpy.Msg.GetCommSize() > 1 :
            dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
            mesh_file_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'
        self._groups = dgpy.dgGroupCollection(mesh_file_name)
        self._mesh_proj = mesh_proj
        

class Region :
    """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
    """
    def __init__(self, mesh, physical_tags='') :
        """keyword arguments:
        
        * mesh
            the slimPre Mesh object
        * physical_tags
            the list of all the tags on which put some data (optional)
            
            If physical_tags are not set, all the mesh nodes are gathered
        """
        self._mesh = mesh
        if slim_private._is_string(physical_tags):
            physical_tags = [physical_tags]
        self._physical_tags = physical_tags
        self._nodes = []
        self._groups = []
        self._node_map = {}
        self._node_ids = []
        self._empty = False
        num_ = 0
        self._type = "interface"
        for iGroup in range(mesh._groups.getNbElementGroups()+mesh._groups.getNbGhostGroups()) :
            group = mesh._groups.getElementGroup(iGroup)
            if ((self._physical_tags != ['']) and (not group.getPhysicalTag() in self._physical_tags)) :
                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" :
            self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
            return

        for iFaceGroup in range(mesh._groups.getNbFaceGroups()) :
            faceGroup = mesh._groups.getFaceGroup(iFaceGroup)
            if (not faceGroup.physicalTag() in self._physical_tags) :
                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)
        if len(self._nodes) == 0:
            self._empty = True

    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 :
    """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)
        Coordinate_system.rotate3D ... 
    """
    def __init__(self, region, data_proj) :
        """keyword arguments:
        
        * region
            the slimPre Region object
        * data_proj
            data projection system.

        """
        self._region = region
        if region._empty:
            return
        print('Projection system for tags:')
        if region._physical_tags == '':
            print('    the entire mesh')
        else:
            if(slim_private._is_string(region._physical_tags)):
                print('    '+region._physical_tags)
            else:
                for tag in region._physical_tags:
                    print('    '+tag)
        mesh_proj = region._mesh._mesh_proj
        if mesh_proj is not None :
            print('  Projection used:')
            print('    mesh proj: %s' % mesh_proj)
            print('    data proj: %s' % data_proj)
            _mesh_proj = dgpy.pj_init_plus(mesh_proj)
            _data_proj = dgpy.pj_init_plus(data_proj)
        elif mesh_proj is None :
            print('  Projection used:')
            print('    mesh proj: no projection. Coordinates are 3D')
            print('    data proj: %s' % data_proj)
            mesh_proj2 = "+proj=latlong +ellps=WGS84"
            _mesh_proj = dgpy.pj_init_plus(mesh_proj2)
            _data_proj = dgpy.pj_init_plus(data_proj)

        xyz = region.coordinates      
        if mesh_proj :
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, xyz)
            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]) 
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, lonlat)
            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)
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])

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

    def rotate(self, u, v):
        """Return the rotated vector data array (u,v). Returned vector will be oriented in the mesh projection system
        
        keyword arguments:
        
        * u
            first component of the vector (in the data coordinate system)
        * v
            second component of the vector (in the data coordinate system)
        """
        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)
        
    def rotate3D(self, F_lon, F_colat, F_r):
        """Return ...
        
        keyword arguments:
        
        * F_lon
            this var ...
        * F_colat
            this var ... 
        * F_r
            this var ...
        """
        lon = self.coordinates[0]
        colat = 0.5*np.pi - self.coordinates[1]
        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)    


class Time :
    """Create the slimPre Time"""

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

        * time_vector
            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
        * initial_time
            initial_time given in the format 'YYYY-MM-DD hh:mm:ss'. If not set, initial_time is '1970-01-01 00:00:00'
        * final_time
            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.
        * time_step
            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.
        """
        self._periodic = periodic
        fmt = '%Y-%m-%d %H:%M:%S'
        if not final_time and not time_step :
            self._time = np.array(time_vector)
            fmt = '%Y-%m-%d %H:%M:%S'
            date_ref  = slim_private.datetime.datetime.strptime('1970-01-01 00:00:00', fmt)
            date_user = slim_private.datetime.datetime.strptime(initial_time, fmt)
            date_shift = (date_user-date_ref)
            time_shift = date_shift.days * 86400 + date_shift.seconds
            self._time += time_shift
        elif not time_vector and initial_time and final_time and time_step:
            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])
            self._time = np.arange(initial_time, final_time+time_step*0.9999, time_step)
        else:
            dgpy.Msg.Fatal("Wrong combination of inputs to create a Time object")


def write_file(file_name, region=None, time=None, data=[]):
    """Create a SLIM netcdf file, standard input for SLIM simulations.

    keyword arguments:

    * file_name
        name of the netcdf file to write
    * region
        slimPre Region object (optional).
    * time
        slimPre Time object (optional).
    * data
        data is a list of data tuples. Each tuple follows the form 
        
        ("data_name", data_array)
        
        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.
    """
    if dgpy.Msg.GetCommSize() > 1 :
        file_name = file_name+"_"+str(dgpy.Msg.GetCommRank())
    oF = dgpy.slimNetCDFWrap(file_name)
    space_dep = region is not None
    time_dep = time is not None
    if space_dep :
        nbNodes = len(region._node_ids)
        nodeDim = oF.create_dimension('node', nbNodes)
        nodeVar = oF.create_variable_int('node', [nodeDim])
    if time_dep :
        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:
            timeUnits = 'seconds since the beginning of the period'
        timeDim = oF.create_dimension('time', nTimes)
        timeVar = oF.create_variable('time', [timeDim], timeUnits, time._periodic)
    nFields = len(data)
    valVar = [0]*nFields
    for iField in range(nFields) :
        if space_dep and time_dep:
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim, nodeDim], None, False)
        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:
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim], None, False)
        else:
            valVar[iField] = oF.create_variable(data[iField][0], [], None, False)
    oF.enddef()
    if space_dep :
        oF.set_variable_int(nodeVar, region._node_ids)        
    if time_dep :
        oF.set_variable(timeVar, timeArray)
    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)


def netcdf_to_msh(mesh_file_name, nc_file_name, variable_name, export_file_name, time = 0):
    """Export a field in the netcdf format to msh format.

    keyword arguments:
    
    * 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)
    """

    if dgpy.Msg.GetCommSize() > 1 :
        dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
        mesh_file_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'

    groups = dgpy.dgGroupCollection(mesh_file_name)
    if dgpy.Msg.GetCommSize() == 1 :
        f = dgpy.slimInputRead(nc_file_name, variable_name)
    else :
        f = dgpy.slimInputRead(nc_file_name+'_'+str(dgpy.Msg.GetCommRank()), variable_name)
    dgpy.function.getTime().set(time)
    groups.exportFunctionMsh(f, export_file_name)

def tpxo_tide(region, time, data_file_name=None): 
    """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)
    
    keyword arguments:
    
    * region
        slimPre Region object on which put the data
    * time
        Time object of the instant at which put the data
    * data_file_name
        name of the netcdf data file where tides are exported (optional)
    """
    lonlat = Coordinate_system(region, data_proj='+proj=latlong +ellps=WGS84')
    h = u = v = np.empty((1,1))
    h_file_name = 'h_tpxo7.2.nc'
    u_file_name = 'u_tpxo7.2.nc'
    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) )
    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')
    if not region._empty:
        lonlat_coordinates = lonlat.coordinates * 180 / np.pi
        n = np.size(time._time)
        m = len(lonlat_coordinates)

        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]

        (u, v) = lonlat.rotate(u,v)

    if data_file_name is not None:
        write_file(data_file_name, region=region, time=time, data=[('h', h), ('u', u), ('v', v)])


    return (h, u, v)

def residual_flow(region, flow, bathymetry, data_file_name=None): 
    """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.
    The physical tags of the region object must be boundary physical tags.
    
    keyword arguments:
    
    * region
        slimPre Region object on which put the data
    * 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
        name of the netcdf data file where the residual flow is exported (optional)
    """

    groups = region._mesh._groups
    integNormals = dgpy.dgFunctionIntegratorInterface(groups, dgpy.function.getNormals())
    nx = 0
    ny = 0
    for tag in region._physical_tags:
        #normals = integNormals.computeVector(tag)
        normals = dgpy.fullMatrixDouble(1,3,True)
        integNormals.compute(tag, normals)
        nx += normals.get(0,0)
        ny += normals.get(0,1)
    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
    sectionFM = dgpy.fullMatrixDouble(1,1,True)
    for tag in region._physical_tags:
        #section += integInterface.compute(tag)
        integInterface.compute(tag, sectionFM)
        section += sectionFM.get(0,0)

    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])
    vals = region._evaluateFunctor(f, 2)
    if region._empty:
        vals = np.empty((1,2))
    if data_file_name:
        write_file(data_file_name, region, None, [("uflow", vals[:,0]), ("vflow", vals[:,1])])
    return vals

def smooth_bathymetry(mesh, bathymetry, output_file_name, coefficient=0.5, transform_p0=True):
    """Smooth the bathymetry in order to locally reduce the noise in the momentum
    
    keyword arguments:
    
    * mesh 
        slimPre mesh object
    * bathymetry
        bathymetry file
        
        Format:
            * 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).
    * 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
        flag to transform the bathymetry to a bathymetry constant per element. This gives a good stability and has negligible impact the results (default: True)
    """
    groups = mesh._groups
    dof = dgpy.dgDofContainer(groups, 1)
    slim_private._load(dof, bathymetry)
    #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) :
                if v >= n:
                    n *= 2
                    data.resize((n))
                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
        #TODO better print
        print(count,emax)
        if emax < .2 :
            break
    #read data
    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) :
                gdata.set(i,ie, data[v])
    
    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":
            dgpy.Msg.Info("The bathymetry is dg and saved in "+output_file_name[:-3]+"/"+slim_private.path.basename(output_file_name[:-3])+".idx")
            output_file_name = output_file_name[:-3]
        exporter = dgpy.dgIdxExporter(bath_P0, output_file_name)
        exporter.exportIdx()
        
def interpolate_on_structured_grid(x, y, ox, oy, dx, dy, data, fill_vallue=-9999.):
    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
    val = np.where(index, fill_vallue, \
          (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
        
def interpolate_from_structured_grid(nodes, structured_file, fill_value=-9999, transpose_matrix=False):
    """
    """
    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:
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((ny, nx))        
        np.save("data.npy",data)
    f.close()
    x = nodes[:,0]
    y = nodes[:,1]
    data = np.transpose(data)
    val = interpolate_on_structured_grid(x, y, ox, oy, dx, dy, data, fill_vallue=fill_value)    
    return val
        
def get_data_from_netcdf(nc_file_name, variable_name):
    """Return the array containing the values of the variable in the nectdf file (SLIM format)

    keyword arguments:
    
    * nc_file_name
        path to the netcdf file (.nc format)
    * variable_name
        name of the variable in the netcdf file
    """
    #TODO use a wrapper instead of netCDF4
    dgpy.Msg.Info("Loading netCDF4 for python, might not be installed on your computer...")
    from netCDF4 import Dataset
    if (dgpy.Msg.GetCommSize() > 1):
        f = Dataset(nc_file_name+'_'+str(dgpy.Msg.GetCommRank()),"r")
    else:
        f = Dataset(nc_file_name,"r")
    val = f.variables[variable_name][:]
    return val

def get_distance(region, physical_tags):
    """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 
    
    keyword arguments:
    
    * region
        slimPre Region object of the nodes on which compute the distance
    * physical_tags
        vector containing the physical tags from which the minimum distance will be computed
    """
    xyz = region.coordinates
    dist = slim_private.distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
    return np.array(dist)

def fetch_ftp(file_name, data_dir=None):
    """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/
    """
    file_name =  slim_private.dgftp.get(file_name)
    if data_dir: 
        if dgpy.Msg.GetCommRank() == 0:
            slim_private.shutil.copyfile(file_name, data_dir+slim_private.path.basename(file_name))
            file_name = data_dir+slim_private.path.basename(file_name)
        dgpy.Msg.Barrier()
    return file_name

def set_global_data_directory(path):
    """Edit the global_data_directory"""
    slim_private.dgftp.global_data_local_directory = path

def make_directory(data_dir):
    """Create data_dir directory if it does not exist yet"""
    if dgpy.Msg.GetCommRank() == 0:
        if not slim_private.path.isdir(data_dir) :
            slim_private.makedirs(data_dir)
    dgpy.Msg.Barrier()

def exit(i):
    """Exit properly the program"""
    dgpy.Msg.Exit(i)