refmap.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_REFMAP_H
#define __H2D_REFMAP_H

#include "../global.h"
#include "../shapeset/precalc.h"
#include "../quadrature/quad_all.h"
#include "shapeset/shapeset_h1_all.h"

namespace Hermes
{
  namespace Hermes2D
  {
    class Element;
    class Mesh;
    namespace Views{
      class Orderizer;
      class Linearizer;
      class Vectorizer;
    };

    class HERMES_API RefMap : public Transformable
    {
    public:
      RefMap();

      ~RefMap();

      void set_quad_2d(Quad2D* quad_2d);

      Quad2D* get_quad_2d() const;

      const Quad1D* get_quad_1d() const;

      virtual void set_active_element(Element* e);

      double3* get_tangent(int edge, int order = -1);

      static void untransform(Element* e, double x, double y, double& xi1, double& xi2);

      static Element* element_on_physical_coordinates(const Mesh* mesh, double x, double y, double* x_reference = NULL, double* y_reference = NULL);

      double* get_phys_x(int order);

      double* get_phys_y(int order);

      bool is_jacobian_const() const;

      double get_const_jacobian() const;

      double* get_jacobian(int order);

      int get_inv_ref_order() const;

      double2x2* get_inv_ref_map(int order);

      double2x2* get_const_inv_ref_map();

      H1ShapesetJacobi ref_map_shapeset;
      PrecalcShapeset ref_map_pss;
    private:
      double3x2* get_second_ref_map(int order);

      void inv_ref_map_at_point(double xi1, double xi2, double& x, double& y, double2x2& m);

      void second_ref_map_at_point(double xi1, double xi2, double& x, double& y, double3x2& mm);

      virtual void push_transform(int son);

      virtual void pop_transform();

      void free();

      void force_transform(uint64_t sub_idx, Trf* ctm)
      {
        this->sub_idx = sub_idx;
        stack[top] = *ctm;
        this->ctm = stack + top;
        update_cur_node();
        if(is_const) calc_const_inv_ref_map();
      }

      Quad2D* quad_2d;

      int num_tables;

      bool is_const;

      int inv_ref_order;

      double const_jacobian;

      double2x2 const_inv_ref_map;

      static const int H2D_MAX_TABLES = g_max_quad + 1 + 4 * g_max_quad + 4;

      struct Node
      {
        double* jacobian[H2D_MAX_TABLES];
        double2x2* inv_ref_map[H2D_MAX_TABLES];
        double3x2* second_ref_map[H2D_MAX_TABLES];
        double* phys_x[H2D_MAX_TABLES];
        double* phys_y[H2D_MAX_TABLES];
        double3* tan[H2D_MAX_NUMBER_EDGES];
      };

      std::map<uint64_t, Node*> nodes;

      Node* cur_node;

      Node* overflow;

      void update_cur_node()
      {
        Node* updated_node = new Node;

        if(sub_idx > H2D_MAX_IDX) {
          delete updated_node;
          cur_node = handle_overflow();
        }
        else {
          if(nodes.insert(std::make_pair(sub_idx, updated_node)).second == false)
            delete updated_node;
          else
            init_node(updated_node);
          cur_node = nodes[sub_idx];
        }
      }

      void calc_inv_ref_map(int order);

      void calc_const_inv_ref_map();

      void calc_second_ref_map(int order);

      static bool is_parallelogram(Element* e);

      void calc_phys_x(int order);

      void calc_phys_y(int order);

      void calc_tangent(int edge, int eo);

      int calc_inv_ref_order();

      void init_node(Node* pp);

      void free_node(Node* node);

      Node* handle_overflow();

      Quad1DStd quad_1d;

      int indices[70];

      int nc;

      double2* coeffs;

      double2  lin_coeffs[H2D_MAX_NUMBER_EDGES];
      template<typename T> friend class MeshFunction;
      template<typename T> friend class DiscreteProblem;
      template<typename T> friend class DiscreteProblemLinear;
      template<typename T> friend class Solution;
      template<typename T> friend class ExactSolution;
      template<typename T> friend class ExactSolutionScalar;
      template<typename T> friend class ExactSolutionVector;
      template<typename T> friend class Adapt;
      template<typename T> friend class KellyTypeAdapt;
      friend class Views::Orderizer;
      friend class Views::Vectorizer;
      friend class Views::Linearizer;
      template<typename T> friend class Global;
      friend class VonMisesFilter;
      template<typename T> friend class Func;
      template<typename T> friend class Geom;
      friend Geom<double>* init_geom_vol(RefMap *rm, const int order);
      friend Geom<double>* init_geom_surf(RefMap *rm, SurfPos* surf_pos, const int order);
      friend Func<double>* init_fn(PrecalcShapeset *fu, RefMap *rm, const int order);
    };
  }
}
#endif