DomainDiscretization_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_DOMAINS_DOMAINDISCRETIZATION_FWD_HPP_
#define MEDUSA_BITS_DOMAINS_DOMAINDISCRETIZATION_FWD_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/types/Range_fwd.hpp>
#include <medusa/bits/utils/memutils.hpp>
#include <medusa/bits/domains/DomainShape_fwd.hpp>
#include <medusa/bits/operators/shape_flags.hpp>
#include <iosfwd>
#include <string>
 
namespace mm {
 
template <typename vec_t, typename OpFamilies>
class RaggedShapeStorage;
 
template <class vec_t>
class DomainDiscretization {
  public:
    typedef vec_t vector_t;  
    typedef typename vec_t::Scalar scalar_t;   
    enum {  dim = vec_t::dim };
 
  protected:
    Range<vec_t> positions_;  
 
    Range<Range<int>> support_;
    Range<int> types_;
 
    Range<int> boundary_map_;
 
    Range<vec_t> normals_;
 
    deep_copy_unique_ptr<DomainShape<vec_t>> shape_;
 
  public:
    explicit DomainDiscretization(const DomainShape<vec_t>& shape);
 
    template <typename hdf5_file>
    static DomainDiscretization<vec_t> load(hdf5_file& file, const std::string& name);
 
    const Range<vec_t>& positions() const { return positions_; }
    const vec_t& pos(int i) const { return positions_[i]; }
    vec_t& pos(int i) { return positions_[i]; }
    scalar_t pos(int i, int j) const { return positions_[i][j]; }
    const Range<Range<int>>& supports() const { return support_; }
    const Range<int>& support(int i) const { return support_[i]; }
    Range<int>& support(int i) { return support_[i]; }
    int support(int i, int j) const { return support_[i][j]; }
    Range<vec_t> supportNodes(int i) const { return positions_[support_[i]]; }
    vec_t supportNode(int i, int j) const { return positions_[support_[i][j]]; }
    scalar_t dr(int i) const { return (positions_[i] - positions_[support_[i][1]]).norm(); }
    int supportSize(int i) const { return support_[i].size(); }
    Range<int> supportSizes() const;
    const Range<int>& types() const { return types_; }
    Range<int>& types() { return types_; }
    int type(int i) const { return types_[i]; }
    int& type(int i) { return types_[i]; }
    const DomainShape<vec_t>& shape() const { return *shape_; }
 
    const Range<int>& bmap() const { return boundary_map_; }
    int bmap(int node) const { return boundary_map_[node]; }
    const Range<vec_t>& normals() const { return normals_; }
    const vec_t& normal(int i) const;
    vec_t& normal(int i);
 
    indexes_t boundary() const { return types_ < 0; }
    indexes_t interior() const { return types_ > 0; }
    indexes_t all() const { return types_ != 0; }
 
    int size() const { return positions_.size(); }
 
    int addInternalNode(const vec_t& point, int type);
 
    int addBoundaryNode(const vec_t& point, int type, const vec_t& normal);
 
  private:
    int addNode(const vec_t& point, int type);
    void chopOff(const DomainDiscretization<vec_t>& d, const Range<int>& relevant);
 
  public:
    indexes_t addNodes(const DomainDiscretization<vec_t>& d);
 
    void changeToBoundary(int i, const vec_t& point, int type, const vec_t& normal);
 
    void changeToInterior(int i, const vec_t& point, int type);
 
    Range<int> addGhostNodes(scalar_t h, int type = 0);
    template <typename func_t>
    Range<int> addGhostNodes(func_t h, int type = 0);
    Range<int> addGhostNodes(scalar_t h, int type, const indexes_t& indexes);
    template <typename func_t>
    Range<int> addGhostNodes(func_t h, int type, const indexes_t& indexes);
 
    void shuffleNodes(const indexes_t& permutation);
 
    template <typename compare_t = std::less<std::pair<vec_t, int>>>
    Range<int> reorderNodes(const compare_t& cmp = compare_t());
 
    DomainDiscretization& add(const DomainDiscretization& d);
    DomainDiscretization& operator+=(const DomainDiscretization& d) { return add(d); }
 
    DomainDiscretization& subtract(const DomainDiscretization& d);
    DomainDiscretization& operator-=(const DomainDiscretization& d) { return subtract(d); }
 
    DomainDiscretization& translate(const vec_t& a);
 
    DomainDiscretization& rotate(const Eigen::Matrix<scalar_t, dim, dim>& Q);
 
    DomainDiscretization& rotate(scalar_t angle);
 
    void assert_is_valid() const;
 
    void clear();
 
    void removeNodes(const Range<int>& to_remove);
 
    void removeInternalNodes() { removeNodes(types_ > 0); }
 
    void removeBoundaryNodes() { removeNodes(types_ < 0); }
 
    bool contains(const vec_t& point) const;
 
    template <typename search_structure_t>
    bool discreteContains(const vec_t& point, const search_structure_t& search) const;
 
    template <typename search_structure_t>
    bool contains(const vec_t& point, const search_structure_t& search) const;
 
    template <typename search_structure_t>
    void makeDiscreteContainsStructure(search_structure_t& search) const;
 
    // TODO(jureslak): optimize order of nodes for matrix fill  -- sort by x coordinate or sth...
    // TODO(jureslak): symmetrize support?
 
    // ENGINE PLUGINS
    #define DOMAIN_PLUGIN(Name) \
        template<typename callable_t, typename... Args> \
        void Name(callable_t& callable, Args&&... args) { \
            return callable(*this, std::forward<Args>(args)...); \
        }
 
    #define DOMAIN_PLUGIN_CONST(Name) \
        template<typename callable_t, typename... Args> \
        void Name(const callable_t& callable, Args&&... args) { \
            callable(*this, std::forward<Args>(args)...); \
        }
 
    DOMAIN_PLUGIN(findSupport)
    DOMAIN_PLUGIN_CONST(findSupport)
    DOMAIN_PLUGIN(fill)
    DOMAIN_PLUGIN_CONST(fill)
    DOMAIN_PLUGIN(relax)
    DOMAIN_PLUGIN_CONST(relax)
 
 
    template <sh::shape_flags mask = sh::all, typename approx_t>
    RaggedShapeStorage<vec_t, typename sh::operator_tuple<mask, vec_t::dim>::type> computeShapes(
            approx_t approx, const indexes_t& indexes = {}) const;
    template <typename OperatorTuple, typename approx_t>
    RaggedShapeStorage<vec_t, OperatorTuple> computeShapes(
            approx_t approx, const indexes_t& indexes = {}) const;
    template <typename OperatorTuple, typename approx_t>
    RaggedShapeStorage<vec_t, OperatorTuple> computeShapes(OperatorTuple operators,
            approx_t approx, const indexes_t& indexes = {}) const;
    template <class V>
    friend std::ostream& operator<<(std::ostream& os, const DomainDiscretization<V>& d);
  private:
    std::ostream& output_report(std::ostream& os = std::cout) const;
};
 
template <class vec_t>
std::ostream& operator<<(std::ostream& os, const DomainDiscretization<vec_t>& d) {
    os << "Domain:\n";
    return d.output_report(os);
}
 
}  // namespace mm
 
#endif  // MEDUSA_BITS_DOMAINS_DOMAINDISCRETIZATION_FWD_HPP_