fe_geom.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-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_ELEMENT_GEOMETRY__
#define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_ELEMENT_GEOMETRY__
 
#include "lib_grid/tools/subset_handler_interface.h"
#include "lib_disc/quadrature/quadrature.h"
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#include "lib_disc/reference_element/reference_mapping.h"
#include "common/util/provider.h"
 
#include <cmath>
 
namespace ug{
 
template <  typename TElem, int TWorldDim,
      typename TTrialSpace, typename TQuadratureRule>
class FEGeometry
{
  public:
    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;
 
    typedef TTrialSpace trial_space_type;
 
    typedef TQuadratureRule quad_rule_type;
 
    static const size_t nsh = trial_space_type::nsh;
 
    static const size_t nip = quad_rule_type::nip;
 
    static const bool staticLocalData = true;
 
  public:
    FEGeometry();
 
    size_t num_ip() const {return nip;}
 
    size_t num_sh() const {return nsh;}
 
    number weight(size_t ip) const {return m_vDetJ[ip] * m_rQuadRule.weight(ip);}
 
    const MathVector<dim>& local_ip(size_t ip) const {return m_rQuadRule.point(ip);}
 
    const MathVector<worldDim>& global_ip(size_t ip) const
    {
      UG_ASSERT(ip < nip, "Wrong ip.");
      return m_vIPGlobal[ip];
    }
 
    const MathVector<dim>* local_ips() const {return m_rQuadRule.points();}
 
    const MathVector<worldDim>* global_ips() const{return &m_vIPGlobal[0];}
 
    number shape(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
      return m_vvShape[ip][sh];
    }
 
    const MathVector<dim>& local_grad(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
      return m_vvGradLocal[ip][sh];
    }
 
    const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
      return m_vvGradGlobal[ip][sh];
    }
 
  public:
    void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
    void update_local(ReferenceObjectID roid, const LFEID& lfeID){
      update_local(roid, lfeID, 2*lfeID.order() + 1);
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner)
    {
      update(pElem, vCorner, LFEID(), m_rQuadRule.order());
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
                const LFEID& lfeID, size_t orderQuad);
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
                const LFEID& lfeID){
      update(pElem, vCorner, lfeID, 2*lfeID.order() + 1);
    }
 
  protected:
    TElem* m_pElem;
 
    const quad_rule_type& m_rQuadRule;
 
    const trial_space_type& m_rTrialSpace;
 
    ReferenceMapping<ref_elem_type, worldDim> m_mapping;
 
    MathVector<worldDim> m_vIPGlobal[nip];
 
    number m_vvShape[nip][nsh];
 
    MathVector<dim> m_vvGradLocal[nip][nsh];
 
    MathVector<worldDim> m_vvGradGlobal[nip][nsh];
 
    MathMatrix<worldDim,dim> m_vJTInv[nip];
 
    number m_vDetJ[nip];
};
 
 
 
template <int TWorldDim, int TRefDim = TWorldDim>
class DimFEGeometry
{
  public:
    static const int dim = TRefDim;
 
    static const int worldDim = TWorldDim;
 
    static const bool staticLocalData = false;
 
  public:
    DimFEGeometry();
 
    DimFEGeometry(size_t order, LFEID lfeid);
 
    DimFEGeometry(ReferenceObjectID roid, size_t order, LFEID lfeid);
 
    size_t num_ip() const {return m_nip;}
 
    size_t num_sh() const {return m_nsh;}
 
    number weight(size_t ip) const
    {
      UG_ASSERT(ip < m_nip, "Wrong index");
      return m_vDetJ[ip] * m_vQuadWeight[ip];
    }
 
    const MathVector<dim>& local_ip(size_t ip) const
    {
      UG_ASSERT(ip < m_nip, "Wrong index");
      return m_vIPLocal[ip];
    }
 
    const MathVector<worldDim>& global_ip(size_t ip) const
    {
      UG_ASSERT(ip < m_vIPGlobal.size(), "Wrong ip.");
      return m_vIPGlobal[ip];
    }
 
    const MathVector<dim>* local_ips() const {return m_vIPLocal;}
 
    const MathVector<worldDim>* global_ips() const{return &m_vIPGlobal[0];}
 
    number shape(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < m_vvShape.size(), "Wrong index");
      UG_ASSERT(sh < m_vvShape[ip].size(), "Wrong index");
      return m_vvShape[ip][sh];
    }
 
    const MathVector<dim>& local_grad(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < m_vvGradLocal.size(), "Wrong index");
      UG_ASSERT(sh < m_vvGradLocal[ip].size(), "Wrong index");
      return m_vvGradLocal[ip][sh];
    }
 
    const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
    {
      UG_ASSERT(ip < m_vvGradGlobal.size(), "Wrong index");
      UG_ASSERT(sh < m_vvGradGlobal[ip].size(), "Wrong index");
      return m_vvGradGlobal[ip][sh];
    }
 
    void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
    void update_local(ReferenceObjectID roid, const LFEID& lfeID){
      update_local(roid, lfeID, 2*lfeID.order() + 1);
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner)
    {
      update(pElem, vCorner, m_lfeID, m_quadOrder);
    }
 
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
          const LFEID& lfeID, size_t orderQuad);
    void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
          const LFEID& lfeID){
      update(pElem, vCorner, lfeID, 2*lfeID.order() + 1);
    }
 
  public:
    class BF
    {
      public:
        BF() {}
 
        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 vDetJ[ip] * 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];}
 
      protected:
        friend class DimFEGeometry<worldDim, dim>;
 
        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
        std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
        std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
 
        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)
    };
 
    void update_boundary_faces(GridObject* pElem,
                   const MathVector<worldDim>* vCornerCoords,
                   size_t orderQuad,
                   const ISubsetHandler* ish);
 
  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;
    }
 
  protected:
    std::map<int, std::vector<BF> > m_mapVectorBF;
    std::vector<BF> m_vEmptyVectorBF;
 
  protected:
    ReferenceObjectID m_roid;
 
    GridObject* m_pElem;
 
    int m_quadOrder;
 
    LFEID m_lfeID;
 
    size_t m_nip;
 
    const MathVector<dim>* m_vIPLocal;
 
    const number* m_vQuadWeight;
 
    std::vector<MathVector<worldDim> > m_vIPGlobal;
 
    std::vector<MathMatrix<worldDim,dim> > m_vJTInv;
 
    std::vector<number> m_vDetJ;
 
    size_t m_nsh;
 
    std::vector<std::vector<number> > m_vvShape;
 
    std::vector<std::vector<MathVector<dim> > > m_vvGradLocal;
 
    std::vector<std::vector<MathVector<worldDim> > > m_vvGradGlobal;
};
 
} // end namespace ug
 
#include "fe_geom_impl.h"
 
#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_ELEMENT_GEOMETRY__ */