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 "../mesh/mesh.h"
#include "../quadrature/quad_all.h"
#include "shapeset/shapeset_h1_all.h"

namespace Hermes
{
  namespace Hermes2D
  {
    class Element;
    class Mesh;

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

      ~RefMap();

      void set_quad_2d(Quad2D* quad_2d);

      Quad2D* get_quad_2d() 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(bool use_MeshHashGrid, MeshSharedPtr mesh, double x, double y, double* x_reference = nullptr, double* y_reference = nullptr);

      static bool is_element_on_physical_coordinates(Element* e, double x, double y, double* x_reference = nullptr, double* y_reference = nullptr);

      double* get_phys_x(int order);

      double* get_phys_y(int order);

      inline bool is_jacobian_const() const
      {
        return is_const;
      }

      inline int get_inv_ref_order() const
      {
        return inv_ref_order;
      }

      inline double get_const_jacobian() const
      {
        return const_jacobian;
      }

      inline double2x2* get_const_inv_ref_map()
      {
        return &const_inv_ref_map;
      }

      inline double* get_jacobian(int order)
      {
        if (this->is_const)
          throw Hermes::Exceptions::Exception("RefMap::get_jacobian() called with a const jacobian.");
        if (order != this->jacobian_calculated)
          this->calc_inv_ref_map(order);
        return this->jacobian;
      }

      inline double2x2* get_inv_ref_map(int order)
      {
        if (this->is_const)
          throw Hermes::Exceptions::Exception("RefMap::get_inv_ref_map() called with a const jacobian.");
        if (order != this->inv_ref_map_calculated)
          this->calc_inv_ref_map(order);
        return this->inv_ref_map;
      }

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

#ifdef H2D_USE_SECOND_DERIVATIVES
      double3x2* get_second_ref_map(int order);

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

      virtual void push_transform(int son);

      virtual void pop_transform();

      void force_transform(uint64_t sub_idx, Trf* ctm);

      static void set_element_iro_cache(Element* element);

    private:
      void reinit_storage();

      H1ShapesetJacobi ref_map_shapeset;
      PrecalcShapesetAssembling ref_map_pss;

      bool is_const;

      double const_jacobian;
      double2x2 const_inv_ref_map;
      double3x2 const_second_ref_map;
      double const_direct_ref_map[2][2];

      double jacobian[H2D_MAX_INTEGRATION_POINTS_COUNT];
      int jacobian_calculated;
      double2x2 inv_ref_map[H2D_MAX_INTEGRATION_POINTS_COUNT];
      int inv_ref_map_calculated;
      double3x2 second_ref_map[H2D_MAX_INTEGRATION_POINTS_COUNT];
      int second_ref_map_calculated;
      double direct_ref_map[2][2][H2D_MAX_INTEGRATION_POINTS_COUNT];
      int direct_ref_map_calculated;
      int inv_ref_order;

      double phys_x[H2D_MAX_INTEGRATION_POINTS_COUNT];
      int phys_x_calculated;
      double phys_y[H2D_MAX_INTEGRATION_POINTS_COUNT];
      int phys_y_calculated;
      double3 tan[H2D_MAX_NUMBER_EDGES][H2D_MAX_INTEGRATION_POINTS_COUNT];
      int tan_calculated[H2D_MAX_NUMBER_EDGES];

      Quad2D* quad_2d;

      void calc_inv_ref_map(int order);

      void calc_const_inv_ref_map();

#ifdef H2D_USE_SECOND_DERIVATIVES
      void calc_second_ref_map(int order);
#endif
      void calc_phys_x(int order);

      void calc_phys_y(int order);

      void calc_tangent(int edge, int eo);

      int calc_inv_ref_order();

      unsigned short indices[70];

      int nc;

      double2* coeffs;

      double2  lin_coeffs[H2D_MAX_NUMBER_EDGES];
    };
  }
}
#endif