mesh.h

Go to the documentation of this file.

 
#ifndef MESH_H
#define MESH_H
 
#include <array>
#include <map>
#include <memory>
#include <string>
#include <vector>
#include <unordered_map>
 
#include <Kokkos_Core.hpp>
#include <toml.hpp>
 
#include "curved.h"
#include "zone.h"
 
enum class MeshType {
    FROM_FILE,
    CARTESIAN,
    CARTESIAN_TRI,
    WEDGE,
    CARTESIAN_TET,
    CARTESIAN_PRISM,
    CARTESIAN_PYRAMID,
    CARTESIAN_MIXED
};
 
static const std::unordered_map<std::string, MeshType> MESH_TYPES = {
    {"file", MeshType::FROM_FILE},
    {"cartesian", MeshType::CARTESIAN},
    {"cartesian_tri", MeshType::CARTESIAN_TRI},
    {"wedge", MeshType::WEDGE},
    {"cartesian_tet", MeshType::CARTESIAN_TET},
    {"cartesian_prism", MeshType::CARTESIAN_PRISM},
    {"cartesian_pyramid", MeshType::CARTESIAN_PYRAMID},
    {"cartesian_mixed", MeshType::CARTESIAN_MIXED}
};
 
static const std::unordered_map<MeshType, std::string> MESH_NAMES = {
    {MeshType::FROM_FILE, "file"},
    {MeshType::CARTESIAN, "cartesian"},
    {MeshType::CARTESIAN_TRI, "cartesian_tri"},
    {MeshType::WEDGE, "wedge"},
    {MeshType::CARTESIAN_TET, "cartesian_tet"},
    {MeshType::CARTESIAN_PRISM, "cartesian_prism"},
    {MeshType::CARTESIAN_PYRAMID, "cartesian_pyramid"},
    {MeshType::CARTESIAN_MIXED, "cartesian_mixed"}
};
 
using Stretching = std::array<rtype, 3>;
 
enum class CellType {
    TRIANGLE,
    QUAD,
    TETRAHEDRON,
    PYRAMID,
    PRISM,
    HEXAHEDRON
};
 
static const std::unordered_map<std::string, CellType> CELL_TYPES = {
    {"triangle", CellType::TRIANGLE},
    {"quad", CellType::QUAD},
    {"tetrahedron", CellType::TETRAHEDRON},
    {"pyramid", CellType::PYRAMID},
    {"prism", CellType::PRISM},
    {"hexahedron", CellType::HEXAHEDRON}
};
 
static const std::unordered_map<CellType, std::string> CELL_NAMES = {
    {CellType::TRIANGLE, "triangle"},
    {CellType::QUAD, "quad"},
    {CellType::TETRAHEDRON, "tetrahedron"},
    {CellType::PYRAMID, "pyramid"},
    {CellType::PRISM, "prism"},
    {CellType::HEXAHEDRON, "hexahedron"}
};
 
const std::vector<std::vector<uint8_t>> & cell_local_faces(uint32_t n_nodes);
 
void cell_tetrahedra(const std::vector<std::array<double, 3>> & nodes,
                     std::vector<std::array<std::array<double, 3>, 4>> & tets);
 
 
struct MeshBlock;
struct PeriodicNodes;
 
class Mesh {
    public:
        Mesh();
 
        ~Mesh();
 
        struct PeriodicPair {
            std::string zone_a, zone_b;
            std::array<rtype, N_DIM> translation;
        };
 
        void init(const toml::value & input);
 
        static std::vector<PeriodicPair> periodic_pairs(const toml::value & input);
 
        MeshType get_type() const;
 
        void set_type(MeshType type);
 
        uint32_t n_face_zones() const;
 
        std::vector<FaceZone> * face_zones();
 
        FaceZone * get_face_zone(const std::string& name);
 
        uint32_t h_n_nodes_of_cell(uint32_t i_cell) const;
 
        uint32_t h_n_faces_of_cell(uint32_t i_cell) const;
 
        uint32_t h_n_nodes_of_face(uint32_t i_face) const;
 
        uint32_t h_node_of_cell(uint32_t i_cell, uint32_t i_node_local) const;
 
        uint32_t h_face_of_cell(uint32_t i_cell, uint32_t i_face_local) const;
 
        uint32_t h_node_of_face(uint32_t i_face, uint32_t i_node_local) const;
 
        CellType h_cell_type(uint32_t i_cell) const;
 
        void h_neighbors_of_cell(uint32_t i_cell,
                                 uint8_t n_order,
                                 std::vector<uint32_t> & neighbors) const;
 
        void compute_cell_centroids();
 
        void compute_cell_volumes();
 
        void compute_face_areas();
 
        void compute_face_normals();
 
        void compute_face_centroids();
 
        void compute_geometry();
 
        void make_axisymmetric();
 
        void h_cell_tetrahedra(uint32_t i_cell, std::vector<std::array<std::array<double, 3>, 4>> & tets) const;
 
