dcell.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 DCELL_HPP_INCLUDED
#define DCELL_HPP_INCLUDED
 
#include "exterior_algebra.hpp"
#include "geometry.hpp"
 
#include <vector>
#include <set>
#include <memory>
 
namespace Manicore {
 
 
  template<size_t dimension>
  class dCell_graph {
    public:
      typedef std::array<std::vector<dCell_graph>,dimension+1> Collection;
 
      dCell_graph(size_t id , 
                  std::vector<size_t> const & map_ids );
 
      dCell_graph(Collection const & ,
                  size_t d , 
                  size_t id , 
                  std::vector<size_t> const & map_ids , 
                  std::vector<size_t> const & bd );
 
      size_t get_id() const {return _bd_ids[_d][0];}
      std::vector<size_t> const & get_map_ids() const {return _map_ids;}
      std::vector<size_t> const & get_boundary(size_t dim ) const {return _bd_ids[dim];}
      std::vector<size_t> const & get_relative_map(size_t dim ) const {return _bd_rel_maps[dim];}
      size_t global_to_local(size_t d_boundary ,size_t i_b ) const;
    private:
      size_t _d;
      std::vector<size_t> _map_ids;
      std::array<std::vector<size_t>,dimension+1> _bd_ids; 
      std::array<std::vector<size_t>,dimension> _bd_rel_maps; // relative ids of the corresponding map in the boundary
      size_t _find_map(size_t id) const;
  };
  
 
  template<size_t dimension,size_t d>
  class dCell_map {
    public:
      static constexpr size_t cell_dim = d;
 
 
      dCell_map(Eigen::Vector<double,dimension> const &center_mass ,
                Eigen::Matrix<double,dimension,d> const &flat_map ,
                double diam ,
                std::vector<Simplex<d>> const &triangulation );
      dCell_map(std::vector<std::unique_ptr<ParametrizedMap<dimension,d>>> & maps, 
                std::vector<std::unique_ptr<ParametrizedDerivedMap<dimension,d>>> & pullback_maps ,
                std::vector<Simplex<d>> const &triangulation );
 
 
      bool is_flat() const {return _is_flat;}
      // Geometry
      std::vector<Simplex<d>> const & get_reference_elem() const {return _triangulation;}
 
      Eigen::Vector<double,dimension> evaluate_I(size_t rel_map_id , 
                                                 Eigen::Vector<double,d> const &x ) const;
 
      Eigen::Vector<double,d> evaluate_J(size_t rel_map_id , 
                                         Eigen::Vector<double,dimension> const &x ) const;
      // Metric
      double evaluate_poly_on_ref(Eigen::Vector<double,d> const &x ,
                                  size_t i_pbasis , 
                                  int r ) const;
 
      double evaluate_poly_pullback(size_t rel_map_id , 
                                    Eigen::Vector<double,dimension> const &x ,
                                    size_t i_pbasis ,
                                    int r ) const;
 
      Eigen::Matrix<double,dimension,d> evaluate_DI(size_t rel_map_id , 
                                                    Eigen::Vector<double,d> const &x ) const;
 
      Eigen::Matrix<double,d,dimension> evaluate_DJ(size_t rel_map_id ,
                                                    Eigen::Vector<double,dimension> const &x ) const;
 
 
      template<size_t l> Eigen::Matrix<double,Dimension::ExtDim(l,d),Dimension::ExtDim(l,dimension)> 
        evaluate_DI_p(size_t rel_map_id , 
                      Eigen::Vector<double,d> const &x ) const; // Return the pullback by I of the basis of l-forms
 
 
      template<size_t l> Eigen::Matrix<double,Dimension::ExtDim(l,dimension),Dimension::ExtDim(l,d)> 
        evaluate_DJ_p(size_t rel_map_id , 
                      Eigen::Vector<double,dimension> const &x ) const; // Return the pullback by J of the basis of l-forms
 
      // Orientation
 
      int get_orientation(size_t rel_map_id , 
                             Eigen::Vector<double,dimension> const &x ,
                             Eigen::Matrix<double,dimension,d-1> const &pM 
                             ) const requires(d>1);
 
    private:
      Eigen::Vector<double,dimension> _center_mass; // Used when _is_flat is true 
      Eigen::Matrix<double,dimension,d> _flat_map; // Used when _is_flat is true and d<dimension
      double _diam; // Used when _is_flat is true
      std::vector<std::unique_ptr<ParametrizedMap<dimension,d>>> _maps; // Used when _is_flat is false
      std::vector<std::unique_ptr<ParametrizedDerivedMap<dimension,d>>> _pullback_maps; // Used when _is_flat is false
      const bool _is_flat; // Element and all its boundaries are flat, it must also correspond to a single map for d>0
 
