fvho_geom.h

Go to the documentation of this file.

/*
 * Copyright (c) 2013-2015:  G-CSC, Goethe University Frankfurt
 * Author: Andreas Vogel
 * 
 * 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.
 */
 
#ifndef __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__
#define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__
 
// extern libraries
#include <cmath>
#include <map>
#include <vector>
 
// other ug4 modules
#include "common/common.h"
 
// library intern includes
#include "lib_grid/tools/subset_handler_interface.h"
#include "lib_disc/reference_element/reference_element.h"
#include "lib_disc/reference_element/reference_element_traits.h"
#include "lib_disc/reference_element/reference_mapping.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/quadrature/gauss/gauss_quad.h"
#include "fv_util.h"
#include "fv_geom_base.h"
 
namespace ug{
 
// FV Geometry for Reference Element Type (all orders)
 
 
template <  int TOrder, typename TElem, int TWorldDim,
      int TQuadOrder = TOrder + 1>
class FVGeometry : public FVGeometryBase
{
  private:
    static const int p = TOrder;
 
  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;
 
  public:
    static const int order = TOrder;
 
    static const size_t numSubElem = Pow<p, dim>::value;
 
 
    typedef fv1_traits<ref_elem_type, worldDim> traits;
 
    typedef typename traits::scvf_type scvf_type;
 
    static const size_t numSCVFPerSubElem = ref_elem_type::numEdges;
 
    static const size_t numSCVF = numSubElem * numSCVFPerSubElem;
 
    static const int quadOrderSCVF = TQuadOrder;
 
    typedef GaussQuadrature<scvf_type, quadOrderSCVF> scvf_quad_rule_type;
 
    static const size_t numSCVFIP = scvf_quad_rule_type::nip * numSCVF;
 
 
    typedef typename traits::scv_type scv_type;
 
    static const size_t numSCVPerSubElem = ref_elem_type::numCorners;
 
    static const size_t numSCV = numSubElem * numSCVPerSubElem;
 
    static const int quadOrderSCV = TQuadOrder;
 
    typedef GaussQuadrature<scv_type, quadOrderSCV> scv_quad_rule_type;
 
    static const size_t numSCVIP = scv_quad_rule_type::nip * numSCV;
 
 
    typedef LagrangeLSFS<ref_elem_type, p> local_shape_fct_set_type;
 
    static const size_t nsh = local_shape_fct_set_type::nsh;
 
    static const bool usesHangingNodes = false;
 
    static const bool staticLocalData = true;
 
  public:
    class SCVF
    {
      public:
        static const size_t nip = scvf_quad_rule_type::nip;
 
