mesh.hpp

Go to the documentation of this file.

/* 
 * Copyright (c) 2023 Marien Hanot <marien-lorenzo.hanot@umontpellier.fr>
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *  
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#ifndef MESH_HPP_INCLUDED
#define MESH_HPP_INCLUDED
 
#include "dcell.hpp"
#include "loader.hpp" // Needed to destruct the loader
 
#include <memory>
 
namespace Manicore {
 
  template<size_t> class Mesh_builder;
 
 
 
  template<size_t dimension>
    class Mesh 
    {
      public:
        size_t n_cells(size_t d ) const {return _cells_graph[d].size();}
        
        // Interface to the graph structure
        std::vector<size_t> const & get_map_ids(size_t d ,
                                                size_t i_cell ) const 
          {return _cells_graph[d][i_cell].get_map_ids();}
 
        std::vector<size_t> const & get_boundary(size_t d_boundary , 
                                                 size_t d ,
                                                 size_t i_cell ) const 
          {return _cells_graph[d][i_cell].get_boundary(d_boundary);}
 
 
        std::vector<size_t> const & get_relative_map(size_t d_boundary , 
                                                     size_t d ,
                                                     size_t i_cell ) const
          {return _cells_graph[d][i_cell].get_relative_map(d_boundary);}
 
 
        size_t global_to_local(size_t d_boundary ,
                               size_t d ,
                               size_t i_b ,
                               size_t i_cell ) const 
          {return _cells_graph[d][i_cell].global_to_local(d_boundary,i_b);}
 
 
        int get_boundary_orientation(size_t d , 
                                     size_t i_cell , 
                                     size_t j_bd ) const; 
 
 
        Eigen::Vector3d get_3D_embedding (size_t m_id ,
                                          Eigen::Vector<double,dimension> const &x ) const 
        {return _maps[m_id]->I(x);}
 
        Eigen::Matrix<double,3,dimension> get_3D_pushforward (size_t m_id ,
                                          Eigen::Vector<double,dimension> const &x ) const 
        {return _maps_pullback[m_id]->DI(x);}
 
        // Access to the geometric structure
 
        template<size_t d> dCell_map<dimension,d> const & get_cell_map(size_t i) const;
 
        // Access to the metric structure
        Eigen::Matrix<double,dimension,dimension> metric (size_t m_id ,
                                                          Eigen::Vector<double,dimension> const &x ) const 
        {return _metric_maps[m_id]->metric(x);} // On the tangent space
 
        Eigen::Matrix<double,dimension,dimension> metric_inv (size_t m_id ,
                                                              Eigen::Vector<double,dimension> const &x ) const 
        {return _metric_maps[m_id]->metric_inv(x);} // On the cotangent space
 
        double volume_form (size_t m_id ,
                            Eigen::Vector<double,dimension> const &x ) const 
        {return _metric_maps[m_id]->volume(x);}
 
        int orientationTopCell(size_t i_cell) const
        {
          return _metric_maps[get_map_ids(dimension,i_cell)[0]]->orientation;
        }
 
 
        template<size_t k,int d> Eigen::Matrix<double,Dimension::ExtDim(d-k,d),Dimension::ExtDim(k,d)>
          getHodge(size_t i_cell , 
                   Eigen::Vector<double,d> const &x ) const; // Compute the Hodge star operator on the i_cell-th d-cell at x
 
      private:
        static constexpr size_t _max_dim = 3;
        Mesh() requires(dimension <= _max_dim) {;}
        friend class Mesh_builder<dimension>;
        typename dCell_graph<dimension>::Collection _cells_graph;
        std::unique_ptr<Maps_loader<dimension>> _loader_ref; // Tie the life of the shared library to the mesh (must be destroyed AFTER every pointer to any mapping 
        std::vector<std::unique_ptr<ParametrizedMap<3,dimension>>> _maps;
        std::vector<std::unique_ptr<ParametrizedDerivedMap<3,dimension>>> _maps_pullback;
        std::vector<std::unique_ptr<ParametrizedMetricMap<dimension>>> _metric_maps;
        // TODO use the generic solution of the class PEC
        std::vector<dCell_map<dimension,0>> _geo0;
        std::vector<dCell_map<dimension,1>> _geo1;
        std::vector<dCell_map<dimension,2>> _geo2;
        std::vector<dCell_map<dimension,3>> _geo3;
    };
 
  // ---------------------------------------------------------------------------------------------------------
  // Implementation
  // ---------------------------------------------------------------------------------------------------------
    template<size_t dimension>
    int Mesh<dimension>::get_boundary_orientation(size_t d, size_t i_cell, size_t j_bd) const
    {
      size_t i_bd_abs = get_boundary(d-1,d,i_cell)[j_bd];
      size_t bd_rel_map = get_relative_map(d-1,d,i_cell)[j_bd];
      auto get_orientation = [&]<size_t _d>(auto&& get_orientation)
      {
        if constexpr(_d == dimension+1) {
          assert(false && "Dimension too high or low");
          return 0;
        } else {
          if (_d == d) {
            auto const & E = get_cell_map<_d-1>(i_bd_abs);
            auto x = Geometry::middleSimplex<_d-1>(E.get_reference_elem()[0]);
            auto pM = E.evaluate_DI(bd_rel_map,x);
    if (d == dimension) {
            return get_cell_map<_d>(i_cell).get_orientation(0,E.evaluate_I(bd_rel_map,x),pM)*orientationTopCell(i_cell);
    } else {
            return get_cell_map<_d>(i_cell).get_orientation(0,E.evaluate_I(bd_rel_map,x),pM);
    }
          } else {
            return get_orientation.template operator()<_d+1>(get_orientation);
          }
        }
      };
      return get_orientation.template operator()<2>(get_orientation); // Start at dim(T) = 2, since the boundary must be at least of dimension 1
    }
 
    // Doxygen wrongly assumes this is an overload, prevent it from duplicating the entry
    template<size_t dimension>
    template<size_t d> 
    dCell_map<dimension,d> const & Mesh<dimension>::get_cell_map(size_t i) const
    {
      static_assert(d <= _max_dim);
      if constexpr (d == 0) return _geo0[i];
      else if constexpr (d == 1) return _geo1[i];
      else if constexpr (d == 2) return _geo2[i];
      else return _geo3[i];
    }
 
    template<size_t dimension>
    template<size_t k,int d> Eigen::Matrix<double,Dimension::ExtDim(d-k,d),Dimension::ExtDim(k,d)>
    Mesh<dimension>::getHodge(size_t i_cell, Eigen::Vector<double,d> const &x) const 
    {
      auto const & F = get_cell_map<d>(i_cell);
      auto const Ix = F.evaluate_I(0,x);
      Eigen::Matrix<double,d,d> invGF;
      double sqrtdetG;
      if constexpr(d == dimension) { // TODO check if this is actually faster
        auto const DJ = F.evaluate_DJ(0,Ix);
        invGF = DJ*metric_inv(get_map_ids(d,i_cell)[0],Ix)*DJ.transpose();
        sqrtdetG = volume_form(get_map_ids(d,i_cell)[0],Ix) * std::abs(F.evaluate_DI(0,x).determinant())*orientationTopCell(i_cell);
      } else {
        auto const DI = F.evaluate_DI(0,x);
        Eigen::Matrix<double,d,d> pullbackMetric = DI.transpose()*metric(get_map_ids(d,i_cell)[0],Ix)*DI;
        invGF = (pullbackMetric).inverse();
        sqrtdetG = std::sqrt(pullbackMetric.determinant());
      }
      auto const partialDet = Compute_pullback<k,d,d>::compute(invGF);
      auto const &complBasis = ComplBasis<k,d>::compute();
      return sqrtdetG*complBasis*partialDet;
      // ComplBasis does not duplicate indices, hence there is no 1/(d-k)! term
    }
 
} // Namespace
 
#endif