curved.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_CURVED_H
#define __H2D_CURVED_H

#include "../global.h"
#include "../shapeset/shapeset_common.h"
#include "../shapeset/precalc.h"

namespace Hermes
{
  namespace Hermes2D
  {
    class Element;
    class RefMap;
    class Quad1DStd;
    class Quad2DStd;
    struct Trf;

    enum CurvType
    {
      ArcType,
      NurbsType
    };

    class HERMES_API Curve
    {
    public:
      Curve(CurvType type);
      virtual ~Curve();
      CurvType type;
    };

    class HERMES_API Arc : public Curve
    {
    public:
      Arc();
      Arc(double angle);
      Arc(const Arc* other);

      double angle;

      static const unsigned char degree = 2;

      static const unsigned char np = 3;
      // there are 6 knots: {0, 0, 0, 1, 1, 1}
      static const unsigned char nk = 6;
      double kv[6];
      double3 pt[3];
    };

    class HERMES_API Nurbs : public Curve
    {
    public:
      Nurbs();
      Nurbs(const Nurbs* other);
      ~Nurbs();

      unsigned char degree;
      unsigned char np;
      double3* pt;
      unsigned char nk;
      double* kv;
    };

    class HERMES_API CurvMap
    {
    public:
      CurvMap();
      CurvMap(const CurvMap* cm);
      ~CurvMap();
      void free();

      bool toplevel;

      // if toplevel=true, this structure belongs to a base mesh element
      // and the array 'nurbs' points to (up to four) NURBS curved edges
      Curve* curves[H2D_MAX_NUMBER_EDGES];

      // if toplevel=false, this structure belongs to a refined element
      // and 'parent' points to the base mesh element CurvMap structure;
      Element* parent;
      uint64_t sub_idx;

      unsigned short order;

      void update_refmap_coeffs(Element* e);

    private:
      PrecalcShapesetAssembling ref_map_pss;

      Trf* ctm;

      unsigned short nc;
      double2* coeffs;

      void get_mid_edge_points(Element* e, double2* pt, unsigned short n);

      static double nurbs_basis_fn(unsigned short i, unsigned short k, double t, double* knot);

      static void nurbs_edge(Element* e, Curve* curve, int edge, double t, double& x,
        double& y);

      static const double2 ref_vert[2][H2D_MAX_NUMBER_VERTICES];

      static void nurbs_edge_0(Element* e, Curve* nurbs, unsigned short edge, double t, double& x, double& y, double& n_x, double& n_y, double& t_x, double& t_y);

      static void calc_ref_map_tri(Element* e, Curve** nurbs, double xi_1, double xi_2, double& x, double& y);
      static void calc_ref_map_quad(Element* e, Curve** nurbs, double xi_1, double xi_2,
        double& x, double& y);

      static void calc_ref_map(Element* e, Curve** nurbs, double xi_1, double xi_2, double2& f);

      void edge_coord(Element* e, unsigned short edge, double t, double2& x) const;
      void calc_edge_projection(Element* e, unsigned short edge, Curve** nurbs, unsigned short order, double2* proj) const;

      void old_projection(Element* e, unsigned short order, double2* proj, double* old[2]);
      void calc_bubble_projection(Element* e, Curve** nurbs, unsigned short order, double2* proj);

      static CurvMap* create_son_curv_map(Element* e, int son);

      template<typename T> friend class Space;
      template<typename T> friend class H1Space;
      template<typename T> friend class L2Space;
      template<typename T> friend class HcurlSpace;
      template<typename T> friend class HdivSpace;
      template<typename T> friend class DiscreteProblem;
      template<typename T> friend class Adapt;
      template<typename T> friend class KellyTypeAdapt;
      friend class RefMap;
      friend class Mesh;
      friend class MeshReader;
      friend class MeshReaderH2D;
      friend class MeshReaderH2DXML;
      friend class MeshReaderH2DBSON;
    };

    class CurvMapStatic
    {
    public:
      CurvMapStatic();
      ~CurvMapStatic();

      void precalculate_cholesky_projection_matrix_edge();
      double** calculate_bubble_projection_matrix(short* indices, ElementMode2D mode);
      void precalculate_cholesky_projection_matrices_bubble();

      double** edge_proj_matrix;
      unsigned short edge_proj_matrix_size;
      double** bubble_proj_matrix_tri;
      double** bubble_proj_matrix_quad;

      double* edge_p;
      double* bubble_tri_p;
      unsigned short tri_bubble_np;
      double* bubble_quad_p;
      unsigned short quad_bubble_np;
    };

    extern CurvMapStatic curvMapStatic;
  }
}
#endif