BoundaryConditions.hpp

Go to the documentation of this file.

// Class to provide description of boundary conditions
//
//
// Author: Jerome Droniou (jerome.droniou@monash.edu)
//
 
/*
*
*   This library was developed around HHO methods, although some parts of it have a more
* general purpose. If you use this code or part of it in a scientific publication, 
* please mention the following book as a reference for the underlying principles
* of HHO schemes:
*
* The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications.
* D. A. Di Pietro and J. Droniou. 2019, 516p. 
* url: https://hal.archives-ouvertes.fr/hal-02151813.
*
*/
 
 
#ifndef _BOUNDARY_CONDITIONS_HPP
#define _BOUNDARY_CONDITIONS_HPP
 
#include <iostream>
#include <string>
#include <mesh.hpp>
 
using namespace HArDCore3D;
 
// ----------------------------------------------------------------------------
//                            Class definition
// ----------------------------------------------------------------------------
 
// @addtogroup TestCases
 
struct Plane {
    Plane(const VectorRd & point, const VectorRd & normal):
      point(point),
      normal(normal)
      {
        // do nothing
      };
    
    VectorRd point; 
    VectorRd normal; 
    
    bool is_in(const VectorRd & x) const;
    inline bool is_in(const Vertex & V) const { return is_in(V.coords()); };
    template<typename MeshEntity>
    bool is_in(const MeshEntity & E) const{
      for (auto V : E.get_vertices()){
        if (! is_in(*V) ) return false;
      }
      return true;
    };
 
    bool is_on_side(const VectorRd & x, const double & side) const;
    inline bool is_on_side(const Vertex & V, const double & side) const { return is_on_side(V.coords(), side); };
    template<typename MeshEntity>
    bool is_on_side(const MeshEntity & E, const double & side) const{
      for (auto V : E.get_vertices()){
        if (! is_on_side(*V, side) ) return false;
      }
      return true;
    };
    
};
 
struct Hull {
    Hull(const std::vector<VectorRd> & nodes, const VectorRd & normal):
      nodes(nodes),
      normal(normal.normalized())
      {
        // Create center and apex
        center = VectorRd::Zero();
        for (size_t i=0; i<nodes.size(); i++){
          center += nodes[i];
        }
        center /= nodes.size();
        apex = center + normal;
      };
    
    std::vector<VectorRd> nodes; 
    VectorRd normal; 
    VectorRd center; 
    VectorRd apex; 
    
    bool is_in(const VectorRd & x) const;
    inline bool is_in(const Vertex & V) const { return is_in(V.coords()); };
    template<typename MeshEntity>
    bool is_in(const MeshEntity & E) const{
      for (auto V : E.get_vertices()){
        if (! is_in(*V) ) return false;
      }
      return true;
    };
 
};
 
 
// Various definitions of planes and convex hulls used in the BCs below
const Plane plane_x_zero(VectorRd::Zero(), VectorRd(1., 0., 0.));
const Plane plane_x_quarter(VectorRd(.25,0.,0.), VectorRd(1.,0.,0.));
const Plane plane_x_half(VectorRd(.5,0.,0.), VectorRd(1.,0.,0.));
const Hull lower_bottom_corner_x_zero({VectorRd(0.,0.,0.), VectorRd(0.,.25,0.), VectorRd(0.,.25,.25), VectorRd(0.,0.,.25)}, VectorRd(1., 0., 0.));
const Hull lower_bottom_corner_x_one({VectorRd(1.,0.,0.), VectorRd(1.,.25,0.), VectorRd(1.,.25,.25), VectorRd(1.,0.,.25)}, VectorRd(1., 0., 0.));
 
 
class BoundaryConditions {
 
public:
  BoundaryConditions(
    const std::string bc_id,  
    Mesh& mesh          
  );
 
    const std::string type(
                            const Face& face    
                          ) const ;
  const std::string type(
                         const Edge& edge    
                         ) const ;
 
  const std::string type(
                         const Vertex& vertex    
                         ) const ;
 
 
  inline const size_t n_dir_faces() const {
    return m_n_dir_faces;
  };
 
  inline const size_t n_dir_edges() const {
    return m_n_dir_edges;
  };
 
  inline const size_t n_dir_vertices() const {
    return m_n_dir_vertices;
  };
 
  inline const std::string id() { return m_bc_id; };
 
  inline const std::string name() const {
    if (m_bc_id == "D"){
      return "Dirichlet";
    }else if (m_bc_id == "N"){
      return "Neumann";
    }else if (m_bc_id == "M0"){
      return "Mixed #0 (Dirichlet at x=0, Neumann elsewhere)";
    }else if (m_bc_id == "M1"){
      return "Mixed #1 (Dirichlet at x=0 a small rectangle, Neumann elsewhere)";
    }else if (m_bc_id == "M2"){
      return "Mixed #2 (Dirichlet between planes x=.25 and x=.5, Neumann elsewhere)";
    }
    std::cout << "Unknown boundary conditions: " << m_bc_id << "\n";
    exit(1);
  };
 
  void reorder_faces(const std::string pos = "end");
 
  void reorder_edges(const std::string pos = "end");
 
  void reorder_vertices(const std::string pos = "end");
 
 
private:
  // Parameters: id of boundary condition, reference to mesh
  const std::string m_bc_id;
  Mesh& m_mesh;
 
  // Number of Dirichlet faces, edges and vertices
  size_t m_n_dir_faces;
  size_t m_n_dir_edges;
  size_t m_n_dir_vertices;
 
};
 
 
 
#endif //_BOUNDARY_CONDITION_HPP