gmsh_io.h

#pragma once

#include <vector>
#include <array>
#include <gmsh.h>
#include <type_traits>

#ifdef USE_MPI
#include "gmsh_mpi.h"
#endif /* USE_MPI */

#include "entity.h"
#include "physical_element.h"

namespace g     = gmsh;
namespace go    = gmsh::option;
namespace gm    = gmsh::model;
namespace gmm   = gmsh::model::mesh;
namespace gmo   = gmsh::model::occ;
namespace gmg   = gmsh::model::geo;
namespace gv    = gmsh::view;

using gvp_t = gmsh::vectorpair;

/* Macros used to annotate GMSH integer inputs */
#define DIMENSION(x) x

std::string quadrature_name(int);
std::string basis_func_name(int);
std::string basis_grad_name(int);

using face_key = std::array<size_t, 3>;

enum class face_type : int
{
    NONE = 0,
    INTERFACE,
    BOUNDARY,
#ifdef USE_MPI
    INTERPROCESS_BOUNDARY
#endif
};

struct boundary_descriptor
{
    face_type   type;
    int         gmsh_entity;        /* Actual GMSH entity */
#ifdef USE_MPI
    int         parent_entity;      /* Parent GMSH entity in partitioned model */
#endif /* USE_MPI */

    bool operator<(const boundary_descriptor& other) const
    {
        return (type < other.type) or
               ( (type == other.type) and (gmsh_entity < other.gmsh_entity) );
    }

    bool operator==(const boundary_descriptor& other) const
    {
        return (type == other.type) and (gmsh_entity == other.gmsh_entity);
    }

    boundary_descriptor()
        : type(face_type::NONE), gmsh_entity(0)
#ifdef USE_MPI
        , parent_entity(-1)
#endif /* USE_MPI */
    {}

    int material_tag(void) const
    {
#ifdef USE_MPI
        if (parent_entity != -1)
            return parent_entity;
#endif
        return gmsh_entity;
    }
};



#define IMPLIES(a, b) ((not(a)) or (b))
#define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))

#ifdef USE_MPI
struct dof_range
{
    int     tag;
    int     dim;
    size_t  base;
    size_t  length;
};

inline std::ostream& operator<<(std::ostream& os, const dof_range& dr)
{
    os << "Tag " << dr.tag << ", dim " << dr.dim << ": (";
    os << dr.base << ", " << dr.length << ")";
    return os;
}
#endif

class model
{
    int                                         geometric_order;
    int                                         approximation_order;
#ifdef USE_MPI
    int                                         num_partitions;
#endif /* USE_MPI */

    /* Map from boundary tag to all the faces of the entity with that tag */
    std::map<int, std::vector<element_key>>     boundary_map;   

    std::vector<entity>                         entities;
    std::vector<boundary_descriptor>            bnd_descriptors;
    element_connectivity<element_key>           conn;

#ifdef USE_MPI
    /* Map from boundary tag to all the partitions it belongs to */
    std::map<int, std::vector<int>>             comm_map;       

    /* Map from the partition number to the entities belonging to that partition */
    std::map<int, std::vector<int>>             partition_map;
    
    /* In a partitioned model, map from 2D entity tag to parent tag */
    std::map<int, int>                          surface_to_parent_map;

    /* Map from 3D entity tag to partition */
    std::map<int, int>                          partition_inverse_map;

    /* Global dof range, only rank0 has this */
    std::map<int, dof_range>                    global_dof_ranges;
#endif /* USE_MPI */

    using entofs_pair = std::pair<size_t, size_t>;
    std::map<size_t, entofs_pair>               etag_to_entity_offset;

    void    map_boundaries(void);
    void    import_gmsh_entities_rank0(void);
    void    update_connectivity(const entity&, size_t);

    void    populate_from_gmsh_rank0(void);

#ifdef USE_MPI
    bool    is_interprocess_boundary(int);
    void    populate_from_gmsh_rankN(void);
    void    import_gmsh_entities_rankN(void);
    int     my_partition(void) const;
    void    make_partition_to_entities_map(void);
#endif /* USE_MPI */

    void    make_boundary_to_faces_map(void);

public:
    model();
    model(int, int);
    model(const char *);
    model(const char *, int, int);
    ~model();

    void build(void);
    void build(double);

    void generate_mesh(void);
    void generate_mesh(double);
    void partition_mesh(int);
    void populate_from_gmsh(void);

    const element_connectivity<element_key>& connectivity() const
    {
        return conn;
    }

    const std::vector<boundary_descriptor>& boundary_descriptors() const
    {
        return bnd_descriptors;
    }

    std::vector<element_key> get_bnd(size_t which)
    {
        return boundary_map.at(which);
    }

    size_t  num_dofs() const;
    size_t  num_fluxes() const;
    size_t  num_cells() const;
    size_t  num_faces() const;
    size_t  num_entities() const;

    std::vector<entity>::const_iterator begin() const;
    std::vector<entity>::const_iterator end() const;
    std::vector<entity>::iterator begin();
    std::vector<entity>::iterator end();

    entity&         entity_at(size_t);
    const entity&   entity_at(size_t) const;

    entofs_pair     lookup_tag(size_t tag) const;
};