        static const size_t numCo = traits::NumCornersOfSCVF;
 
      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 number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];;}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip];}
 
        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 FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
 
      // 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)
 
      // ordering is:
      // 1D: edgeMidPoint
      // 2D: edgeMidPoint, CenterOfElement
      // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
        MathVector<dim> vLocPos[numCo]; // local corners of scvf
        MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
 
      // scvf part
        MathVector<dim> vLocalIP[nip]; // local integration point
        MathVector<worldDim> vGlobalIP[nip]; // global integration point
        const number* vWeight; // weight at ip
 
      // shapes and derivatives
        number vvShape[nip][nsh]; // shapes at ip
        MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
        MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
        MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
        number vDetJ[nip]; // Jacobian det at ip
        number vDetJMap[nip]; // Jacobian det at ip
    };
 
    class SCV
    {
      public:
        static const size_t nip = scv_quad_rule_type::nip;
 
        static const size_t numCo = traits::NumCornersOfSCV;
 
      public:
        SCV() {};
 
        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];}
 
        inline size_t node_id() const {return nodeId;}
 
        inline size_t num_ip() const {return nip;}
 
        inline number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip];}
 
      private:
      //  let outer class access private members
        friend class FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
 
      //  node id of associated node
        size_t nodeId;
 
      //  scv part
        MathVector<dim> vLocalIP[nip]; // local integration point
        MathVector<worldDim> vGlobalIP[nip]; // global intergration point
        const number* vWeight; // weight at ip
 
      //  local and global positions of this element
        MathVector<dim> vLocPos[numCo]; // local position of node
        MathVector<worldDim> vGloPos[numCo]; // global position of node
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scv
 
      // shapes and derivatives
        number vvShape[nip][nsh]; // shapes at ip
        MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
        MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
        MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
        number vDetJ[nip]; // Jacobian det at ip
        number vDetJMap[nip]; // Jacobian det at ip for scv integral map
    };
 
    class BF
    {
      public:
        static const size_t nip = scvf_quad_rule_type::nip;
 
        static const size_t numCo = traits::NumCornersOfSCVF;
 
      public:
        BF() {}
 
        inline size_t node_id() const {return nodeId;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
 
        inline number volume() const {return Vol;}
 
        inline size_t num_ip() const {return nip;}
 
        inline number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip];}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip];}
 
        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 FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
 
      //  id of scv this bf belongs to
        size_t nodeId;
 
      // ordering is:
      // 1D: edgeMidPoint
      // 2D: edgeMidPoint, CenterOfElement
      // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
        MathVector<dim> vLocPos[numCo]; // local corners of bf
        MathVector<worldDim> vGloPos[numCo]; // global corners of bf
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
 
      //  scvf part
        MathVector<dim> vLocalIP[nip]; // local integration point
        MathVector<worldDim> vGlobalIP[nip]; // global integration point
        const number* vWeight; // weight at ip
        MathVector<worldDim> Normal; // normal (incl. area)
        number Vol; // volume of bf
 
      //  shapes and derivatives
        number vvShape[nip][nsh]; // shapes at ip
        MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
        MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
        MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
        number vDetJ[nip]; // Jacobian det at ip
    };
 
 
  public:
    FVGeometry();
 
    void update_local_data();
 
    void update_local(ReferenceObjectID roid,
                      const LFEID& lfeID = LFEID(LFEID::LAGRANGE, worldDim, 1),
                      size_t orderQuad = TQuadOrder);
 
    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);
 
    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 numSCVFIP;}
 
    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 numSCVIP;}
 
    const MathVector<dim>* scv_local_ips() const {return m_vLocSCV_IP;}
 
    const MathVector<worldDim>* scv_global_ips() const {return m_vGlobSCV_IP;}
 
 
  protected:
  //  global and local ips on SCVF
    MathVector<worldDim> m_vGlobSCVF_IP[numSCVFIP];
    MathVector<dim> m_vLocSCVF_IP[numSCVFIP];
 
  //  global and local ips on SCVF
    MathVector<worldDim> m_vGlobSCV_IP[numSCVIP];
    MathVector<dim> m_vLocSCV_IP[numSCVIP];
 
  protected:
  //  maximum number of geom objects in a dimension
    static const int maxMid = numSCVF + 1;
 
  //  subelement
    struct SubElement
    {
    //  shape number of corners of subelement
      size_t vDoFID[ref_elem_type::numCorners];
 
    //  local and global geom object midpoints for each dimension
    //  (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
      MathVector<dim> vvLocMid[dim+1][maxMid];
      MathVector<worldDim> vvGloMid[dim+1][maxMid];
 
    //  flag is subelement has boundary sides
      bool isBndElem;
 
    //  -1 is no bnd side, >= 0 corresponding side of whole element
      std::vector<int> vElemBndSide;
    };
 
    SubElement m_vSubElem[numSubElem];
 
  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 si) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) return 0;
      else return (*it).second.size();
    }
 
    inline const BF& bf(int si, size_t i) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) UG_THROW("FVGeom: No BndSubset"<<si);
      return (*it).second[i];
    }
 
    inline const std::vector<BF>& bf(int si) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
      return (*it).second;
    }
 
    void reset_curr_elem() {m_pElem = NULL;}
 
  protected:
    std::map<int, std::vector<BF> > m_mapVectorBF;
    std::vector<BF> m_vEmptyVectorBF;
 
  private:
    TElem* m_pElem;
 
    MathVector<dim> m_vLocCorner[ref_elem_type::numCorners];
 
    SCVF m_vSCVF[numSCVF];
 
    SCV m_vSCV[numSCV];
 
    ReferenceMapping<ref_elem_type, worldDim> m_rMapping;
 
    const ref_elem_type& m_rRefElem;
 
    const local_shape_fct_set_type& m_rTrialSpace;
 
    const scvf_quad_rule_type& m_rSCVFQuadRule;
 
    const scv_quad_rule_type& m_rSCVQuadRule;
};
 
