forms.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_FORMS_H
#define __H2D_FORMS_H

#include "global.h"
#include "quadrature/quad.h"
#include "function/function.h"
#include "function/exact_solution.h"
#include "mesh/refmap.h"
#include "mesh/traverse.h"

namespace Hermes
{
  namespace Hermes2D
  {
    static const char* ERR_UNDEFINED_NEIGHBORING_ELEMENTS =
      "Neighboring elements are not defined and so are not function traces on their interface. "
      "Did you forget setting H2D_ANY_INNER_EDGE in add_matrix/vector_form?";

#pragma region Geometry

    template<typename T>
    class HERMES_API Geom
    {
    public:
      T x[H2D_MAX_INTEGRATION_POINTS_COUNT];
      T y[H2D_MAX_INTEGRATION_POINTS_COUNT];

      int elem_marker;
    };

    template<typename T>
    class HERMES_API GeomVol : public Geom < T >
    {
    public:
      int id;

      T get_diam_approximation(unsigned char n);
      T get_area(unsigned char n, double* wt);
    };

    template<typename T>
    class HERMES_API GeomSurf : public Geom < T >
    {
    public:
      T nx[H2D_MAX_INTEGRATION_POINTS_COUNT];
      T ny[H2D_MAX_INTEGRATION_POINTS_COUNT];
      T tx[H2D_MAX_INTEGRATION_POINTS_COUNT];
      T ty[H2D_MAX_INTEGRATION_POINTS_COUNT];

      unsigned char isurf;

      int edge_marker;

      bool orientation;
    };

    template<>
    class HERMES_API GeomVol < Hermes::Ord >
    {
    public:
      GeomVol()
      {
        x[0] = y[0] = Hermes::Ord(1);
      }
      Hermes::Ord x[1];
      Hermes::Ord y[1];

      int elem_marker;

      Hermes::Ord get_diam_approximation(unsigned char n) { return Hermes::Ord(1); };
      Hermes::Ord get_area(unsigned char n, double* wt) { return Hermes::Ord(1); };
    };

    template<>
    class HERMES_API GeomSurf < Hermes::Ord >
    {
    public:
      GeomSurf()
      {
        x[0] = y[0] = tx[0] = ty[0] = nx[0] = ny[0] = Hermes::Ord(1);
      }
      Hermes::Ord x[1];
      Hermes::Ord y[1];
      Hermes::Ord tx[1];
      Hermes::Ord ty[1];
      Hermes::Ord nx[1];
      Hermes::Ord ny[1];

      unsigned char np;
    };

    template<typename T>
    class HERMES_API InterfaceGeom
    {
    public:
      InterfaceGeom(GeomSurf<T>* geom, Element* central_el, Element* neighb_el);
      T* x;
      T* y;
      T* nx;
      T* ny;
      T* tx;
      T* ty;
      Element* central_el;
      Element* neighb_el;

      int elem_marker;
      int edge_marker;

      bool orientation;

      unsigned char isurf;

      T get_diam_approximation(unsigned char n);
      T get_area(unsigned char n, double* wt);

    private:
      GeomSurf<T>* wrapped_geom;
      template<typename Scalar> friend class KellyTypeAdapt;
    };

    template<>
    class HERMES_API InterfaceGeom < Hermes::Ord >
    {
    public:
      InterfaceGeom()
      {
        x[0] = y[0] = tx[0] = ty[0] = nx[0] = ny[0] = Hermes::Ord(1);
      }
      Hermes::Ord x[1];
      Hermes::Ord y[1];
      Hermes::Ord tx[1];
      Hermes::Ord ty[1];
      Hermes::Ord nx[1];
      Hermes::Ord ny[1];

      Hermes::Ord get_diam_approximation(unsigned char n) { return Hermes::Ord(1); };
      Hermes::Ord get_area(unsigned char n, double* wt) { return Hermes::Ord(1); };
    };

    HERMES_API GeomVol<double>* init_geom_vol(RefMap *rm, const int order);
    HERMES_API void init_geom_vol_allocated(GeomVol<double>& geom, RefMap *rm, const int order);
    HERMES_API GeomSurf<double>* init_geom_surf(RefMap *rm, unsigned char isurf, int marker, const int order, double3*& tan);
    HERMES_API void init_geom_surf_allocated(GeomSurf<double>& geom, RefMap *rm, unsigned char isurf, int marker, const int order, double3*& tan);
#pragma endregion

#pragma region Func
    template<typename Scalar>
    class HERMES_API Func
    {
    };

    template<>
    class HERMES_API Func < double >
    {
    public:
      Func();

      Func(int np, int nc);
      union
      {
        double val[H2D_MAX_INTEGRATION_POINTS_COUNT];
        double val0[H2D_MAX_INTEGRATION_POINTS_COUNT];
      };

      union
      {
        double dx[H2D_MAX_INTEGRATION_POINTS_COUNT];
        double val1[H2D_MAX_INTEGRATION_POINTS_COUNT];
      };

