fvcr_geom.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2015:  G-CSC, Goethe University Frankfurt
 * Author: Christian Wehner
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * 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 Lesser General Public License for more details.
 */
 
/*
 * Node centered finite volume geometry for Crouzeix-Raviart-Elements
 */
 
#ifndef __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__
#define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__
 
// extern libraries
#include <cmath>
#include <map>
#include <vector>
 
// other ug4 modules
#include "common/common.h"
 
// library intern includes
#include "lib_disc/reference_element/reference_element.h"
#include "lib_disc/reference_element/reference_mapping.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#include "lib_disc/reference_element/reference_element_traits.h"
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/local_finite_element/crouzeix-raviart/crouzeix_raviart.h"
#include "lib_disc/quadrature/gauss/gauss_quad.h"
#include "lib_disc/common/geometry_util.h"
#include "lib_disc/domain_util_impl.h"
 
#include "fv_geom_base.h"
#include "fv_util.h"
#include "hfvcr_geom.h"
 
namespace ug{
 
/* crfv traits */
 
template <int TWorldDim,int nrfaceco> struct crfv_traits
{
    typedef void scv_type;
    typedef void face_type;
  static const size_t maxNumSCVF;
  static const size_t maxNumSCV;
  static const size_t maxNSH;
  static const size_t maxNumCo;
};
 
template <> struct crfv_traits<1, 1>
{
    typedef ReferenceEdge scv_type;
    typedef ReferenceVertex face_type;
  static const size_t maxNumSCVF = 1;
  static const size_t maxNumSCV = 2;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 2;
};
 
template <> struct crfv_traits<1, 2>
{
  typedef ReferenceEdge scv_type;
    typedef ReferenceVertex face_type;
  static const size_t maxNumSCVF = 1;
  static const size_t maxNumSCV = 2;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 2;
};
 
template <> struct crfv_traits<2, 2>
{
    typedef ReferenceTriangle scv_type;
    typedef ReferenceEdge face_type;
  static const size_t maxNumSCVF = 4;
  static const size_t maxNumSCV = 4;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 4;
};
 
template <> struct crfv_traits<2, 3>
{
    typedef ReferenceTriangle scv_type;
    typedef ReferenceEdge face_type;
  static const size_t maxNumSCVF = 4;
  static const size_t maxNumSCV = 4;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 4;
};
 
template <> struct crfv_traits<3, 3>
{
  typedef ReferenceTetrahedron scv_type;
    typedef ReferenceTriangle face_type;
  static const size_t maxNumSCVF = 10;
  static const size_t maxNumSCV = 6;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 8;
};
 
template <> struct crfv_traits<3, 4>
{
    typedef ReferencePyramid scv_type;
    typedef ReferenceQuadrilateral face_type;
  static const size_t maxNumSCVF = 10;
  static const size_t maxNumSCV = 6;
  static const size_t maxNSH = maxNumSCV;
  static const size_t maxNumCo = 8;
};
 
// Dim-dependent Finite Volume Geometry
 
//
//  basics:
//  Unknowns are centered in element vertices(1d)/edges(2d)/faces(3d)
//  scv i corners are given by face i corners + element barycenter
//  scvf i corners are given by edge i + element barycenter
//
//  important differences to standard linear Lagrange FV class:
//
//  new scv class parameter "Normal" is normal on associated face  (this is computed for
//  Navier-Stokes-Crouzeix-Raviart-staggered grid continuity equation discretization)
//  midpoint computations and related parameters eliminated (not necessary for scv computation)
//  parameter nodeID refers to node located in edge/face for the unknown, not to corner node
//  ip of scv is not corner but associated edge/face midpoint
//  location of unknowns (edge/face barycenter) is stored in m_vLocUnkCoords, m_vGlobUnkCoords
//
template <int TDim, int TWorldDim = TDim>
class DimCRFVGeometry : public FVGeometryBase
{
  public:
    typedef crfv_traits<TDim, TDim> traits;
    typedef hcrfv_traits<TDim, TDim> htraits;
 
  public:
    static const int dim = TDim;
 
    static const int worldDim = TWorldDim;
 
    static const bool usesHangingNodes = false;
 
    static const bool staticLocalData = false;
 
  public:
    static const int order = 1;
 
    static const size_t maxNumSCV = htraits::maxNumSCV;
 
    typedef typename crfv_traits<TWorldDim,TWorldDim>::scv_type scv_type0;
      typedef typename crfv_traits<TWorldDim,TWorldDim>::face_type face_type0;
      typedef typename crfv_traits<TWorldDim,TWorldDim+1>::scv_type scv_type1;
      typedef typename crfv_traits<TWorldDim,TWorldDim+1>::face_type face_type1;
 
