curved.h

Go to the documentation of this file.

 
#ifndef CURVED_H
#define CURVED_H
 
#include <array>
#include <cstdint>
#include <map>
#include <string>
#include <utility>
#include <vector>
 
#include <toml.hpp>
 
#include "common_typedef.h"
 
class Mesh;
struct MeshBlock;
 
namespace curved {
 
using Vec3 = std::array<double, 3>;
 
struct Shape {
    enum class Kind : uint8_t { CIRCLE, SPHERE, CYLINDER };
    Kind kind = Kind::CIRCLE;
    std::string zone;
    Vec3 center = {0.0, 0.0, 0.0};
    Vec3 axis = {0.0, 0.0, 1.0};  // cylinder: unit vector
    double radius = 1.0;
 
    void project(const Vec3 & x, Vec3 & p, double * J) const;
};
 
struct SurfaceFace {
    std::vector<uint64_t> corners;
    int32_t shape = -1;
    std::vector<Vec3> nodes;
};
 
struct Surfaces {
    std::vector<Shape> shapes;
    std::vector<SurfaceFace> faces;
 
    bool empty() const { return faces.empty(); }
};
 
std::vector<Shape> read_shapes(const toml::value & input);
 
bool enabled(const toml::value & input);
 
uint32_t element_order(uint32_t n_corners, uint32_t n_nodes);
 
Surfaces surfaces_of_mesh(const Mesh & mesh, const toml::value & input);
 
Surfaces surfaces_of_block(const MeshBlock & block, const toml::value & input);
 
class Geometry {
    public:
        Geometry(const Mesh & mesh, const Surfaces & surfaces);
 
        bool face_is_curved(uint32_t f) const { return face_index[f] >= 0; }
        bool cell_is_curved(uint32_t c) const { return curved_cell[c] != 0; }
        uint32_t n_curved_faces() const { return uint32_t(faces.size()); }
 
        double max_shape_offset() const;
 
        void eval(uint32_t f, const double * xi, Vec3 & x, Vec3 * dx) const;
 
        Vec3 area_vector(uint32_t f) const;
 
        void face_rule(uint32_t f, const std::vector<std::array<double, 2>> & ref, const std::vector<double> & ref_w,
                       std::vector<Vec3> & x, std::vector<Vec3> & n, std::vector<double> & w) const;
 
        void cell_rule(uint32_t c, int degree, std::vector<Vec3> & x, std::vector<double> & w) const;
 
    private:
        struct Edge {
            int32_t shape = -1;
            std::vector<Vec3> nodes;  // Lagrange nodes, equispaced from the lower global node id to the higher
        };
        struct Face {
            int32_t shape = -1;                        // analytic surface face: its own projection
            std::array<int32_t, 4> edge = {-1, -1, -1, -1};  // local edge k (node k to k + 1), or -1 if straight
            std::array<uint8_t, 4> flip = {0, 0, 0, 0};      // edge k runs from its higher global node id
            std::vector<Vec3> bubble;                  // interior corrections of high-order elements
        };
 
        void edge_eval(const Face & face, size_t k, const Vec3 & a, const Vec3 & b, double t, Vec3 & x, Vec3 & dx) const;
        void face_eval(const Face & face, const std::vector<Vec3> & v, const double * xi, Vec3 & x, Vec3 * dx) const;
        std::vector<Vec3> face_nodes(uint32_t f) const;
 
        const Mesh & mesh;
        std::vector<Shape> shapes;
        std::vector<Edge> edges;
        std::vector<Face> faces;
        std::vector<int32_t> face_index;  // per local face: index into faces, or -1
        std::vector<uint8_t> curved_cell;
};
 
}  // namespace curved
 
#endif  // CURVED_H