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/solution.h"
#include "mesh/refmap.h"
#include "mesh/traverse.h"
#include <complex>

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?";

    template<typename Scalar> class OGProjection;

    template<typename T>
    class HERMES_API Func
    {
    public:

      Func(int num_gip, int num_comps);

      T *val;            
      T *dx, *dy;        
      T *laplace;        
      T *val0, *val1;    
      T *dx0, *dx1;      
      T *dy0, *dy1;      
      T *curl;           
      T *div;            

      virtual T& get_val_central(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_val_neighbor(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_dx_central(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_dx_neighbor(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_dy_central(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS);  return * new T; }
      virtual T& get_dy_neighbor(int k) const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_laplace_central(int k) { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }
      virtual T& get_laplace_neighbor(int k) { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return * new T; }

      virtual void free_ord();

      virtual void free_fn();

      virtual ~Func() { };

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

      int get_num_gip() const;

    protected:
      const int num_gip; 
      const int nc;      


      void subtract(T* attribute, T* other_attribute);


      void add(Func<T>* func);

      friend Func<Hermes::Ord>* init_fn_ord(const int order);
      friend Func<double>* init_fn(PrecalcShapeset *fu, RefMap *rm, const int order);
      template<typename Scalar> friend Func<Scalar>* init_fn(MeshFunction<Scalar>*fu, const int order);
      template<typename Scalar> friend Func<Scalar>* init_fn(Solution<Scalar>*fu, const int order);

      template<typename Scalar> friend class DiscontinuousFunc;
      template<typename Scalar> friend class Adapt;
      template<typename Scalar> friend class KellyTypeAdapt;
      template<typename Scalar> friend class DiscreteProblem;
      template<typename Scalar> friend class DiscreteProblemLinear;
      template<typename Scalar> friend class BasicKellyAdapt;
      template<typename Scalar> friend class Hermes::Hermes2D::OGProjection;
      friend class ErrorEstimatorFormKelly;
    };


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

      virtual T& get_val_central(int k) const;
      virtual T& get_val_neighbor(int k) const;
      virtual T& get_dx_central(int k) const;
      virtual T& get_dx_neighbor(int k) const;
      virtual T& get_dy_central(int k) const;
      virtual T& get_dy_neighbor(int k) const;
      virtual T& get_laplace_central(int k);
      virtual T& get_laplace_neighbor(int k);

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

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

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

      virtual void free_fn();

      virtual void free_ord();

      bool reverse_neighbor_side; 

      static T zero;              
      template<typename Scalar> friend class DiscreteProblem;
      template<typename Scalar> friend class DiscreteProblemLinear;
      template<typename Scalar> friend class KellyTypeAdapt;
      template<typename Scalar> friend class Adapt;
      template<typename Scalar> friend class NeighborSearch;
    };

    template<typename T>
    class HERMES_API Geom
    {
    public:
      Geom();

      T diam;           
      T area;           
      T *x, *y;         
      T *nx, *ny;       


      T *tx, *ty;       
      int id;           
      int isurf;        

      virtual int get_neighbor_marker() const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return -1; }
      virtual int get_neighbor_id()     const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return -1; }
      virtual T   get_neighbor_diam()   const { throw Hermes::Exceptions::Exception(ERR_UNDEFINED_NEIGHBORING_ELEMENTS); return  T(); }

      virtual ~Geom() {};

      virtual void free();
      virtual void free_ord() {};
      int elem_marker;       
      int edge_marker;       

    protected:
      int orientation;  



      friend Geom<Hermes::Ord>* init_geom_ord();
      friend Geom<double>* init_geom_vol(RefMap *rm, const int order);
      friend Geom<double>* init_geom_surf(RefMap *rm, int isurf, int marker, const int order, double3*& tan);

      template<typename Scalar> friend class DiscreteProblem;
      template<typename Scalar> friend class DiscreteProblemLinear;
      template<typename Scalar> friend class InterfaceGeom;
      template<typename Scalar> friend class KellyTypeAdapt;
      template<typename Scalar> friend class BasicKellyAdapt;
      friend class ErrorEstimatorFormKelly;
      template<typename Scalar> friend class Adapt;
    };

    template<typename T>
    class HERMES_API InterfaceGeom : public Geom<T>
    {
    public:
      int neighb_id;
      T   neighb_diam;
      int get_neighbor_marker() const;
      int get_neighbor_id()  const;
      T get_neighbor_diam() const;

    private:
      InterfaceGeom(Geom<T>* geom, int n_marker, int n_id, T n_diam);
      Geom<T>* wrapped_geom;

      virtual void free();
      virtual void free_ord();

      int neighb_marker;
      template<typename Scalar> friend class DiscreteProblem;
      template<typename Scalar> friend class DiscreteProblemLinear;
      template<typename Scalar> friend class KellyTypeAdapt;
      template<typename Scalar> friend class Adapt;
    };

    HERMES_API Geom<Hermes::Ord>* init_geom_ord();
    HERMES_API Geom<double>* init_geom_vol(RefMap *rm, const int order);
    HERMES_API Geom<double>* init_geom_surf(RefMap *rm, int isurf, int marker, const int order, double3*& tan);

    HERMES_API Func<Hermes::Ord>* init_fn_ord(const int order);
    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>* init_fn(Solution<Scalar>*fu, const int order);
  }
}
#endif