    static const size_t maxNumSCVF = htraits::maxNumSCVF;
 
    static const size_t maxNSH = htraits::maxNSH;
 
    static const size_t maxNumCo = traits::maxNumCo;
 
    static const size_t nip = 1;
 
  // local/global element barycenter
    MathVector<dim> localBary;
    MathVector<worldDim> globalBary;
 
  public:
 
    class SCVF
    {
      public:
        static const size_t numCo = dim;
 
      public:
        SCVF() {}
 
        inline size_t from() const {return From;}
 
        inline size_t to() const {return To;}
 
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return localIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return globalIP;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return numSH;}
 
        inline number shape(size_t sh) const {return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
        inline size_t num_corners() const {return numCo;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
      
      public:
        inline void operator = (const SCVF &s)
        {
          From = s.from();
          To   = s.to();
          localIP = s.local_ip();
          globalIP = s.global_ip();
          Normal = s.normal();
          JtInv = s.JTInv();
          detj = s.detJ();
          numSH = s.num_sh();
          for (size_t i=0;i<numSH;i++){
            vLocalGrad[i] = s.global_grad(i);
            vGlobalGrad[i] = s.global_grad(i);
            vShape[i] = s.shape(i);
          }
          for (size_t i=0;i<numCo;i++){
            vLocPos[i] = s.local_corner(i);
            vGloPos[i] = s.global_corner(i);
          }
        }
 
      private:
      //  let outer class access private members
        friend class DimCRFVGeometry<dim, worldDim>;
 
      // This scvf separates the scv with the ids given in "from" and "to"
      // The computed normal points in direction from->to
        size_t From, To;
 
      //  The normal on the SCVF pointing (from -> to)
        MathVector<worldDim> Normal; // normal (incl. area)
 
        MathVector<dim> vLocPos[numCo]; // local corners of scvf
        MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
 
      // scvf part
        MathVector<dim> localIP; // local integration point
        MathVector<worldDim> globalIP; // global integration point
 
      // shapes and derivatives
        size_t numSH;
        number vShape[maxNSH]; // shapes at ip
        MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
 
    class SCV
    {
      public:
        static const size_t maxNumCo = 5;
 
      public:
        SCV() : numCorners(maxNumCo) {};
 
        inline number volume() const {return Vol;}
 
        inline size_t num_corners() const {return numCorners;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
 
        inline size_t node_id() const {return nodeID;}
        
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return vLocIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return vGlobIP;}
        
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return numSH;}
 
        inline number shape(size_t sh) const {return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
      private:
      //  let outer class access private members
        friend class DimCRFVGeometry<dim, worldDim>;
        
      //  The normal on the associated face pointing outward
        MathVector<worldDim> Normal; // normal (incl. area)
 
      //  node id of associated node
        size_t nodeID;
 
      //  volume of scv
        number Vol;
 
      //  number of corners of this element
        int numCorners;
 
      //  local and global positions of this element
        MathVector<dim> vLocPos[maxNumCo]; // local position of node
        MathVector<worldDim> vGloPos[maxNumCo]; // global position of node
        MathVector<dim> vLocIP;
        MathVector<dim> vGlobIP;
 
      // shapes and derivatives
        size_t numSH;
        number vShape[maxNSH]; // shapes at ip
        MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
 
    class BF
    {
      public:
        static const size_t maxNumCo=4;
 
      public:
        BF() {}
 
        inline size_t node_id() const {return nodeID;}
 
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return localIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return globalIP;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
 
        inline number volume() const {return Vol;}
 
        inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return numSH;}
 
        inline number shape(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
        inline size_t num_corners() const {return numCo;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
 
      private:
        friend class DimCRFVGeometry<dim, worldDim>;
 
      //  id of scv this bf belongs to
        size_t nodeID;
 
        size_t numCo;
 
        MathVector<dim> vLocPos[maxNumCo]; // local corners of bf
        MathVector<worldDim> vGloPos[maxNumCo]; // global corners of bf
 
      //  scvf part
        MathVector<dim> localIP; // local integration point
        MathVector<worldDim> globalIP; // global integration point
        MathVector<worldDim> Normal; // normal (incl. area)
        number Vol; // volume of bf
 
      //  shapes and derivatives
        size_t numSH;
        number vShape[maxNSH]; // shapes at ip
        MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
    
    class CONSTRAINED_DOF
    {
      public:
        static const size_t maxNumConstrainingDofs = 4;
        inline size_t constraining_dofs_index(size_t i) const{
          return cDofInd[i];
        }
        inline number constraining_dofs_weight(size_t i) const{
          return cDofWeights[i];
        }
        inline size_t index() const{
          return i;
        }
        inline size_t num_constraining_dofs() const{
          return numConstrainingDofs;
        }
      private:
        //  let outer class access private members
        friend class DimCRFVGeometry<dim, worldDim>;
 
        //  constraining dofs indices
        size_t cDofInd[maxNumConstrainingDofs];
        //  weights
        number cDofWeights[maxNumConstrainingDofs];
        //  local index of dof in element
        size_t i;
        //  nr of constraining dofs
        size_t numConstrainingDofs;
    };
  
    class HandledEdge
    {
    public:
      size_t index;
      size_t associatedSCV[2];
      size_t scvfIndex;
      // indicates if the handled side is from or to
      bool from;
    };
 
  public:
    DimCRFVGeometry() : m_pElem(NULL), m_roid(ROID_UNKNOWN) {};
 
    void update_local_data();
 
    void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
                const ISubsetHandler* ish = NULL);
          
    void update_hanging(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
                const ISubsetHandler* ish = NULL,bool keepSCV=false,bool keepSCVF=false);
          
    void update_geometric_data(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
                const ISubsetHandler* ish = NULL);
 
    void update_boundary_faces(GridObject* elem,
                               const MathVector<worldDim>* vCornerCoords,
                               const ISubsetHandler* ish = NULL);
 
    inline size_t num_scvf() const {return m_numSCVF;};
 
    inline const SCVF& scvf(size_t i) const
      {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
 
    inline size_t num_scv() const {return m_numSCV;}
 
    inline const SCV& scv(size_t i) const
      {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
      
    inline size_t num_constrained_dofs() const {return m_numConstrainedDofs;}
 
    inline const CONSTRAINED_DOF& constrained_dof(size_t i) const
      {UG_ASSERT(i < m_numConstrainedDofs, "Invalid Index."); return m_vCD[i];}
  
    inline size_t num_constrained_scvf(){return m_numConstrainedSCVF;}
 
    inline size_t num_sh() const {return m_nsh;};
 
  public:
    size_t num_scvf_ips() const {return m_numSCVF;}
 
    const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
 
    const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
 
    size_t num_scv_ips() const {return m_numSCV;}
 
    const MathVector<dim>* scv_local_ips() const {return m_vLocUnkCoords;}
 
    const MathVector<worldDim>* scv_global_ips() const {return m_vGlobUnkCoords;}
    
    const MathVector<dim> local_bary() const {return localBary;}
    
    const MathVector<worldDim> global_bary() const {return globalBary;}
 
  protected:
  //  global and local ips on SCVF
    MathVector<worldDim> m_vGlobSCVF_IP[maxNumSCVF];
    MathVector<dim> m_vLocSCVF_IP[maxNumSCVF];
  // coord of location for unknowns in faces (edge/face barycenter)
    MathVector<worldDim> m_vGlobUnkCoords[maxNumSCV];
    MathVector<dim> m_vLocUnkCoords[maxNumSCV];
 
  public:
    inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
 
    inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
 
    inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
 
    inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
 
    inline size_t num_bf() const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      size_t num = 0;
      for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
        num += (*it).second.size();
      return num;
    }
 
    inline size_t num_bf(int subsetIndex) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      if(it == m_mapVectorBF.end()) return 0;
      else return (*it).second.size();
    }
 
    inline const BF& bf(int subsetIndex, size_t i) const
    {
      UG_ASSERT(i < num_bf(subsetIndex), "Invalid index.");
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
      return (*it).second[i];
    }
 
    inline const std::vector<BF>& bf(int subsetIndex) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
      return (*it).second;
    }
 
  protected:
    std::map<int, std::vector<BF> > m_mapVectorBF;
    
    CONSTRAINED_DOF m_vCD[maxNumSCV];
    
    std::vector<HandledEdge> handledEdges;
    
  private:
    GridObject* m_pElem;
 
    ReferenceObjectID m_roid;
    
    DimReferenceElement<dim> m_rRefElem;
 
    size_t m_numSCV;
 
    size_t m_numSCVF;
 
    size_t m_nsh;
 
    MathVector<dim> m_ipCoord[maxNumSCVF];
    MathVector<dim> cornerCoord[maxNumCo];
 
    SCVF m_vSCVF[maxNumSCVF];
  
    SCVF m_vConstrainedSCVF[maxNumSCVF];
 
    SCV m_vSCV[maxNumSCV];
    
    size_t m_numDofs;
    
    size_t m_numConstrainedDofs;
  
    size_t m_numConstrainedSCVF;
 
    static const size_t deleted = 117;
};
 
// Dimension-indipendent Finite Volume Geometry
 
template <  typename TElem, int TWorldDim>
class CRFVGeometry : public FVGeometryBase
{
  public:
    typedef TElem elem_type;
 
    typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
 
    static const int dim = ref_elem_type::dim;
 
    static const int worldDim = TWorldDim;
 
    static const bool usesHangingNodes = false;
 
    static const bool staticLocalData = true;
 
    typedef CrouzeixRaviartLSFS<ref_elem_type> local_shape_fct_set_type;
 
    static const size_t nsh = local_shape_fct_set_type::nsh;
 
    static const size_t numSCV = nsh;
 
    static const size_t numSCVF = ref_elem_type::numEdges;
    
    static const size_t maxNumSCVF = ref_elem_type::numEdges;
 
  public:
    static const int order = 1;
 
    typedef typename crfv_traits<TWorldDim,TWorldDim>::scv_type scv_type0;
      typedef typename crfv_traits<TWorldDim,TWorldDim>::face_type face_type0;
      typedef typename crfv_traits<TWorldDim,TWorldDim+1>::scv_type scv_type1;
      typedef typename crfv_traits<TWorldDim,TWorldDim+1>::face_type face_type1;
 
    static const size_t nip = 1;
 
  // local/global element barycenter
    MathVector<dim> localBary;
    MathVector<worldDim> globalBary;
 
  public:
 
    class SCVF
    {
      public:
        static const size_t numCo = dim;
 
      public:
        SCVF() {}
 
        inline size_t from() const {return From;}
 
        inline size_t to() const {return To;}
 
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return localIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return globalIP;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t sh) const {return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
        inline size_t num_corners() const {return numCo;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
 
      private:
      //  let outer class access private members
        friend class CRFVGeometry<TElem, TWorldDim>;
 
      // This scvf separates the scv with the ids given in "from" and "to"
      // The computed normal points in direction from->to
        size_t From, To;
 
      //  The normal on the SCVF pointing (from -> to)
        MathVector<worldDim> Normal; // normal (incl. area)
 
        MathVector<dim> vLocPos[numCo]; // local corners of scvf
        MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
 
      // scvf part
        MathVector<dim> localIP; // local integration point
        MathVector<worldDim> globalIP; // global integration point
 
      // shapes and derivatives
        number vShape[nsh]; // shapes at ip
        MathVector<dim> vLocalGrad[nsh]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
 
    class SCV
    {
      public:
        static const size_t maxNumCo = 5;
 
      public:
        SCV() : numCorners(maxNumCo) {};
 
        inline number volume() const {return Vol;}
 
        inline size_t num_corners() const {return numCorners;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
 
        inline size_t node_id() const {return nodeID;}
 
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return vLocIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return vGlobIP;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t sh) const {return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
      private:
      //  let outer class access private members
        friend class CRFVGeometry<TElem, TWorldDim>;
 
      //  The normal on the associated face pointing outward
        MathVector<worldDim> Normal; // normal (incl. area)
 
      //  node id of associated node
        size_t nodeID;
 
      //  volume of scv
        number Vol;
 
      //  number of corners of this element
        int numCorners;
 
      //  local and global positions of this element
        MathVector<dim> vLocPos[maxNumCo]; // local position of node
        MathVector<worldDim> vGloPos[maxNumCo]; // global position of node
 
        MathVector<dim> vLocIP; // local position of node
        MathVector<worldDim> vGlobIP; // global position of node
 
      // shapes and derivatives
        number vShape[nsh]; // shapes at ip
        MathVector<dim> vLocalGrad[nsh]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
 
    class BF
    {
      public:
        static const size_t maxNumCo=4;
 
      public:
        BF() {}
 
        inline size_t node_id() const {return nodeID;}
 
        inline size_t num_ip() const {return nip;}
 
        inline const MathVector<dim>& local_ip() const {return localIP;}
 
        inline const MathVector<worldDim>& global_ip() const {return globalIP;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
 
        inline number volume() const {return Vol;}
 
        inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
 
        inline number detJ() const {return detj;}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
 
        inline const number* shape_vector() const {return vShape;}
 
        inline const MathVector<dim>& local_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
 
        inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
 
        inline const MathVector<worldDim>& global_grad(size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
 
        inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
 
        inline size_t num_corners() const {return numCo;}
 
        inline const MathVector<dim>& local_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
 
        inline const MathVector<worldDim>& global_corner(size_t co) const
          {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
 
      private:
        friend class CRFVGeometry<TElem, TWorldDim>;
 
      //  id of scv this bf belongs to
        size_t nodeID;
 
        size_t numCo;
 
        MathVector<dim> vLocPos[maxNumCo]; // local corners of bf
        MathVector<worldDim> vGloPos[maxNumCo]; // global corners of bf
 
      //  scvf part
        MathVector<dim> localIP; // local integration point
        MathVector<worldDim> globalIP; // global integration point
        MathVector<worldDim> Normal; // normal (incl. area)
        number Vol; // volume of bf
 
      //  shapes and derivatives
        number vShape[nsh]; // shapes at ip
        MathVector<dim> vLocalGrad[nsh]; // local grad at ip
        MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
        MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
        number detj; // Jacobian det at ip
    };
 
  public:
    CRFVGeometry();
 
    void update_local_data();
 
    void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
                const ISubsetHandler* ish = NULL);
 
    void update_boundary_faces(GridObject* elem,
                               const MathVector<worldDim>* vCornerCoords,
                               const ISubsetHandler* ish = NULL);
                   
    const MathVector<worldDim>* corners() const {return m_vCo;}             
 
    inline size_t num_scvf() const {return numSCVF;};
 
    inline const SCVF& scvf(size_t i) const
      {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
 
    inline size_t num_scv() const {return numSCV;}
 
    inline const SCV& scv(size_t i) const
      {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
 
    inline size_t num_sh() const {return nsh;};
 
  public:
    size_t num_scvf_ips() const {return numSCVF;}
 
    const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
 
    const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
 
    size_t num_scv_ips() const {return numSCV;}
 
    const MathVector<dim>* scv_local_ips() const {return m_vLocUnkCoords;}
 
    const MathVector<worldDim>* scv_global_ips() const {return m_vGlobUnkCoords;}
 
    const MathVector<dim> local_bary() const {return localBary;}
 
    const MathVector<worldDim> global_bary() const {return globalBary;}
 
  protected:
  //  global and local ips on SCVF
    MathVector<worldDim> m_vGlobSCVF_IP[numSCVF];
    MathVector<dim> m_vLocSCVF_IP[numSCVF];
   // coord of location for unknowns in faces (edge/face barycenter)
    MathVector<worldDim> m_vGlobUnkCoords[numSCV];
    MathVector<dim> m_vLocUnkCoords[numSCV];
    
    static const size_t numMaxCo = 8;
   // corner coordinates
    MathVector<worldDim> m_vCo[numMaxCo];
 
  public:
    inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
 
    inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
 
    inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
 
    inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
 
    inline size_t num_bf() const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      size_t num = 0;
      for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
        num += (*it).second.size();
      return num;
    }
 
    inline size_t num_bf(int subsetIndex) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      if(it == m_mapVectorBF.end()) return 0;
      else return (*it).second.size();
    }
 
    inline const BF& bf(int subsetIndex, size_t i) const
    {
      UG_ASSERT(i < num_bf(subsetIndex), "Invalid index.");
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
      return (*it).second[i];
    }
 
    inline const std::vector<BF>& bf(int subsetIndex) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(subsetIndex);
      UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
      return (*it).second;
    }
 
  protected:
    std::map<int, std::vector<BF> > m_mapVectorBF;
 
  private:
    MathVector<dim> m_ipCoord[numSCVF];
 
    SCVF m_vSCVF[numSCVF];
 
    SCV m_vSCV[numSCV];
 
    TElem* m_pElem;
 
    ReferenceMapping<ref_elem_type, worldDim> m_mapping;
 
    const ref_elem_type& m_rRefElem;
 
    const local_shape_fct_set_type& m_rTrialSpace;
    
  public:
    class CONSTRAINED_DOF
    {
      public:
        inline size_t constraining_dofs_index(size_t i) const {
          return 0;
        }
        inline size_t constraining_dofs_weight(size_t i) const{
          return 0;
        }
        inline size_t index() const{
          return 0;
        }
        inline size_t num_constraining_dofs() const{
          return 0;
        }
    };
    inline size_t num_constrained_dofs() const {return 0;}
    inline const CONSTRAINED_DOF& constrained_dof(size_t i) const
      {return cd;}
  private:
    CONSTRAINED_DOF cd;
};
 
 
}
 
#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__ */