function.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_FUNCTION_H
#define __H2D_FUNCTION_H

#include "transformable.h"
#include "../quadrature/quad.h"
#include "exceptions.h"

namespace Hermes
{
  namespace Hermes2D
  {
    struct SurfPos;
    class PrecalcShapeset;
    namespace RefinementSelectors{
      template<typename Scalar> class Selector;
      template<typename Scalar> class HOnlySelector;
      template<typename Scalar> class POnlySelector;
      template<typename Scalar> class OptimumSelector;
      template<typename Scalar> class ProjBasedSelector;
      template<typename Scalar> class L2ProjBasedSelector;
      template<typename Scalar> class H1ProjBasedSelector;
      template<typename Scalar> class HcurlProjBasedSelector;
    };
    template<typename Scalar> class Geom;

    namespace Views{
      class Orderizer;
      class Linearizer;
      class Vectorizer;
    };

    enum
    {
      H2D_FN_VAL_0 = 0x0001, H2D_FN_VAL_1 = 0x0040, // Function values
      H2D_FN_DX_0  = 0x0002, H2D_FN_DX_1  = 0x0080, // First derivative
      H2D_FN_DY_0  = 0x0004, H2D_FN_DY_1  = 0x0100, // First derivative
      H2D_FN_DXX_0 = 0x0008, H2D_FN_DXX_1 = 0x0200, // Second derivative
      H2D_FN_DYY_0 = 0x0010, H2D_FN_DYY_1 = 0x0400, // Second derivative
      H2D_FN_DXY_0 = 0x0020, H2D_FN_DXY_1 = 0x0800  // Second mixed derivative
    };

    const int H2D_FN_VAL = H2D_FN_VAL_0 | H2D_FN_VAL_1;
    const int H2D_FN_DX  = H2D_FN_DX_0  | H2D_FN_DX_1;
    const int H2D_FN_DY  = H2D_FN_DY_0  | H2D_FN_DY_1;
    const int H2D_FN_DXX = H2D_FN_DXX_0 | H2D_FN_DXX_1;
    const int H2D_FN_DYY = H2D_FN_DYY_0 | H2D_FN_DYY_1;
    const int H2D_FN_DXY = H2D_FN_DXY_0 | H2D_FN_DXY_1;

    const int H2D_FN_DEFAULT = H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY;            
    const int H2D_FN_ALL = H2D_FN_DEFAULT | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY; 

    const int H2D_FN_COMPONENT_0 = H2D_FN_VAL_0 | H2D_FN_DX_0 | H2D_FN_DY_0 | H2D_FN_DXX_0 | H2D_FN_DYY_0 | H2D_FN_DXY_0;
    const int H2D_FN_COMPONENT_1 = H2D_FN_VAL_1 | H2D_FN_DX_1 | H2D_FN_DY_1 | H2D_FN_DXX_1 | H2D_FN_DYY_1 | H2D_FN_DXY_1;

    template<typename Scalar>
    class HERMES_API Function : public Transformable
    {
    public:

      Function();

      virtual ~Function();

      int get_num_components() const;

      Scalar* get_fn_values(int component = 0);

      Scalar* get_dx_values(int component = 0);

      Scalar* get_dy_values(int component = 0);

      void get_dx_dy_values(Scalar*& dx, Scalar*& dy, int component = 0);

      Scalar* get_dxx_values(int component = 0);

      Scalar* get_dyy_values(int component = 0);

      Scalar* get_dxy_values(int component = 0);

      Quad2D* get_quad_2d() const;

      void set_quad_order(unsigned int order, int mask = H2D_FN_DEFAULT);

      Scalar* get_values(int a, int b);

      int get_fn_order() const;

    protected:
      virtual void set_quad_2d(Quad2D* quad_2d);

      virtual void free() = 0;

      struct Node
      {
        int mask;           

        int size;           

        Scalar* values[H2D_MAX_SOLUTION_COMPONENTS][6]; 

        Scalar data[1];       

      private:
        Node(const Node& org) {}; 

        Node& operator=(const Node& other) { return *this; }; 
      };

      virtual int get_edge_fn_order(int edge) const;

      virtual void precalculate(int order, int mask) = 0;

      int order;          

      int num_components; 

      std::map<uint64_t, LightArray<Node*>*>* sub_tables;

      LightArray<Node*>* nodes;

      Node* cur_node;

      LightArray<Node*>* overflow_nodes;

      void update_nodes_ptr();

      void force_transform(uint64_t sub_idx, Trf* ctm);

      Quad2D* quads[8]; 

      int cur_quad;     

      int total_mem;    

      int max_mem;      

      Node* new_node(int mask, int num_points); 

      virtual void  handle_overflow_idx() = 0;

      void replace_cur_node(Node* node);

      static void check_params(int component, Node* cur_node, int num_components);

      static void check_table(int component, Node* cur_node, int n, const char* msg);

      static int idx2mask[6][2];  
      template<typename T> friend class KellyTypeAdapt;
      template<typename T> friend class RefinementSelectors::H1ProjBasedSelector;
      template<typename T> friend class RefinementSelectors::L2ProjBasedSelector;
      template<typename T> friend class RefinementSelectors::HcurlProjBasedSelector;
      template<typename T> friend class Adapt;
      friend class Views::Orderizer;
      friend class Views::Vectorizer;
      friend class Views::Linearizer;

      template<typename T> friend class DiscontinuousFunc;
      template<typename T> friend class DiscreteProblem;
      template<typename T> friend class DiscreteProblemLinear;
      template<typename T> friend class Global;
      friend class CurvMap;

      template<typename T> friend class Func;
      template<typename T> friend class Geom;
      template<typename T> friend class Filter;
      template<typename T> friend class SimpleFilter;
      template<typename T> friend class DXDYFilter;
      friend class ComplexFilter;
      friend class VonMisesFilter;
      friend HERMES_API Geom<double>* init_geom_vol(RefMap *rm, const int order);
      friend HERMES_API Geom<double>* init_geom_surf(RefMap *rm, SurfPos* surf_pos, const int order);
      friend HERMES_API Func<double>* init_fn(PrecalcShapeset *fu, RefMap *rm, const int order);
      template<typename T> friend HERMES_API Func<T>* init_fn(MeshFunction<T>*fu, const int order);
    };
  }
}
#endif