err_est_data.h

Go to the documentation of this file.

/*
 * Copyright (c) 2014-2015:  G-CSC, Goethe University Frankfurt
 * Authors: Dmitry Logashenko, Markus Breit
 * 
 * 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.
 */
 
/*
 *  Data shared by element discretizations for a-posteriori error estimation
 */
 
#ifndef __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__
#define __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__
 
// extern headers
#include <vector>
#include <string>
#include <limits>
 
// intern headers
#include "lib_grid/tools/surface_view.h"
#include "lib_grid/algorithms/multi_grid_util.h"
#include "lib_disc/function_spaces/integrate.h"
 
#ifdef UG_PARALLEL
  #include "lib_grid/parallelization/util/compol_attachment_reduce.h"
  #include "lib_grid/parallelization/util/compol_copy_attachment.h"
#endif
 
#include <boost/mpl/for_each.hpp>
 
 
namespace ug{
 
 
template <typename TDomain>
class IErrEstData
{
  public:
    static const int dim = TDomain::dim;
  
    IErrEstData () : m_consider(true), m_scale(1.0) {};
    
    virtual ~IErrEstData () {};
 
    virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl) = 0;
    
    virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain) = 0;
 
    virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]) = 0;
    
    virtual void release_err_est_data () = 0;
 
    bool consider_me() const {return m_consider;};
 
 
    void set_consider_me(bool b) {m_consider = b;};
 
    void set_scaling_factor(number scale) {m_scale = scale;};
 
 
    number scaling_factor() {return m_scale;}
 
  private:
  bool m_consider;
  number m_scale;
};
 
 
template <typename TDomain>
class SideFluxErrEstData : public IErrEstData<TDomain>
{
public:
    typedef TDomain domain_type;
    
    static const int dim = TDomain::dim;
    
    typedef typename domain_traits<dim>::side_type side_type;
  
public:
    SideFluxErrEstData() : IErrEstData<TDomain>() {};
 
    virtual ~SideFluxErrEstData() {};
 
//  Functions to access data
 
    number& operator()
    (
      side_type* pSide 
    )
    {
      return m_aaFluxJump[pSide];
    };
    
    ConstSmartPtr<SurfaceView>& surface_view() {return m_spSV;};
 
//  Interface virtual functions inherited from IErrEstData
 
    virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
    
    virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain);
 
    virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]) {return 0;};
    
    virtual void release_err_est_data ();
  
private:
    ANumber m_aFluxJumpOverSide;
    
    MultiGrid::AttachmentAccessor<side_type, ANumber> m_aaFluxJump;
    
    ConstSmartPtr<SurfaceView> m_spSV;
    
    GridLevel m_errEstGL;
};
 
 
 
 
template <typename TDomain>
class SideAndElemErrEstData : public IErrEstData<TDomain>
{
public:
    typedef TDomain domain_type;
 
    static const int dim = TDomain::dim;
 
    typedef typename domain_traits<dim>::side_type side_type;
    typedef typename domain_traits<dim>::element_type elem_type;
 
    typedef Attachment<std::vector<number> > attachment_type;
 
    typedef SideAndElemErrEstData<TDomain> this_type;
 
    static const int MAX_NUM_SIDES = 8;
 
public:
    SideAndElemErrEstData(size_t _sideOrder, size_t _elemOrder, const char* subsets);
    SideAndElemErrEstData(size_t _sideOrder, size_t _elemOrder,
                std::vector<std::string> subsets = std::vector<std::string>(0));
 
    virtual ~SideAndElemErrEstData() {};
 
  //  Functions to access data
 
    size_t side_order() const {return sideOrder;}
 
    size_t elem_order() const {return elemOrder;}
 
    number& operator()
    (
      side_type* pSide,   
      size_t ip   
    );
 
    number& operator()
    (
      elem_type* pElem,   
      size_t ip   
    );
 
    template <int refDim>
    const MathVector<refDim>* side_local_ips(const ReferenceObjectID roid);
 
    template <int refDim>
    const MathVector<refDim>* elem_local_ips(const ReferenceObjectID roid);
 
