FictitiousGrid.h

#pragma once
 
#include <sofa/version.h>
#include <sofa/core/objectmodel/BaseObject.h>
#include <sofa/defaulttype/Vec.h>
#include <sofa/defaulttype/VecTypes.h>
#include <sofa/core/topology/BaseMeshTopology.h>
 
#include <Caribou/Geometry/RectangularQuad.h>
#include <Caribou/Geometry/RectangularHexahedron.h>
#include <Caribou/Topology/Grid/Grid.h>
#include <Caribou/config.h>
 
#include <SofaCaribou/Topology/IsoSurface.h>
 
#include <memory>
#include <exception>
#include <bitset>
#include <functional>
#include <sstream>
 
#if (defined(SOFA_VERSION) && SOFA_VERSION < 201299)
namespace sofa { using Index = unsigned int; }
#endif
 
namespace SofaCaribou::topology {
 
using namespace sofa::core::objectmodel;
 
using sofa::defaulttype::Vec2Types;
using sofa::defaulttype::Vec3Types;
 
template <typename DataTypes>
class FictitiousGrid : public virtual BaseObject
{
public:
 
    SOFA_CLASS(SOFA_TEMPLATE(FictitiousGrid, DataTypes), BaseObject);
 
    static constexpr unsigned char Dimension = DataTypes::spatial_dimensions;
 
    // --------------------
    // Caribou data aliases
    // --------------------
    using Index = std::size_t;
    using Int   = INTEGER_TYPE;
    using Float = FLOATING_POINT_TYPE;
 
    // -----------------
    // Sofa data aliases
    // -----------------
    using SofaFloat = typename DataTypes::Real;
    using SofaVecInt = sofa::defaulttype::Vec<Dimension, UNSIGNED_INTEGER_TYPE>;
    using SofaVecFloat = sofa::defaulttype::Vec<Dimension, SofaFloat>;
    using Coord = typename DataTypes::Coord;
    using SofaVecCoord = sofa::helper::vector<Coord>;
    using ElementId = sofa::Index;
    using VecElementId = sofa::helper::vector<ElementId>;
    using SofaHexahedron = sofa::core::topology::BaseMeshTopology::Hexahedron;
    using SofaQuad = sofa::core::topology::BaseMeshTopology::Quad;
    using SofaTriangle = sofa::core::topology::BaseMeshTopology::Triangle;
    using SofaEdge = sofa::core::topology::BaseMeshTopology::Edge;
 
    // -----------------
    // Grid data aliases
    // -----------------
    using GridType = caribou::topology::Grid<Dimension>;
    using NodeIndex = typename GridType::NodeIndex;
    using CellIndex = typename GridType::CellIndex;
    using Dimensions = typename GridType::Dimensions;
    using Subdivisions = typename GridType::Subdivisions;
    using LocalCoordinates = typename GridType::LocalCoordinates;
    using WorldCoordinates = typename GridType::WorldCoordinates;
    using GridCoordinates = typename GridType::GridCoordinates;
    using CellSet = typename GridType::CellSet;
    using CellElement = typename GridType::Element;
 
    // -----------------
    // Structures
    // -----------------
    enum class Type : INTEGER_TYPE {
        Undefined = std::numeric_limits<INTEGER_TYPE>::lowest(),
        Inside =   (unsigned) 1 << (unsigned) 0,
        Outside =  (unsigned) 1 << (unsigned) 1,
        Boundary = (unsigned) 1 << (unsigned) 2
    };
 
    struct CellData {
        CellData(const Type & t, const Float& w, const int & r, const bool &b)
        : type(t), weight(w), region_id(r), boundary_of_region(b) {}
        Type type = Type::Undefined;
        Float weight = 0.; // 1 for full cell, 1/4 (resp. 1/8) for the first level of subdivision in 2D (resp. 3D), etc.
        int region_id = -1;
        bool boundary_of_region = false; // True if the leaf-cell makes the boundary between its region and another one
    };
 
    struct Cell {
        Cell * parent = nullptr;
        CellIndex index = 0; // Index relative to the parent cell
        std::unique_ptr<CellData> data; // Data is only stored on leaf cells
        std::unique_ptr<std::array<Cell,(unsigned) 1 << Dimension>> childs;
 
        inline bool is_leaf() const {return childs.get() == nullptr;}
    };
 
    struct Region {
        Type type = Type::Undefined;
        std::vector<Cell*> cells; // Leaf cells filling the region
    };
 
    // -------
    // Aliases
    // -------
 
    template <typename ObjectType>
    using Link = SingleLink<FictitiousGrid<DataTypes>, ObjectType, BaseLink::FLAG_STRONGLINK>;
 
    // ----------------
    // Public functions
    // ----------------
    FictitiousGrid();
 
    virtual void
    create_grid();
 
    inline UNSIGNED_INTEGER_TYPE
    number_of_cells() const {
        return p_cell_index_in_grid.size();
    }
 
    inline UNSIGNED_INTEGER_TYPE
    number_of_nodes() const {
        return p_node_index_in_grid.size();
    }
 
    inline UNSIGNED_INTEGER_TYPE
    number_of_subdivisions() const {
        return d_number_of_subdivision.getValue();
    }
 
    std::vector<Cell *>
    get_neighbors(const Cell * cell) const;
 
    std::vector<Cell *>
    get_neighbors(const Cell * cell, UNSIGNED_INTEGER_TYPE axis, INTEGER_TYPE direction) const;
 
    std::vector<std::pair<LocalCoordinates, FLOATING_POINT_TYPE>>
    get_gauss_nodes_of_cell(const CellIndex & sparse_cell_index) const;
 
    std::vector<std::pair<LocalCoordinates, FLOATING_POINT_TYPE>>
    get_gauss_nodes_of_cell(const CellIndex & sparse_cell_index, const UNSIGNED_INTEGER_TYPE level) const;
 
