weakform.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 Licenserix
// along with Hermes2D.  If not, see <http://www.gnu.org/licenses/>.

#ifndef __H2D_WEAKFORM_H
#define __H2D_WEAKFORM_H

#include "../function/exact_solution.h"

namespace Hermes
{
  namespace Hermes2D
  {
#pragma region forward-declarations
    class RefMap;
    template<typename Scalar> class DiscreteProblem;
    template<typename Scalar> class DiscreteProblemSelectiveAssembler;
    template<typename Scalar> class DiscreteProblemIntegrationOrderCalculator;
    template<typename Scalar> class RungeKutta;
    template<typename Scalar> class Space;
    template<typename Scalar> class MeshFunction;
    struct SurfPos;

    class Element;
    class Shapeset;
    template<typename T> class Func;
    template<typename T> class DiscontinuousFunc;
    template<typename T> class InterfaceGeom;
    template<typename T> class GeomVol;
    template<typename T> class GeomSurf;

    template<typename Scalar> class Form;
    template<typename Scalar> class OGProjection;
    template<typename Scalar> class MatrixFormVol;
    template<typename Scalar> class VectorFormVol;
    template<typename Scalar> class MatrixFormSurf;
    template<typename Scalar> class VectorFormSurf;
    template<typename Scalar> class MatrixFormDG;
    template<typename Scalar> class VectorFormDG;
#pragma endregion

    template<typename Scalar>
    class HERMES_API WeakFormSharedPtr : public std::tr1::shared_ptr < Hermes::Hermes2D::WeakForm<Scalar> >
    {
    public:
      WeakFormSharedPtr(Hermes::Hermes2D::WeakForm<Scalar>* ptr = nullptr);

      WeakFormSharedPtr(const WeakFormSharedPtr<Scalar>& other);

      void operator=(const WeakFormSharedPtr<Scalar>& other);
    };

    template<typename Scalar>
    class HERMES_API WeakForm : public Hermes::Mixins::IntegrableWithGlobalOrder, public Hermes::Mixins::Loggable
    {
    public:
      WeakForm(unsigned int neq = 1, bool mat_free = false);

      virtual ~WeakForm();

      enum FormIntegrationDimension
      {
        FormVol = 0,
        FormSurf = 1
      };

      enum FormEquationSide
      {
        MatrixForm = 0,
        VectorForm = 1
      };

      void add_matrix_form(MatrixFormVol<Scalar>* mfv);

      void add_matrix_form_surf(MatrixFormSurf<Scalar>* mfs);

      void add_matrix_form_DG(MatrixFormDG<Scalar>* mfDG);

      void add_vector_form(VectorFormVol<Scalar>* vfv);

      void add_vector_form_surf(VectorFormSurf<Scalar>* vfs);

      void add_vector_form_DG(VectorFormDG<Scalar>* vfDG);

      virtual void set_active_state(Element** e);

      virtual void set_active_edge_state(Element** e, unsigned char isurf);

      virtual void set_active_DG_state(Element** e, unsigned char isurf);

      inline unsigned char get_neq() const { return neq; }

      bool is_matrix_free() const { return is_matfree; }

      void set_current_time(double time);

      void set_current_time_step(double time_step);

      virtual double get_current_time() const;

      virtual double get_current_time_step() const;

      void set_u_ext_fn(UExtFunctionSharedPtr<Scalar> ext);

      void set_ext(MeshFunctionSharedPtr<Scalar> ext);

      void set_u_ext_fn(std::vector<UExtFunctionSharedPtr<Scalar> > ext);

      void set_ext(std::vector<MeshFunctionSharedPtr<Scalar> > ext);

      std::vector<MeshFunctionSharedPtr<Scalar> > get_ext() const;

      virtual WeakForm* clone() const;

      // Checks presence of DG forms.
      bool is_DG() const;

      std::vector<Form<Scalar> *> get_forms() const;
      std::vector<MatrixFormVol<Scalar> *> get_mfvol() const;
      std::vector<MatrixFormSurf<Scalar> *> get_mfsurf() const;
      std::vector<MatrixFormDG<Scalar> *> get_mfDG() const;
      std::vector<VectorFormVol<Scalar> *> get_vfvol() const;
      std::vector<VectorFormSurf<Scalar> *> get_vfsurf() const;
      std::vector<VectorFormDG<Scalar> *> get_vfDG() const;

