optimum_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_OPTIMUM_SELECTOR_H
#define __H2D_REFINEMENT_OPTIMUM_SELECTOR_H

#include <ostream>
#include "order_permutator.h"
#include "selector.h"
#include "../shapeset/shapeset.h"

namespace Hermes
{
  namespace Hermes2D
  {
    namespace RefinementSelectors {
      enum CandList {
        H2D_NONE,  
        H2D_P_ISO, 
        H2D_P_ANISO, 
        H2D_H_ISO, 
        H2D_H_ANISO, 
        H2D_HP_ISO, 
        H2D_HP_ANISO_H, 
        H2D_HP_ANISO_P, 
        H2D_HP_ANISO 
      };

      enum SelOption {
        H2D_PREFER_SYMMETRIC_MESH, 
        H2D_APPLY_CONV_EXP_DOF 
      };


      extern HERMES_API const char* get_cand_list_str(const CandList cand_list);


      extern HERMES_API bool is_hp(const CandList cand_list);


      extern HERMES_API bool is_p_aniso(const CandList cand_list);


      template<typename Scalar>
      class HERMES_API OptimumSelector : public Selector<Scalar> {
      protected: //options
        bool opt_symmetric_mesh; 
        bool opt_apply_exp_dof; 
      public:

        virtual void set_option(const SelOption option, bool enable);

        struct Cand {
          double error; 
          int dofs;  
          int split; 
          int p[H2D_MAX_ELEMENT_SONS]; 
          double score; 


          Cand(const int split, const int order_elems[H2D_MAX_ELEMENT_SONS])
            : dofs(-1), split(split), score(0) {
              p[0] = order_elems[0];
              p[1] = order_elems[1];
              p[2] = order_elems[2];
              p[3] = order_elems[3];
          };


          Cand(const int split, const int order_elem0, const int order_elem1 = 0, const int order_elem2 = 0, const int order_elem3 = 0)
            : dofs(-1), split(split), score(0) {
              p[0] = order_elem0;
              p[1] = order_elem1;
              p[2] = order_elem2;
              p[3] = order_elem3;
          };


          int get_num_elems() const {
            switch (split) {
            case H2D_REFINEMENT_H: return 4;
            case H2D_REFINEMENT_P: return 1;
            case H2D_REFINEMENT_ANISO_H:
            case H2D_REFINEMENT_ANISO_V:
              return 2;
            default:
              throw Hermes::Exceptions::Exception("Invalid refinement type %d.", split);
              return -1;
              break;
            }
          }
        };


        const Hermes::vector<Cand>& get_candidates() const { return candidates; };

        int ****num_shapes;

      protected: //candidates
        struct CandsInfo {
          bool uniform_orders; 
          int min_quad_order; 
          int max_quad_order; 

          CandsInfo() : uniform_orders(true), min_quad_order(-1), max_quad_order(-1) {};


          bool is_empty() const { return (min_quad_order < 0 || max_quad_order < 0); };
        };

        CandList cand_list; 
        double conv_exp; 
        Hermes::vector<Cand> candidates; 


        void update_cands_info(CandsInfo& info_h, CandsInfo& info_p, CandsInfo& info_aniso) const;


        void append_candidates_split(const int start_quad_order, const int last_order, const int split, bool iso_p);


        virtual void create_candidates(Element* e, int quad_order, int max_ha_quad_order, int max_p_quad_order);


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


        virtual void select_best_candidate(Element* e, const double avg_error, const double dev_error, int* selected_cand, int* selected_h_cand);


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


        virtual void evaluate_cands_dof(Element* e, Solution<Scalar>* rsln);


        virtual void evaluate_cands_score(Element* e);

      private:

        static bool compare_cand_score(const Cand& a, const Cand& b);

      protected: //orders and their range
        int current_max_order; 
        int current_min_order; 


        virtual void set_current_order_range(Element* element) = 0;

      protected: //shape functions
        enum ShapeType {
          H2DST_VERTEX = 0x01, 
          H2DST_HORIZ_EDGE = 0x02, 
          H2DST_VERT_EDGE = 0x04, 
          H2DST_TRI_EDGE = 0x08, 
          H2DST_BUBBLE = 0x10 
        };

        enum ShapeTypeInt {
          H2DSI_VERTEX,
          H2DSI_HORIZ_EDGE,
          H2DSI_VERT_EDGE,
          H2DSI_TRI_EDGE,
          H2DSI_BUBBLE,
          H2DSI_ANY
        };


        struct ShapeInx {
          int order_h; 
          int order_v; 
          int inx; 
          ShapeType type; 


          ShapeInx(int order_h, int order_v, int inx, ShapeType type) : order_h(order_h), order_v(order_v), inx(inx), type(type) {};
        };

        class Range
        {
        protected:
          int lower_bound;    
          int upper_bound;    
          bool empty_range; 
          Range();
          Range(const int& lower_bound, const int& upper_bound);
          bool empty() const;
          const int& lower() const;
          const int& upper() const;
          bool is_in_closed(const Range& range) const;
          bool is_in_closed(const int& value) const;
          bool is_in_open(const int& value) const;
          void enlarge_to_include(const int& value);

          static Range make_envelope(const Range& a, const Range& b);
          template<typename T> friend class OptimumSelector;
          template<typename T> friend class ProjBasedSelector;
          template<typename T> friend class H1ProjBasedSelector;
          template<typename T> friend class L2ProjBasedSelector;
          template<typename T> friend class HcurlProjBasedSelector;
        };

        Shapeset *shapeset; 

        Hermes::vector<ShapeInx> shape_indices[H2D_NUM_MODES]; 
        int max_shape_inx[H2D_NUM_MODES]; 
        int next_order_shape[H2D_NUM_MODES][H2DRS_MAX_ORDER+1]; 
        bool has_vertex_shape[H2D_NUM_MODES]; 
        bool has_edge_shape[H2D_NUM_MODES]; 
        bool has_bubble_shape[H2D_NUM_MODES]; 


        void add_bubble_shape_index(int order_h, int order_v, std::map<int, bool>& used_shape_index, Hermes::vector<ShapeInx>& indices, ElementMode2D mode);


        void build_shape_indices(const ElementMode2D mode, const Range& vertex_order, const Range& edge_bubble_order);


        int calc_num_shapes(int mode, int order_h, int order_v, int allowed_type_mask);


        OptimumSelector(CandList cand_list, double conv_exp, int max_order, Shapeset* shapeset, const Range& vertex_order, const Range& edge_bubble_order);

      public:
        virtual ~OptimumSelector();
      protected:

        virtual bool select_refinement(Element* element, int quad_order, Solution<Scalar>* rsln, ElementToRefine& refinement); 


        virtual void generate_shared_mesh_orders(const Element* element, const int orig_quad_order, const int refinement, int tgt_quad_orders[H2D_MAX_ELEMENT_SONS], const int* suggested_quad_orders); 
      };
    }
  }
}

#endif