        void compute_cell_neighbors();
 
        void copy_host_to_device();
 
        void copy_device_to_host();
 
        void init_cart(uint32_t nx, uint32_t ny, rtype Lx, rtype Ly);
 
        void init_cart_tri(uint32_t nx, uint32_t ny, rtype Lx, rtype Ly);
 
        void init_cart_3d(uint32_t nx, uint32_t ny, uint32_t nz, rtype Lx, rtype Ly, rtype Lz, MeshType kind,
                          const std::vector<PeriodicPair> & periodic = {}, const Stretching & stretching = {});
 
        struct BoundaryFace {
            std::vector<uint32_t> nodes;
            std::string zone;
        };
 
        void init_from_connectivity(const std::vector<std::array<rtype, N_DIM>> & nodes,
                                    const std::vector<std::vector<uint32_t>> & cells,
                                    const std::vector<BoundaryFace> & boundary_faces,
                                    const std::string & unlisted_zone = "unassigned",
                                    const std::vector<PeriodicPair> & periodic = {});
 
        void init_from_connectivity(const std::vector<std::array<rtype, N_DIM>> & nodes,
                                    const std::vector<std::vector<uint32_t>> & cells,
                                    const std::vector<BoundaryFace> & boundary_faces,
                                    const std::string & unlisted_zone, const PeriodicNodes & classes);
 
        void init_from_block(const MeshBlock & block, const std::vector<PeriodicPair> & periodic = {});
 
        void apply_curved(const curved::Surfaces & surfaces);
 
        const std::vector<FaceZone> & h_face_zones() const { return m_face_zones; }
 
        bool is_periodic() const { return !periodic_translations.empty(); }
 
        rtype h_face_offset(uint32_t f, uint32_t c, int i) const {
            return (h_cells_of_face(f, 1) == int32_t(c)) ? h_shifts(h_face_shift(f), i) : rtype(0);
        }
 
        void init_file(const std::string & filename, const std::vector<PeriodicPair> & periodic = {});
 
        void init_wedge(uint32_t nx, uint32_t ny, rtype Lx, rtype Ly);
 
        uint64_t h_global_cell(uint32_t i_cell) const {
            return h_global_cell_id.empty() ? i_cell : h_global_cell_id[i_cell];
        }
 
        uint32_t n_cells, n_nodes, n_faces;
        // High-order boundary elements of the mesh file (global node ids), for curved::surfaces_of_mesh
        std::vector<curved::SurfaceFace> high_order_faces;
        // Curved geometry, or null for straight-sided meshes; it refers to this mesh
        std::shared_ptr<const curved::Geometry> curved_geometry;
        // Cells [0, n_owned()) are owned by this rank; the rest are halo cells
        uint32_t n_owned_cells = 0;
        uint32_t n_owned() const { return n_owned_cells ? n_owned_cells : n_cells; }
        // Cells [0, n_reconstructed()) need face values: owned cells and halo layer 1
        uint32_t n_reconstructed_cells = 0;
        uint32_t n_reconstructed() const { return n_reconstructed_cells ? n_reconstructed_cells : n_cells; }
        // Cells [0, n_complete()) have all their vertex neighbors: all but the outermost halo layer
        uint32_t n_complete_cells = 0;
        uint32_t n_complete() const { return n_complete_cells ? n_complete_cells : n_cells; }
        // Distributed runs: global id of every local cell and node, and the global counts (0 otherwise)
        std::vector<uint64_t> h_global_cell_id;
        uint64_t n_global_cells = 0;
        std::vector<uint64_t> h_global_node_id;
        uint64_t n_global_nodes = 0;
        // Volume of each cell and area of each face in the flow's geometry: the
        // revolved ones per radian in axisymmetric runs, else the same views as
        // cell_volume and face_area
        bool axisymmetric = false;
        Kokkos::View<rtype *> cell_measure;
        Kokkos::View<rtype *> face_measure;
        Kokkos::View<rtype *>::host_mirror_type h_cell_measure;
        Kokkos::View<rtype *>::host_mirror_type h_face_measure;
        // Axisymmetric runs: r-weighted second central moments of each cell
        // [xx, xy, yy] about its center (empty otherwise)
        Kokkos::View<rtype *[3]> cell_covariance;
        Kokkos::View<rtype *[3]>::host_mirror_type h_cell_covariance;
        Kokkos::View<rtype *[N_DIM]> node_coords;
        Kokkos::View<rtype *[N_DIM]> cell_coords;
        Kokkos::View<rtype *> cell_volume;
        Kokkos::View<rtype *> face_area;
        Kokkos::View<rtype *[N_DIM]> face_normals;
        Kokkos::View<rtype *[N_DIM]> face_coords;
        Kokkos::View<uint32_t *> nodes_of_cell;
        Kokkos::View<uint32_t *> offsets_nodes_of_cell;
        Kokkos::View<uint32_t *> faces_of_cell;
        Kokkos::View<uint32_t *> offsets_faces_of_cell;
        Kokkos::View<uint32_t *> nodes_of_face;
        Kokkos::View<uint32_t *> offsets_nodes_of_face;
        Kokkos::View<int32_t *[2]> cells_of_face;
        Kokkos::View<uint32_t *> cells_of_cell;            // Vertex neighbors (CSR)
        Kokkos::View<uint32_t *> offsets_cells_of_cell;
        // Periodic translations: row 0 of shifts is zero. Cell 1 of face f, seen
        // from cell 0, is at cell_coords(c1) + shifts(face_shift(f)); vertex
        // neighbor k of a cell is at cell_coords(cells_of_cell(k)) + shifts(cells_of_cell_shift(k))
        Kokkos::View<rtype *[N_DIM]> shifts;
        Kokkos::View<uint8_t *> face_shift;
        Kokkos::View<uint8_t *> cells_of_cell_shift;
 
