entity.h

#pragma once

#include <iostream>

#include "reference_element.h"
#include "physical_element.h"
#include "connectivity.h"
#include "entity_data.h"

using scalar_function = std::function<double(const point_3d&)>;

class entity
{
    int     dim;
    int     tag;
    int     elemType;
    int     elemType_2D;
    int     g_order;
    int     a_order;

    size_t                              m_dof_base;
    size_t                              m_flux_base;
    size_t                              m_index_base;

    entity_ordering                     cur_elem_ordering;

    std::vector<reference_element>      reference_cells;
    std::vector<physical_element>       physical_cells;

    std::vector<reference_element>      reference_faces;
    std::vector<physical_element>       physical_faces;

    std::vector<size_t>                 faceTags;
    std::vector<size_t>                 faceNodesTags;

    using dofs_f2c_t = std::map<size_t, std::vector<size_t>>;
    void populate_lifting_matrices_planar(entity_data_cpu&, const dofs_f2c_t&) const;
    void populate_lifting_matrices_curved(entity_data_cpu&, const dofs_f2c_t&) const;
    void populate_lifting_matrices(entity_data_cpu&) const;

    void populate_differentiation_matrices_planar(entity_data_cpu&) const;
    void populate_differentiation_matrices_curved(entity_data_cpu&) const;
    void populate_differentiation_matrices(entity_data_cpu&) const;

    void populate_jacobians(entity_data_cpu&) const;
    void populate_normals(entity_data_cpu&) const;

    void recompute_tag_mappings(void);

public:
    entity() = delete;
    entity(const entity&) = delete;
    entity(entity&&) = default;
    entity& operator=(const entity&) = delete;
    entity& operator=(entity&&) = delete;

    entity(const entity_params& ep);

    /* add method to adjust offset in system!! */

    double                      measure(void) const;
    vecxd                       project(const scalar_function&) const;
    void                        project(const scalar_function&, vecxd&) const;

    size_t                      num_cells(void) const;
    size_t                      num_cell_orientations(void) const;
    std::vector<size_t>         num_cells_per_orientation(void) const;

    const physical_element&     cell(size_t) const;
    const reference_element&    cell_refelem(size_t) const;
    const reference_element&    cell_refelem(const physical_element&) const;

    size_t                      cell_global_index(size_t) const;
    size_t                      cell_global_index_by_gmsh(size_t) const;
    size_t                      cell_local_dof_offset(size_t) const;
    size_t                      cell_global_dof_offset(size_t) const;

    size_t                      num_faces(void) const;

    std::vector<size_t>         face_tags() const;

    const physical_element&     face(size_t) const;
    const reference_element&    face_refelem(size_t) const;
    const reference_element&    face_refelem(const physical_element&) const;

    size_t                      num_dofs(void) const;
    size_t                      num_fluxes(void) const;
    void                        base(size_t, size_t, size_t);

    size_t                      dof_base(void) const;
    size_t                      flux_base(void) const;
    size_t                      index_base(void) const;

    matxd                       mass_matrix(size_t) const;
    
    void                        sort_by_orientation(void);
    void                        sort_by_gmsh(void);

    int                         gmsh_tag(void) const;
    int                         gmsh_elem_type(void) const;
    entity_ordering             current_elem_ordering(void) const;

    void                        populate_entity_data(entity_data_cpu&) const;

    constexpr size_t            num_faces_per_elem(void) const { return 4; }


    std::vector<physical_element>::const_iterator  begin() const { return physical_cells.begin(); }
    std::vector<physical_element>::const_iterator    end() const { return physical_cells.end(); }

};

std::pair<std::vector<point_3d>, std::vector<quadrature_point_3d>>
integrate(const entity& e, int order);