      void delete_all();

    protected:
      std::vector<MeshFunctionSharedPtr<Scalar> > ext;
      std::vector<UExtFunctionSharedPtr<Scalar> > u_ext_fn;

      double current_time;
      double current_time_step;

      unsigned int neq;

      unsigned int original_neq;

      bool is_matfree;

      std::vector<Form<Scalar> *> forms;

      std::vector<MatrixFormVol<Scalar> *> mfvol;

      std::vector<MatrixFormSurf<Scalar> *> mfsurf;

      std::vector<MatrixFormDG<Scalar> *> mfDG;

      std::vector<VectorFormVol<Scalar> *> vfvol;

      std::vector<VectorFormSurf<Scalar> *> vfsurf;

      std::vector<VectorFormDG<Scalar> *> vfDG;

      bool** get_blocks(bool force_diagonal_blocks) const;

      friend class DiscreteProblem < Scalar > ;
      friend class Form < Scalar > ;
      friend class DiscreteProblemDGAssembler < Scalar > ;
      friend class DiscreteProblemThreadAssembler < Scalar > ;
      friend class DiscreteProblemIntegrationOrderCalculator < Scalar > ;
      friend class DiscreteProblemSelectiveAssembler < Scalar > ;
      friend class RungeKutta < Scalar > ;
      friend class OGProjection < Scalar > ;
      friend class Hermes::Preconditioners::Precond < Scalar > ;

      // Internal.
      virtual void cloneMembers(const WeakFormSharedPtr<Scalar>& other_wf);
      // Internal.
      void cloneMemberExtFunctions(std::vector<MeshFunctionSharedPtr<Scalar> > source_ext, std::vector<MeshFunctionSharedPtr<Scalar> >& cloned_ext);

      // Internal - processes markers, translates from strings to ints.
      template<typename FormType>
      void processFormMarkers(const std::vector<SpaceSharedPtr<Scalar> > spaces, bool surface, std::vector<FormType> forms_to_process);
      void processFormMarkers(const std::vector<SpaceSharedPtr<Scalar> > spaces);

    private:
      void free_ext();
    };

    template<typename Scalar>
    class HERMES_API Form
    {
    public:
      Form(int i = 0);
      virtual ~Form();

      void set_area(std::string area);
      void set_areas(std::vector<std::string> areas);
      std::vector<std::string> getAreas() const;

      void set_ext(MeshFunctionSharedPtr<Scalar> ext);
      void set_u_ext_fn(UExtFunctionSharedPtr<Scalar> ext);

      void set_ext(std::vector<MeshFunctionSharedPtr<Scalar> > ext);
      void set_u_ext_fn(std::vector<UExtFunctionSharedPtr<Scalar> > ext);
      std::vector<MeshFunctionSharedPtr<Scalar> > get_ext() const;

      void setScalingFactor(double scalingFactor);

      unsigned int i;

    protected:
      void set_weakform(WeakForm<Scalar>* wf);

      std::vector<std::string> areas;

      std::vector<int> areas_internal;

      bool assembleEverywhere;

      int u_ext_offset;

      unsigned int previous_iteration_space_index;

      std::vector<MeshFunctionSharedPtr<Scalar> > ext;
      std::vector<UExtFunctionSharedPtr<Scalar> > u_ext_fn;

      double get_current_stage_time() const;

      WeakForm<Scalar>* wf;
    private:
      double stage_time;
      void set_uExtOffset(int u_ext_offset);
      double scaling_factor;
      void set_current_stage_time(double time);
      void copy_base(Form<Scalar>* other_form);
      friend class WeakForm < Scalar > ;
      friend class RungeKutta < Scalar > ;
      friend class DiscreteProblem < Scalar > ;
      friend class DiscreteProblemDGAssembler < Scalar > ;
      friend class DiscreteProblemIntegrationOrderCalculator < Scalar > ;
      friend class DiscreteProblemSelectiveAssembler < Scalar > ;
      friend class DiscreteProblemThreadAssembler < Scalar > ;
    };