      std::vector<Simplex<d>> _triangulation; // Used to compute integrals
  };
 
 
  template<size_t dimension>
  class dCell_map<dimension,0> { // Specialize for vertices
    public:
      dCell_map(std::vector<Eigen::Vector<double,dimension>> const &coords ) : _locs(coords) {};
 
      Eigen::Vector<double,dimension> const & coord(size_t map_rel_id ) const {return _locs[map_rel_id];}
 
    private:
      std::vector<Eigen::Vector<double,dimension>> _locs;
  };
 
  // ---------------------------------------------------------------------------------------------------------
  // Implementation
  // ---------------------------------------------------------------------------------------------------------
  // ---------------------------------------------------------------------------------------------------------
  // Graph
 
  template<size_t dimension>
  dCell_graph<dimension>::dCell_graph(size_t id, std::vector<size_t> const & map_ids)
  : _d(0), _map_ids(map_ids) {
    _bd_ids[_d] = std::vector<size_t>{id};
  }
 
  template<size_t dimension>
  dCell_graph<dimension>::dCell_graph(Collection const &coll,size_t d, size_t id, std::vector<size_t> const & map_ids, std::vector<size_t> const & bd) 
  : _d(d), _map_ids(map_ids) {
    assert(d > 0);
    // Insert map of faces
    _bd_ids[d] = std::vector<size_t>{id};
    _bd_ids[d-1] = bd;
    for (size_t i = 0; i < bd.size(); ++i) {
      _bd_rel_maps[d-1].push_back(coll[d-1][bd[i]]._find_map(_map_ids[0]));
    }
    // Iterate lower dimensions
    for (size_t dim = d-1; dim > 0; --dim) {
      std::set<size_t> clist;
      for (size_t i = 0; i < _bd_ids[dim].size(); ++i) {
        std::vector<size_t> const & cb = coll[dim][_bd_ids[dim][i]].get_boundary(dim-1);
        for (size_t j = 0; j < cb.size(); ++j) {
          clist.insert(cb[j]);
        }
      }
      for (auto const &p : clist) {
        _bd_ids[dim-1].push_back(p);
        _bd_rel_maps[dim-1].push_back(coll[dim-1][p]._find_map(_map_ids[0]));
      }
    }
  }
  
      // This is relative to the principal map of the object (to _map_ids[0])
  template<size_t dimension>
  size_t dCell_graph<dimension>::_find_map(size_t id) const {
    for (size_t i = 0; i < _map_ids.size(); ++i) {
      if (_map_ids[i] == id) return i;
    }
    return _map_ids.size();
  }
 
  template<size_t dimension>
  size_t dCell_graph<dimension>::global_to_local(size_t d_boundary,size_t i_b) const
  {
    for (size_t i = 0; i < _bd_ids[d_boundary].size(); ++ i) {
      if (_bd_ids[d_boundary][i] == i_b) return i;
    }
    assert(false);
    return _bd_ids[d_boundary].size();
  }
 
  // ---------------------------------------------------------------------------------------------------------
  // Map
 
  template<size_t dimension,size_t d>
  dCell_map<dimension,d>::dCell_map(Eigen::Vector<double,dimension> const &center_mass,
                                    Eigen::Matrix<double,dimension,d> const &flat_map,
                                    double diam, std::vector<Simplex<d>> const &triangulation)
    : _center_mass(center_mass), _flat_map(flat_map), _diam(diam), _is_flat(true), _triangulation(triangulation) 
  {}
 
  template<size_t dimension,size_t d>
  dCell_map<dimension,d>::dCell_map(std::vector<std::unique_ptr<ParametrizedMap<dimension,d>>> &maps,
                std::vector<std::unique_ptr<ParametrizedDerivedMap<dimension,d>>> & pullback_maps,
                std::vector<Simplex<d>> const &triangulation)
  : _maps(std::move(maps)), _pullback_maps(std::move(pullback_maps)), _is_flat(false), _triangulation(triangulation)
  {}
 
