BChandlers.hpp

Go to the documentation of this file.

// Free functions to help handle BC and DOFs/indices
//
// 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 BCHANDLERS_HPP
#define BCHANDLERS_HPP
 
#include <BoundaryConditions/BoundaryConditions.hpp>
#include <globaldofspace.hpp>
#include <mesh.hpp>
 
using namespace HArDCore2D;
 
// ----------------------------------------------------------------------------
//                            Class definition
// ----------------------------------------------------------------------------
 
void setBCLabels(const BoundaryConditions &BC, GlobalDOFSpace &globaldofspace) {
  
  for (Vertex *v : globaldofspace.mesh().get_vertices()){
    int label = 0;
    if (BC.type(*v)=="dir"){
      label = 1;
    }else if (BC.type(*v)=="neu"){
      label = -1;
    }
    size_t offset_v = globaldofspace.globalOffset(*v);
    for (size_t i=offset_v; i < offset_v + globaldofspace.numLocalDofsVertex(); i++){
      globaldofspace.setLabelDOF(i, label);
    }          
  }
 
  for (Edge *e : globaldofspace.mesh().get_edges()){
    int label = 0;
    if (BC.type(*e)=="dir"){
      label = 1;
    }else if (BC.type(*e)=="neu"){
      label = -1;
    }
    size_t offset_e = globaldofspace.globalOffset(*e);
    for (size_t i=offset_e; i < offset_e + globaldofspace.numLocalDofsEdge(); i++){
      globaldofspace.setLabelDOF(i, label);
    }          
  }
 
}
 
 
 
int offsetIndex(const std::vector<size_t> &c, const int &i){
  assert( c.size() >= 2 && c.size()%2 == 0);
  
  int val = -1;
  int idx = -1;
  int offset = 0;
  // Look for index in c such that c[idx]>i and add offsets (every other pair of indices)
  do {
    idx++;
    if (idx%2 == 1){
      // c.size() even so no risk of having idx=c.size() here
      offset += c[idx]-c[idx-1];
    } 
  }
  while (idx<int(c.size()) && int(c[idx])<=i);
  
  if (idx%2 == 0){
    val = i - offset;
  }
  
  return val;
  
}
 
Eigen::ArrayXi create_mapDOF(const std::vector<size_t> &c, const size_t N){
  assert( c.size() >= 2 && c.size()%2 == 0);
  
  Eigen::ArrayXi dofs = Eigen::ArrayXi::LinSpaced(N, 0, N-1);
  Eigen::ArrayXi map = -Eigen::ArrayXi::Ones(N);
  
  size_t offset = 0;
  map.head(c[0]) = dofs.head(c[0]);
  for (size_t i = 1; i < c.size()-1; i += 2){
    offset += c[i] - c[i-1];
    map.segment(c[i], c[i+1]-c[i]) = dofs.segment(c[i], c[i+1]-c[i]) - offset;
  }
  offset += c[c.size()-1]-c[c.size()-2];
  map.tail(N-c[c.size()-1]) = dofs.tail(N-c[c.size()-1]) - offset;
  
  return map;
 
}
 
 
 
template<typename VecType>
VecType replaceSectionsVector(
                              const VecType &V,
                              const VecType &Z,
                              const std::vector<std::pair<size_t,size_t>> &sec
                              )
{
  assert ( sec[sec.size()-1].first + sec[sec.size()-1].second <= size_t(V.rows()) );
  
  VecType val = V;
  size_t posZ = 0;
  for (size_t i=0; i < sec.size(); i++){
    assert( posZ + sec[i].second <= size_t(Z.rows()) );
    
    val.segment(sec[i].first, sec[i].second) = Z.segment(posZ, sec[i].second);
    posZ += sec[i].second;
  }
  return val;
}
 
 
#endif