// Dimension FV Geometry (all orders)   DIM FV
 
 
template <int TWorldDim, int TDim = TWorldDim>
class DimFVGeometry : public FVGeometryBase
{
  public:
    typedef fv1_dim_traits<TDim, TWorldDim> traits;
 
    static const int dim = TDim;
 
    static const int worldDim = TWorldDim;
 
  public:
    typedef typename traits::scvf_type scvf_type;
 
    typedef typename traits::scv_type scv_type;
 
    static const bool usesHangingNodes = false;
 
    static const bool staticLocalData = false;
 
  public:
    class SCVF
    {
      public:
        static const size_t numCo = traits::NumCornersOfSCVF;
 
      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 number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathVector<worldDim>& normal() const {return Normal;}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvShape[ip][0];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvLocalGrad[ip][0];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvGlobalGrad[ip][0];}
 
        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 DimFVGeometry<worldDim, dim>;
 
      // 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)
 
      // ordering is:
      // 1D: edgeMidPoint
      // 2D: edgeMidPoint, CenterOfElement
      // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
        MathVector<dim> vLocPos[numCo]; // local corners of scvf
        MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
 
      // scvf part
        size_t nip;
        std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
        std::vector<MathVector<worldDim> > vGlobalIP; // global integration point (size: nip)
        const number* vWeight; // weight at ip
 