        Kokkos::View<rtype *[N_DIM]>::host_mirror_type h_node_coords;
        Kokkos::View<rtype *[N_DIM]>::host_mirror_type h_cell_coords;
        Kokkos::View<rtype *>::host_mirror_type h_cell_volume;
        Kokkos::View<rtype *>::host_mirror_type h_face_area;
        Kokkos::View<rtype *[N_DIM]>::host_mirror_type h_face_normals;
        Kokkos::View<rtype *[N_DIM]>::host_mirror_type h_face_coords;
        Kokkos::View<uint32_t *>::host_mirror_type h_nodes_of_cell;
        Kokkos::View<uint32_t *>::host_mirror_type h_offsets_nodes_of_cell;
        Kokkos::View<uint32_t *>::host_mirror_type h_faces_of_cell;
        Kokkos::View<uint32_t *>::host_mirror_type h_offsets_faces_of_cell;
        Kokkos::View<uint32_t *>::host_mirror_type h_nodes_of_face;
        Kokkos::View<uint32_t *>::host_mirror_type h_offsets_nodes_of_face;
        Kokkos::View<int32_t *[2]>::host_mirror_type h_cells_of_face;
        Kokkos::View<uint32_t *>::host_mirror_type h_cells_of_cell;
        Kokkos::View<uint32_t *>::host_mirror_type h_offsets_cells_of_cell;
        Kokkos::View<rtype *[N_DIM]>::host_mirror_type h_shifts;
        Kokkos::View<uint8_t *>::host_mirror_type h_face_shift;
        Kokkos::View<uint8_t *>::host_mirror_type h_cells_of_cell_shift;
 
        // Periodicity (host): the translation of each periodic pair, the zones
        // the pairs joined, and per node its key (the lowest node id of its
        // periodic class) and lattice offset, x_n = x_key + sum_j L_j T_j
        // (empty for meshes without periodic pairs)
        std::vector<std::array<rtype, N_DIM>> periodic_translations;
        std::vector<std::string> periodic_zones;
        std::vector<uint32_t> h_node_key;
        std::vector<std::array<int8_t, 3>> h_node_lattice;
        // Lattice offset of every row of shifts
        std::vector<std::array<int8_t, 3>> shift_lattice;
 
        uint8_t shift_index(const std::array<int8_t, 3> & lattice);
    protected:
    private:
        void init_box(uint32_t nx, uint32_t ny, rtype Lx, rtype Ly, bool triangles, bool wedge);
        // Local cells in increasing global id: the first cell to visit a face becomes its cell 0
        std::vector<uint32_t> cells_by_global_id() const;
        static void orient_cells_3d(const std::vector<std::array<rtype, N_DIM>> & nodes,
                                    const std::vector<std::vector<uint32_t>> & cells,
                                    std::vector<uint32_t> & offsets, std::vector<uint32_t> & cell_nodes);
 
        void allocate_and_fill(const std::vector<std::array<rtype, N_DIM>> & nodes,
                               const std::vector<uint32_t> & cell_node_offsets,
                               const std::vector<uint32_t> & cell_nodes,
                               const std::vector<uint32_t> & cell_face_offsets,
                               const std::vector<uint32_t> & cell_faces,
                               const std::vector<uint32_t> & face_node_offsets,
                               const std::vector<uint32_t> & face_nodes,
                               const std::vector<std::array<int32_t, 2>> & face_cells,
                               const std::vector<uint8_t> & face_shifts,
                               const std::vector<uint32_t> & interior,
                               const std::map<std::string, std::vector<uint32_t>> & zone_faces);
 
        void h_neighbors_of_cell_helper(uint32_t i_cell,
                                        uint8_t n_order,
                                        std::vector<uint32_t> & neighbors) const;
 
        MeshType type = MeshType::FROM_FILE;
        std::vector<CellZone> m_cell_zones;
        std::vector<FaceZone> m_face_zones;
};
 
#endif // MESH_H