    inline CellElement
    get_cell_element(const CellIndex & sparse_cell_index) const {
        const auto cell_index = p_cell_index_in_grid[sparse_cell_index];
        return std::move(p_grid->cell_at(cell_index));
    }
 
    inline const SofaHexahedron &
    get_node_indices_of(const CellIndex & sparse_cell_index) const {
        return d_hexahedrons.getValue().at(sparse_cell_index);
    }
 
    inline Type
    get_type_at(const WorldCoordinates & p) const;
 
    inline std::map<FLOATING_POINT_TYPE, std::vector<CellIndex>>
    cell_volume_ratio_distribution(UNSIGNED_INTEGER_TYPE number_of_decimals=0) const;
 
    // ---------------------
    // SOFA METHOD OVERRIDES
    // ---------------------
 
    void init() override;
    void draw(const sofa::core::visual::VisualParams* vparams) override;
 
    void computeBBox(const sofa::core::ExecParams*, bool onlyVisible) override {
        if( !onlyVisible )
            return;
 
        if (Dimension == 2) {
            const Float min[3] = {
                d_min.getValue()[0], d_min.getValue()[1], -1
            };
            const Float max[3] = {
                d_max.getValue()[0], d_max.getValue()[1], +1
            };
            this->f_bbox.setValue(sofa::defaulttype::TBoundingBox<Float>(min, max));
        } else {
            this->f_bbox.setValue(sofa::defaulttype::TBoundingBox<Float>(
                d_min.getValue().array(),d_max.getValue().array()));
        }
    }
 
    std::string getTemplateName() const override {
        return templateName(this);
    }
 
    static std::string templateName(const FictitiousGrid<DataTypes>* = nullptr) {
        return DataTypes::Name();
    }
 
private:
    virtual void tag_intersected_cells_from_implicit_surface();
    virtual void tag_intersected_cells();
    virtual void tag_outside_cells();
    virtual void tag_inside_cells();
    virtual void subdivide_intersected_cells();
    virtual void create_regions_from_same_type_cells();
    virtual void create_sparse_grid();
    virtual void populate_drawing_vectors();
    virtual void validate_grid();
 
    std::array<CellElement, (unsigned) 1 << Dimension> get_subcells_elements(const CellElement & e) const;
    std::vector<Cell *> get_leaf_cells(const Cell & c) const {return std::move(get_leaf_cells(&c));}
    std::vector<Cell *> get_leaf_cells(const Cell * c) const;
    inline FLOATING_POINT_TYPE get_cell_weight(const Cell & c) const;
 
private:
    // ------------------
    // Input data members
    // ------------------
    Data<SofaVecInt> d_n;
    Data<SofaVecFloat> d_min;
    Data<SofaVecFloat> d_max;
    Data<UNSIGNED_INTEGER_TYPE> d_number_of_subdivision;
    Data<Float> d_volume_threshold;
    Link<IsoSurface<DataTypes>> d_iso_surface;
    Data<bool> d_draw_boundary_cells;
    Data<bool> d_draw_outside_cells;
    Data<bool> d_draw_inside_cells;
    Data<double> d_draw_scale;
 
    Data< SofaVecCoord > d_surface_positions;
 
    Data<sofa::helper::vector<SofaEdge> > d_surface_edges;
 
    Data<sofa::helper::vector<SofaTriangle> > d_surface_triangles;
 
    // -------------------
    // Output data members
    // -------------------
    Data< SofaVecCoord > d_positions;
 
    Data < sofa::helper::vector<SofaQuad> > d_quads;
 
    Data < sofa::helper::vector<SofaHexahedron> > d_hexahedrons;
 
 
    // ---------------
    // Private members
    // ---------------
    std::unique_ptr<GridType> p_grid;
 
    std::vector<Type> p_cells_types;
 
    std::vector<std::vector<Index>> p_triangles_of_cell;
 
    std::vector<Cell> p_cells;
 
    std::vector<Region> p_regions;
 
    std::vector<INTEGER_TYPE> p_node_index_in_sparse_grid;
 
    std::vector<UNSIGNED_INTEGER_TYPE> p_node_index_in_grid;
 
    std::vector<INTEGER_TYPE> p_cell_index_in_sparse_grid;
 
    std::vector<UNSIGNED_INTEGER_TYPE> p_cell_index_in_grid;
 
    std::vector<sofa::defaulttype::Vector3> p_drawing_nodes_vector;
 
    std::vector<sofa::defaulttype::Vector3> p_drawing_edges_vector;
 
    std::vector<std::vector<sofa::defaulttype::Vector3>> p_drawing_subdivided_edges_vector;
 
    std::vector<std::vector<sofa::defaulttype::Vector3>> p_drawing_cells_vector;
 
    // ----------------------
    // Private static members
    // ----------------------
    static const GridCoordinates subcell_coordinates[(unsigned) 1 << Dimension];
 
};
 
template<> const FictitiousGrid<Vec2Types>::GridCoordinates FictitiousGrid<Vec2Types>::subcell_coordinates[4];
template<>  const FictitiousGrid<Vec3Types>::GridCoordinates FictitiousGrid<Vec3Types>::subcell_coordinates[8];
 
template<> void FictitiousGrid<Vec2Types>::tag_intersected_cells ();
template<> void FictitiousGrid<Vec3Types>::tag_intersected_cells ();
 
template<> void FictitiousGrid<Vec2Types>::subdivide_intersected_cells ();
template<> void FictitiousGrid<Vec3Types>::subdivide_intersected_cells ();
 
extern template class FictitiousGrid<Vec2Types>;
extern template class FictitiousGrid<Vec3Types>;
 
} // namespace SofaCaribou::topology