      // shapes and derivatives
        size_t nsh;
        std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
        std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
        std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
        std::vector<number> vDetJMap; // Jacobian det at ip (size: nip)
    };
 
    class SCV
    {
      public:
        static const size_t numCo = traits::NumCornersOfSCV;
 
      public:
        SCV() {};
 
        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];}
 
        inline size_t node_id() const {return nodeId;}
 
        inline size_t num_ip() const {return nip;}
 
        inline number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvShape[ip][0];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvLocalGrad[ip][0];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvGlobalGrad[ip][0];}
 
      private:
      //  let outer class access private members
        friend class DimFVGeometry<worldDim, dim>;
 
      //  node id of associated node
        size_t nodeId;
 
      //  local and global positions of this element
        MathVector<dim> vLocPos[numCo]; // local position of node
        MathVector<worldDim> vGloPos[numCo]; // global position of node
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scv
 
      //  scv part
        size_t nip;
        std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
        std::vector<MathVector<worldDim> > vGlobalIP; // global intergration point (size: nip)
        const number* vWeight; // weight at ip
 
      // shapes and derivatives
        size_t nsh;
        std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
        std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
        std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
        std::vector<number> vDetJMap; // Jacobian det at ip (size: nip)
    };
 
    class BF
    {
      public:
        static const size_t numCo = traits::NumCornersOfSCVF;
 
      public:
        BF() {}
 
        inline size_t node_id() const {return nodeId;}
 
        inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
 
        inline number volume() const {return Vol;}
 
        inline size_t num_ip() const {return nip;}
 
        inline number weight(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip];}
 
        inline const MathVector<dim>& local_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
 
        inline const MathVector<worldDim>& global_ip(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
 
        inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
 
        inline number detJ(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
 
        inline size_t num_sh() const {return nsh;}
 
        inline number shape(size_t ip, size_t sh) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
 
        inline const number* shape_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvShape[ip][0];}
 
        inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
 
        inline const MathVector<dim>* local_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvLocalGrad[ip][0];}
 
        inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
          {UG_ASSERT(sh < num_sh(), "Invalid index");
           UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
 
        inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
          {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvGlobalGrad[ip][0];}
 
        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 DimFVGeometry<worldDim, dim>;
 
      //  id of scv this bf belongs to
        size_t nodeId;
 
      // ordering is:
      // 1D: edgeMidPoint
      // 2D: edgeMidPoint, CenterOfElement
      // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
        MathVector<dim> vLocPos[numCo]; // local corners of bf
        MathVector<worldDim> vGloPos[numCo]; // global corners of bf
        MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
 
      //  scvf part
        size_t nip;
        std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
        std::vector<MathVector<worldDim> > vGlobalIP; // global integration point (size: nip)
        const number* vWeight; // weight at ip
        MathVector<worldDim> Normal; // normal (incl. area)
        number Vol; // volume of bf
 
      //  shapes and derivatives
        size_t nsh;
        std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
        std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
        std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
        std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
    };
 
 
  public:
    DimFVGeometry();
 
    void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
    void update_local(ReferenceObjectID roid, const LFEID& lfeID)
    {
      update_local(roid, lfeID, lfeID.order() + 1);
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords,
                const ISubsetHandler* ish = NULL)
    {
      update(pElem, vCornerCoords, m_lfeID, m_quadOrderSCV, ish);
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords,
                const LFEID& lfeID, size_t orderQuad,
                const ISubsetHandler* ish = NULL);
 
    void update_boundary_faces(GridObject* pElem,
                               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_sh() const {return m_nsh;};
 
  public:
    size_t num_scvf_ips() const {return m_numSCVFIP;}
 
    const MathVector<dim>* scvf_local_ips() const {return &(m_vLocSCVF_IP[0]);}
 
    const MathVector<worldDim>* scvf_global_ips() const {return &(m_vGlobSCVF_IP[0]);}
 
    size_t num_scv_ips() const {return m_numSCVIP;}
 
    const MathVector<dim>* scv_local_ips() const {return &(m_vLocSCV_IP[0]);}
 
    const MathVector<worldDim>* scv_global_ips() const {return &(m_vGlobSCV_IP[0]);}
 
 
  protected:
  //  global and local ips on SCVF (size: numSCVFIP)
    std::vector<MathVector<worldDim> > m_vGlobSCVF_IP;
    std::vector<MathVector<dim> > m_vLocSCVF_IP;
 
  //  global and local ips on SCVF
    std::vector<MathVector<worldDim> > m_vGlobSCV_IP;
    std::vector<MathVector<dim> > m_vLocSCV_IP;
 
  protected:
  //  miximum number of geom objects in a dimension
    static const int maxMid = traits::maxNumSCVF +1;
 
  //  subelement
    struct SubElement
    {
    //  shape number of corners of subelement
      size_t vDoFID[traits::maxNumSCV];
 
    //  local and global geom object midpoints for each dimension
    //  (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
      MathVector<dim> vvLocMid[dim+1][maxMid];
      MathVector<worldDim> vvGloMid[dim+1][maxMid];
 
    //  flag is subelement has boundary sides
      bool isBndElem;
 
    //  -1 is no bnd side, >= 0 corresponding side of whole element
      std::vector<int> vElemBndSide;
    };
 
    std::vector<SubElement> m_vSubElem;
 
  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 si) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) return 0;
      else return (*it).second.size();
    }
 
    inline const BF& bf(int si, size_t i) const
    {
      UG_ASSERT(i < num_bf(si), "Invalid index.");
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) UG_THROW("DimFVGeom: No BndSubset: "<<si);
      return (*it).second[i];
    }
 
    inline const std::vector<BF>& bf(int si) const
    {
      typename std::map<int, std::vector<BF> >::const_iterator it;
      it = m_mapVectorBF.find(si);
      if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
      return (*it).second;
    }
 
    void reset_curr_elem() {m_pElem = NULL;}
 
  protected:
    std::map<int, std::vector<BF> > m_mapVectorBF;
    std::vector<BF> m_vEmptyVectorBF;
 
  private:
    GridObject* m_pElem;
 
    ReferenceObjectID m_roid;
 
    int m_orderShape;
 
    LFEID m_lfeID;
 
    size_t m_numSubElem;
 
    size_t m_nsh;
 
    size_t m_numSCVF;
 
    size_t m_numSCVFPerSubElem;
 
    std::vector<SCVF> m_vSCVF;
 
    int m_quadOrderSCVF;
 
    size_t m_numSCVFIP;
 
 
    size_t m_numSCV;
 
    size_t m_numSCVPerSubElem;
 
    std::vector<SCV> m_vSCV;
 
    int m_quadOrderSCV;
 
    size_t m_numSCVIP;
};
 
}
 
#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__ */