  template<size_t dimension,size_t d>
  Eigen::Vector<double,dimension> dCell_map<dimension,d>::evaluate_I(size_t rel_map_id, Eigen::Vector<double,d> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return x*_diam  + _center_mass;
      } else {
        return _flat_map*x*_diam + _center_mass;
      }
    } else {
      assert(rel_map_id < _maps.size());
      return _maps[rel_map_id]->I(x);
    }
  }
 
  template<size_t dimension,size_t d>
  Eigen::Vector<double,d> dCell_map<dimension,d>::evaluate_J(size_t rel_map_id, Eigen::Vector<double,dimension> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return (x - _center_mass)/_diam;
      } else {
        return _flat_map.transpose()*(x - _center_mass)/_diam;
      }
    } else {
      assert(rel_map_id < _maps.size());
      return _maps[rel_map_id]->J(x);
    }
  }
 
  // Evaluations
  template<size_t dimension,size_t d>
  double dCell_map<dimension,d>::evaluate_poly_on_ref(Eigen::Vector<double,d> const &x,size_t i_pbasis,int r) const 
  {
    assert((i_pbasis < Dimension::PolyDim(r,d)) && "Index out of range");
    auto const & power = Monomial_powers<d>::complete(r)[i_pbasis];
    double rv = std::pow(x(0),power[0]);
    for (size_t i = 1; i < d; ++i) {
      rv *= std::pow(x(i),power[i]);
    }
    return rv;
  }
  
  template<size_t dimension,size_t d>
  double dCell_map<dimension,d>::evaluate_poly_pullback(size_t rel_map_id, Eigen::Vector<double,dimension> const &x,size_t i_pbasis,int r) const 
  {
    return evaluate_poly_on_ref(evaluate_J(rel_map_id,x),i_pbasis,r);
  }
 
  template<size_t dimension,size_t d>
  Eigen::Matrix<double,dimension,d> dCell_map<dimension,d>::evaluate_DI(size_t rel_map_id, Eigen::Vector<double,d> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return Eigen::Matrix<double,d,d>::Identity()*_diam;
      } else {
        return _flat_map*_diam;
      }
    } else {
      return _pullback_maps[rel_map_id]->DI(x);
    }
  }
 
  template<size_t dimension,size_t d>
  Eigen::Matrix<double,d,dimension> dCell_map<dimension,d>::evaluate_DJ(size_t rel_map_id, Eigen::Vector<double,dimension> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return Eigen::Matrix<double,d,d>::Identity()/_diam;
      } else {
        return _flat_map/_diam;
      }
    } else {
      return _pullback_maps[rel_map_id]->DJ(x);
    }
  }
 
  template<size_t dimension,size_t d>
    template<size_t l>
    Eigen::Matrix<double,Dimension::ExtDim(l,d),Dimension::ExtDim(l,dimension)> dCell_map<dimension,d>::evaluate_DI_p(size_t rel_map_id, Eigen::Vector<double,d> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return Eigen::Matrix<double,Dimension::ExtDim(l,d),Dimension::ExtDim(l,d)>::Identity()*std::pow(_diam,l);
      } else {
        return Compute_pullback<l,d,dimension>::compute(_flat_map*_diam);
      }
    } else {
      return Compute_pullback<l,d,dimension>::compute(_pullback_maps[rel_map_id]->DI(x));
    }
  }
 
  template<size_t dimension,size_t d>
    template<size_t l>
    Eigen::Matrix<double,Dimension::ExtDim(l,dimension),Dimension::ExtDim(l,d)> dCell_map<dimension,d>::evaluate_DJ_p(size_t rel_map_id, Eigen::Vector<double,dimension> const &x) const
  {
    if (_is_flat) {
      if constexpr(dimension == d) {
        return Eigen::Matrix<double,Dimension::ExtDim(l,d),Dimension::ExtDim(l,d)>::Identity()*std::pow(1/_diam,l);
      } else {
        return Compute_pullback<l,dimension,d>::compute(_flat_map.transpose()/_diam);
      }
    } else {
      return Compute_pullback<l,dimension,d>::compute(_pullback_maps[rel_map_id]->DJ(x));
    }
  }
 
  template<size_t dimension,size_t d>
  int dCell_map<dimension,d>::get_orientation(size_t rel_map_id, 
      Eigen::Vector<double,dimension> const &x,
      Eigen::Matrix<double,dimension,d-1> const &pM) const requires(d>1) 
  {
    auto const Jx = evaluate_J(rel_map_id,x);
    for (auto const &S : _triangulation) {
      if (Geometry::inside<d>(Jx,S)) {
        auto mid = Geometry::middleSimplex<d>(S);
        auto DI = evaluate_DI(rel_map_id,Jx);
        Eigen::Matrix<double,dimension,d> cmpM;
        cmpM.rightCols(d-1) = pM;
        cmpM.leftCols(1) = DI*(Jx - mid).normalized();
        double tmp = (evaluate_DJ(rel_map_id,x)*cmpM).determinant(); 
        assert(std::abs(tmp) > 1e-10 && "Volume form is 0");
        return (tmp > 0.)? 1 : -1;
      }
    }
    assert(false && "Point outside of element");
    return 0.;
  }
 
} // Namespace
 
#endif