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 "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 std::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)
          {
            free_with_check(values);
            values = nullptr;

            if (other.values == nullptr)
              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 std::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 std::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, TrfShape& ortho_svals, ElementMode2D mode) {};

      protected:

        ProjBasedSelector(CandList cand_list, int max_order, Shapeset* shapeset, const Range& vertex_order, const 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];

        double error_weight_h;
        double error_weight_p;
        double error_weight_aniso;


        virtual void evaluate_cands_error(std::vector<Cand>& candidates, Element* e, MeshFunction<Scalar>* rsln);


        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, MeshFunction<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, int* sons, std::vector<TrfShapeExp>** sub_nonortho_svals, std::vector<TrfShapeExp>** sub_ortho_svals, const typename OptimumSelector<Scalar>::CandsInfo& info, CandElemProjError errors_squared, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]);

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


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

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

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


        virtual void precalc_ref_solution(int inx_son, MeshFunction<Scalar>* rsln, Element* element, int intr_gip_order, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]) = 0;

        virtual void free_ref_solution_data(int inx_son, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]) = 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, int son, const ElemSubTrf& sub_trf, const ElemSubShapeFunc& sub_shape, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]) = 0;


        virtual double evaluate_error_squared_subdomain(Element* sub_elem, const ElemGIP& sub_gip, int son, const ElemSubTrf& sub_trf, const ElemProj& elem_proj, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS]) = 0;
      };
    }
  }
}
#endif