      union
      {
        double dy[H2D_MAX_INTEGRATION_POINTS_COUNT];
        double curl[H2D_MAX_INTEGRATION_POINTS_COUNT];
      };

      union
      {
        double laplace[H2D_MAX_INTEGRATION_POINTS_COUNT];
        double div[H2D_MAX_INTEGRATION_POINTS_COUNT];
      };

      int np;
      int nc;

      void subtract(Func<double>* func);
      void subtract(double* attribute, double* other_attribute);

      void add(Func<double>* func);
      void add(double* attribute, double* other_attribute);
    };

    template<>
    class HERMES_API Func < std::complex<double> >
    {
    public:
      Func();

      Func(int np, int nc);

      std::complex<double> val[H2D_MAX_INTEGRATION_POINTS_COUNT];
      std::complex<double> val0[H2D_MAX_INTEGRATION_POINTS_COUNT];

      std::complex<double> dx[H2D_MAX_INTEGRATION_POINTS_COUNT];
      std::complex<double> val1[H2D_MAX_INTEGRATION_POINTS_COUNT];

      std::complex<double> dy[H2D_MAX_INTEGRATION_POINTS_COUNT];
      std::complex<double> curl[H2D_MAX_INTEGRATION_POINTS_COUNT];

      std::complex<double> laplace[H2D_MAX_INTEGRATION_POINTS_COUNT];
      std::complex<double> div[H2D_MAX_INTEGRATION_POINTS_COUNT];

      int np;
      int nc;

      void subtract(Func<std::complex<double> >* func);
      void subtract(std::complex<double> * attribute, std::complex<double> * other_attribute);

      void add(Func<std::complex<double> >* func);
      void add(std::complex<double> * attribute, std::complex<double> * other_attribute);
    };

    template<>
    class HERMES_API Func < Ord >
    {
    public:
      Func(const int order);
      const int order;
      Ord val;
      Ord dx;
      Ord dy;
      Ord laplace;

      Ord val0;
      Ord val1;
      Ord curl;
      Ord div;
      void subtract(Func<Ord>* func);
    };

    template<typename T>
    class HERMES_API DiscontinuousFunc : public Func < T >
    {
    public:
      Func<T> *fn_central;
      Func<T> *fn_neighbor;

      T *val;
      T *dx, *dy;
      T *val_neighbor;
      T *dx_neighbor, *dy_neighbor;

      DiscontinuousFunc(Func<T>* fn, bool support_on_neighbor, bool reverse = false);

      DiscontinuousFunc(Func<T>* fn_c, Func<T>* fn_n, bool reverse = false);

      virtual ~DiscontinuousFunc();
      void free();

      using Func<T>::subtract;
      void subtract(const DiscontinuousFunc<T>& func);

      bool reverse_neighbor_side;
      static T zero;
    };

    template<>
    class HERMES_API DiscontinuousFunc<Ord> : public Func < Ord >
    {
    public:
      Func<Ord> *fn_central;
      Func<Ord> *fn_neighbor;

      Ord val;
      Ord dx, dy;
      Ord val_neighbor;
      Ord dx_neighbor, dy_neighbor;

      DiscontinuousFunc(Func<Ord>* fn, bool support_on_neighbor, bool reverse = false);

      DiscontinuousFunc(Func<Ord>* fn_c, Func<Ord>* fn_n, bool reverse = false);

      using Func<Ord>::subtract;
      void subtract(const DiscontinuousFunc<Ord>& func);

      bool reverse_neighbor_side;
      static Ord zero;
    };

    HERMES_API Func<double>* init_fn(PrecalcShapeset *fu, RefMap *rm, const int order);
    template<typename Scalar>
    HERMES_API Func<Scalar>* init_fn(MeshFunction<Scalar>* fu, const int order);

    template<typename Scalar>
    HERMES_API Func<Scalar>* preallocate_fn(pj_pool_t* memoryPool = nullptr);

    HERMES_API void init_fn_preallocated(Func<double>* u, PrecalcShapeset *fu, RefMap *rm, const int order);
    template<typename Scalar>
    HERMES_API void init_fn_preallocated(Func<Scalar>* u, MeshFunction<Scalar>* fu, const int order);
    template<typename Scalar>
    HERMES_API void init_fn_preallocated(Func<Scalar>* u, UExtFunction<Scalar>* fu, Func<Scalar>** ext, Func<Scalar>** u_ext, const int order, Geom<double>* geometry, ElementMode2D mode);

    template<typename Scalar>
    HERMES_API Func<Scalar>* init_zero_fn(ElementMode2D mode, int order, Quad2D* quad_2d = nullptr, int nc = 1);

    template<typename Scalar>
    HERMES_API Func<Scalar>* init_fn(UExtFunction<Scalar>* fu, Func<Scalar>** ext, Func<Scalar>** u_ext, const int order, Geom<double>* geometry, ElementMode2D mode);
#pragma endregion
  }
}
#endif