    template<typename Scalar>
    class HERMES_API MatrixForm : public Form < Scalar >
    {
    public:
      MatrixForm(unsigned int i, unsigned int j);

      virtual ~MatrixForm();

      unsigned int j;

      SymFlag sym;

    protected:
      friend class DiscreteProblem < Scalar > ;
    };

    template<typename Scalar>
    class HERMES_API MatrixFormVol : public MatrixForm < Scalar >
    {
    public:
      MatrixFormVol(unsigned int i, unsigned int j);

      void setSymFlag(SymFlag sym);
      SymFlag getSymFlag() const;

      virtual ~MatrixFormVol();

      virtual Scalar value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *u, Func<double> *v,
        GeomVol<double> *e, Func<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
        GeomVol<Hermes::Ord> *e, Func<Ord> **ext) const;

      virtual MatrixFormVol* clone() const;
    };

    template<typename Scalar>
    class HERMES_API MatrixFormSurf : public MatrixForm < Scalar >
    {
    public:
      MatrixFormSurf(unsigned int i, unsigned int j);

      virtual ~MatrixFormSurf();

      virtual Scalar value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *u, Func<double> *v,
        GeomSurf<double> *e, Func<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
        GeomSurf<Hermes::Ord> *e, Func<Ord> **ext) const;

      virtual MatrixFormSurf* clone() const;
    };

    template<typename Scalar>
    class HERMES_API MatrixFormDG : public Form < Scalar >
    {
    public:
      MatrixFormDG(unsigned int i, unsigned int j);

      virtual ~MatrixFormDG();

      unsigned int j;

      virtual Scalar value(int n, double *wt, DiscontinuousFunc<Scalar> **u_ext, DiscontinuousFunc<double> *u, DiscontinuousFunc<double> *v,
        InterfaceGeom<double> *e, DiscontinuousFunc<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, DiscontinuousFunc<Hermes::Ord> **u_ext, DiscontinuousFunc<Hermes::Ord> *u, DiscontinuousFunc<Hermes::Ord> *v,
        InterfaceGeom<Hermes::Ord> *e, DiscontinuousFunc<Ord> **ext) const;

      virtual MatrixFormDG* clone() const;
    protected:
      friend class DiscreteProblem < Scalar > ;
    };

    template<typename Scalar>
    class VectorForm : public Form < Scalar >
    {
    public:
      VectorForm(unsigned int i);

      virtual ~VectorForm();

    protected:
      friend class DiscreteProblem < Scalar > ;
    };

    template<typename Scalar>
    class VectorFormVol : public VectorForm < Scalar >
    {
    public:
      VectorFormVol(unsigned int i);

      virtual ~VectorFormVol();

      virtual Scalar value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *v,
        GeomVol<double> *e, Func<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v, GeomVol<Hermes::Ord> *e,
        Func<Ord> **ext) const;

      virtual VectorFormVol* clone() const;
    };

    template<typename Scalar>
    class VectorFormSurf : public VectorForm < Scalar >
    {
    public:
      VectorFormSurf(unsigned int i);

      virtual ~VectorFormSurf();

      virtual Scalar value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *v,
        GeomSurf<double> *e, Func<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v, GeomSurf<Hermes::Ord> *e,
        Func<Ord> **ext) const;

      virtual VectorFormSurf* clone() const;
    };

    template<typename Scalar>
    class VectorFormDG : public Form < Scalar >
    {
    public:
      VectorFormDG(unsigned int i);

      virtual ~VectorFormDG();

      virtual Scalar value(int n, double *wt, DiscontinuousFunc<Scalar> **u_ext, Func<double> *v,
        InterfaceGeom<double> *e, DiscontinuousFunc<Scalar> **ext) const;

      virtual Hermes::Ord ord(int n, double *wt, DiscontinuousFunc<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v, InterfaceGeom<Hermes::Ord> *e,
        DiscontinuousFunc<Ord> **ext) const;

      virtual VectorFormDG* clone() const;
    protected:
      friend class DiscreteProblem < Scalar > ;
    };
  }
}
#endif