proj_based_selector.h

// This file is part of Hermes2D.
//
// Hermes2D is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// Hermes2D 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D.  If not, see <http://www.gnu.org/licenses/>.

#ifndef __H2D_REFINEMENT_PROJ_BASED_SELECTOR_H
#define __H2D_REFINEMENT_PROJ_BASED_SELECTOR_H

#include "../global.h"
#include "optimum_selector.h"

using namespace Hermes::Algebra::DenseMatrixOperations;
namespace Hermes
{
  namespace Hermes2D
  {
    namespace RefinementSelectors {

      typedef double CandElemProjError[H2DRS_MAX_ORDER + 2][H2DRS_MAX_ORDER + 2];


      template<typename Scalar>
      class HERMES_API ProjBasedSelector : public OptimumSelector<Scalar> {
      protected:
        class TrfShapeExp;

      public: //API
        virtual ~ProjBasedSelector();


        void set_error_weights(double weight_h = H2DRS_DEFAULT_ERR_WEIGHT_H, double weight_p = H2DRS_DEFAULT_ERR_WEIGHT_P, double weight_aniso = H2DRS_DEFAULT_ERR_WEIGHT_ANISO);

        double get_error_weight_h() const;
        double get_error_weight_p() const;
        double get_error_weight_aniso() const;

        typedef Hermes::vector<TrfShapeExp> TrfShape[H2D_TRF_NUM];

        bool* cached_shape_vals_valid; 
        TrfShape* cached_shape_ortho_vals; 
        TrfShape* cached_shape_vals; 

      protected: //evaluated shape basis

        class TrfShapeExp {
        public:
          TrfShapeExp();

          virtual ~TrfShapeExp();

          const TrfShapeExp& operator = (const TrfShapeExp& other)
          {
            delete [] values; values = NULL;
            if(other.values == NULL)
              throw Exceptions::Exception("Unable to assign a non-empty values. Use references instead.");
            return *this;
          }
        private:
          int num_gip; 
          int num_expansion; 
          double** values; 


          void allocate(int num_expansion, int num_gip);


          inline double* operator[](int inx_expansion);


          inline bool empty() const;

          template<typename T> friend class ProjBasedSelector;
          template<typename T> friend class L2ProjBasedSelector;
          template<typename T> friend class H1ProjBasedSelector;
          template<typename T> friend class HcurlProjBasedSelector;
          template<typename T> friend class Adapt;
        };


        virtual void precalc_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const Hermes::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, TrfShape& svals, ElementMode2D mode) {};


        virtual void precalc_ortho_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const Hermes::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, TrfShape& ortho_svals, ElementMode2D mode) {};

      protected:

        ProjBasedSelector(CandList cand_list, double conv_exp, int max_order, Shapeset* shapeset, const typename OptimumSelector<Scalar>::Range& vertex_order, const typename OptimumSelector<Scalar>::Range& edge_bubble_order);

      protected: //internal logic

        bool warn_uniform_orders;

      protected: //error evaluation
        static const int H2DRS_VALCACHE_INVALID = 0; 
        static const int H2DRS_VALCACHE_VALID = 1; 
        static const int H2DRS_VALCACHE_USER = 2; 

        template<typename T>
        struct ValueCacheItem {

          bool is_valid() const { return state != H2DRS_VALCACHE_INVALID; };


          void mark(int new_state = H2DRS_VALCACHE_VALID) { state = new_state; };


          void set(T new_value) { value = new_value; };


          T get() const { return value; };


          ValueCacheItem(const T& value = 0, const int state = H2DRS_VALCACHE_INVALID) : value(value), state(state) {}; 
        private:
          T value; 
          int state; 
        };

        typedef double** ProjMatrixCache[H2DRS_MAX_ORDER + 2][H2DRS_MAX_ORDER + 2];


        ProjMatrixCache proj_matrix_cache[H2D_NUM_MODES];


        Hermes::vector< ValueCacheItem<Scalar> > nonortho_rhs_cache;
        Hermes::vector< ValueCacheItem<Scalar> > ortho_rhs_cache;

        double error_weight_h; 
        double error_weight_p; 
        double error_weight_aniso; 


        virtual void evaluate_cands_error(Element* e, Solution<Scalar>* rsln, double* avg_error, double* dev_error);


        virtual void calc_projection_errors(Element* e, const typename OptimumSelector<Scalar>::CandsInfo& info_h, const typename OptimumSelector<Scalar>::CandsInfo& info_p, const typename OptimumSelector<Scalar>::CandsInfo& info_aniso, Solution<Scalar>* rsln, CandElemProjError herr[H2D_MAX_ELEMENT_SONS], CandElemProjError perr, CandElemProjError anisoerr[H2D_MAX_ELEMENT_SONS]);


        void calc_error_cand_element(const ElementMode2D mode, double3* gip_points, int num_gip_points, const int num_sub, Element** sub_domains, Trf** sub_trfs, Scalar*** sub_rvals, Hermes::vector<TrfShapeExp>** sub_nonortho_svals, Hermes::vector<TrfShapeExp>** sub_ortho_svals, const typename OptimumSelector<Scalar>::CandsInfo& info, CandElemProjError errors_squared);

      protected: //projection
        struct ElemProj {
          int* shape_inxs; 
          int num_shapes; 
          Hermes::vector<TrfShapeExp>& svals; 
          Scalar* shape_coeffs; 
          int max_quad_order; 
        };


        struct ElemGIP {
          double3* gip_points; 
          int num_gip_points; 
          Scalar** rvals; 
        };

        struct ElemSubTrf {
          Trf* trf; 
          double coef_mx; 
          double coef_my; 
        };

        struct ElemSubShapeFunc {
          int inx; 
          TrfShapeExp& svals; 
        };


        virtual Scalar** precalc_ref_solution(int inx_son, Solution<Scalar>* rsln, Element* element, int intr_gip_order) = 0;


        virtual double** build_projection_matrix(double3* gip_points, int num_gip_points, const int* shape_inx, const int num_shapes, ElementMode2D) = 0;


        virtual Scalar evaluate_rhs_subdomain(Element* sub_elem, const ElemGIP& sub_gip, const ElemSubTrf& sub_trf, const ElemSubShapeFunc& sub_shape) = 0;


        virtual double evaluate_error_squared_subdomain(Element* sub_elem, const ElemGIP& sub_gip, const ElemSubTrf& sub_trf, const ElemProj& elem_proj) = 0;
      };
    }
  }
}
#endif