    MathVector<TDomain::dim>* all_side_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    MathVector<TDomain::dim>* side_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    MathVector<TDomain::dim>* elem_global_ips(GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    size_t num_side_ips(const side_type* pSide);
 
    size_t num_side_ips(const ReferenceObjectID roid);
 
    size_t first_side_ips(const ReferenceObjectID roid, const size_t side);
 
    size_t num_all_side_ips(const ReferenceObjectID roid);
 
    size_t num_elem_ips(const ReferenceObjectID roid);
 
    size_t side_ip_index(const ReferenceObjectID roid, const size_t side, const size_t ip);
 
    ConstSmartPtr<SurfaceView>& surface_view () {return m_spSV;};
 
  //  virtual functions inherited from IErrEstData
    virtual void alloc_err_est_data (ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
 
    virtual void summarize_err_est_data (SmartPtr<TDomain> spDomain);
 
    virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    virtual void release_err_est_data ();
 
    void set_type(int type) {m_type = (type<=H1_ERROR_TYPE) ? type : H1_ERROR_TYPE;}
 
protected:
    void init_quadrature();
 
    template<int refDim>
    struct GetQuadRules
    {
        GetQuadRules(QuadratureRule<refDim>** ppQuadRule, size_t quadOrder) :
          m_ppQuadRule(ppQuadRule), m_quadOrder(quadOrder) {}
        QuadratureRule<refDim>** m_ppQuadRule;
        size_t m_quadOrder;
        template< typename TElem > void operator()(TElem&)
        {
          const ReferenceObjectID roid = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
          m_ppQuadRule[roid] =
            const_cast<QuadratureRule<refDim>*>(&QuadratureRuleProvider<refDim>::get(roid, m_quadOrder));
        }
    };
 
private:
    size_t sideOrder;
    size_t elemOrder;
 
    std::vector<std::string> m_vSs;
    SubsetGroup m_ssg;
 
    QuadratureRule<dim-1>* quadRuleSide[NUM_REFERENCE_OBJECTS];
    QuadratureRule<dim>* quadRuleElem[NUM_REFERENCE_OBJECTS];
 
    std::vector<MathVector<TDomain::dim> > m_SideIPcoords[NUM_REFERENCE_OBJECTS];
 
    std::vector<MathVector<TDomain::dim> > m_sideGlobalIPcoords;
    std::vector<MathVector<TDomain::dim> > m_singleSideGlobalIPcoords;  // not the most elegant solution...
    std::vector<MathVector<TDomain::dim> > m_elemGlobalIPcoords;
 
    size_t m_sideIPsStartIndex[NUM_REFERENCE_OBJECTS][MAX_NUM_SIDES];
 
    attachment_type m_aSide;
 
    attachment_type m_aElem;
 
    MultiGrid::AttachmentAccessor<side_type, attachment_type > m_aaSide;
 
    MultiGrid::AttachmentAccessor<elem_type, attachment_type > m_aaElem;
 
    ConstSmartPtr<SurfaceView> m_spSV;
 
    GridLevel m_errEstGL;
 
    enum type {L2_ERROR_TYPE=0, H1_ERROR_TYPE};
    int m_type;
};
 
 
 
template <typename TDomain, typename TErrEstData>
class MultipleErrEstData : public IErrEstData<TDomain>
{
  public:
    static const int dim = TDomain::dim;
 
    MultipleErrEstData(ConstSmartPtr<ApproximationSpace<TDomain> > approx)
    : IErrEstData<TDomain>(), m_spApprox(approx),
      m_fctGrp(m_spApprox->dof_distribution_info())
    {
      this->set_consider_me(false);
    }
 
    virtual ~MultipleErrEstData() {};
 
    virtual void add(SmartPtr<TErrEstData> spEed, const char* fct)
    {
      // check that fct is not already contained in fctGrp
      size_t uid = m_spApprox->fct_id_by_name(fct);
      if (m_fctGrp.contains(uid))
      {
        UG_THROW("Error estimator for function '" << fct << "' can not be added\n"
             "as another error estimator object for the same function is already\n"
             "contained here.")
      }
 
      // add to function group
      try
      {
        m_fctGrp.add(fct);
      }
      UG_CATCH_THROW("Error estimator data object for function '"
              << fct << "' could not be added.");
 
      m_vEed.push_back(spEed.get());
    }
 
    size_t num() const {return m_vEed.size();};
 
    TErrEstData* get(size_t uid)
    {
      if (!m_fctGrp.contains(uid))
      {
        std::string name = m_spApprox->name(uid);
        UG_THROW("Trying to access error estimator data object "
             "for unique function index " << uid << "(aka '" << name << "')\n"
             "which is not present in this collection.")
      }
 
      return m_vEed[m_fctGrp.local_index(uid)];
    }
 
  //  inherited from IErrEstData
    virtual void alloc_err_est_data(ConstSmartPtr<SurfaceView> spSV, const GridLevel& gl);
 
    virtual void summarize_err_est_data(SmartPtr<TDomain> spDomain);
 
    virtual number get_elem_error_indicator(GridObject* elem, const MathVector<dim> vCornerCoords[]);
 
    virtual void release_err_est_data();
 
  protected:
    std::vector<TErrEstData*> m_vEed;
 
    ConstSmartPtr<ApproximationSpace<TDomain> > m_spApprox;
 
    FunctionGroup m_fctGrp;
};
 
 
template <typename TDomain>
class MultipleSideAndElemErrEstData
  : public MultipleErrEstData<TDomain, SideAndElemErrEstData<TDomain> >
{
  public:
    static const int dim = TDomain::dim;
 
    typedef typename SideAndElemErrEstData<TDomain>::side_type side_type;
    typedef typename SideAndElemErrEstData<TDomain>::elem_type elem_type;
 
    MultipleSideAndElemErrEstData(ConstSmartPtr<ApproximationSpace<TDomain> > approx)
      : MultipleErrEstData<TDomain, SideAndElemErrEstData<TDomain> >(approx),
        m_bEqSideOrder(false), m_bEqElemOrder(false) {};
 
    virtual ~MultipleSideAndElemErrEstData() {};
 
    virtual void add(SmartPtr<SideAndElemErrEstData<TDomain> > spEed, const char* fct);
 
    bool equal_side_order() const {return m_bEqSideOrder;}
 
    bool equal_elem_order() const {return m_bEqElemOrder;}
 
  protected:
      void check_equal_order();
 
      void check_equal_side_order();
 
      void check_equal_elem_order();
 
  private:
    bool m_bEqSideOrder;
    bool m_bEqElemOrder;
 
};
 
} // end of namespace ug
 
#include "err_est_data_impl.h"
 
#endif /* __H__UG__LIB_DISC__SPATIAL_DISC__ELEM_DISC__ERR_EST_DATA__ */
 
/* End of File */