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/solution.h"
#include <string>

namespace Hermes
{
  namespace Hermes2D
  {
    enum GeomType
    {
      HERMES_PLANAR = 0,         // Planar problem.
      HERMES_AXISYM_X = 1,       // Axisymmetric problem where x-axis is the axis of symmetry.
      HERMES_AXISYM_Y = 2        // Axisymmetric problem where y-axis is the axis of symmetry.
    };

    class RefMap;
    template<typename Scalar> class DiscreteProblem;
    template<typename Scalar> class DiscreteProblemLinear;
    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 Geom;

    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;

    enum SymFlag
    {
      HERMES_ANTISYM = -1,
      HERMES_NONSYM = 0,
      HERMES_SYM = 1
    };

    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);

      ~WeakForm();

      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);

      unsigned int 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;

      Hermes::vector<Form<Scalar> *> get_forms() const;

      Hermes::vector<MatrixFormVol<Scalar> *> get_mfvol() const;
      Hermes::vector<MatrixFormSurf<Scalar> *> get_mfsurf() const;
      Hermes::vector<MatrixFormDG<Scalar> *> get_mfDG() const;

      Hermes::vector<VectorFormVol<Scalar> *> get_vfvol() const;
      Hermes::vector<VectorFormSurf<Scalar> *> get_vfsurf() const;
      Hermes::vector<VectorFormDG<Scalar> *> get_vfDG() const;

      void delete_all();

      void set_ext(MeshFunction<Scalar>* ext);
      void set_ext(Hermes::vector<MeshFunction<Scalar>*> ext);
      Hermes::vector<MeshFunction<Scalar>*> get_ext() const;

      virtual WeakForm* clone() const;

    protected:
      Hermes::vector<MeshFunction<Scalar>*> ext;

      double current_time;
      double current_time_step;

      unsigned int neq;

      bool is_matfree;

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

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

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

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

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

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

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

      bool** get_blocks(bool force_diagonal_blocks) const;

      friend class DiscreteProblem<Scalar>;
      friend class DiscreteProblemLinear<Scalar>;
      friend class RungeKutta<Scalar>;
      friend class OGProjection<Scalar>;
      friend class Hermes::Preconditioners::Precond<Scalar>;

      bool warned_nonOverride;

      // internal.
    private:
      void free_ext();
      void cloneMembers(const WeakForm<Scalar>* otherWf);
    };

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

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

      void set_ext(MeshFunction<Scalar>* ext);
      void set_ext(Hermes::vector<MeshFunction<Scalar>*> ext);
      Hermes::vector<MeshFunction<Scalar>*> get_ext() const;

    protected:
      inline void set_weakform(WeakForm<Scalar>* wf) { this->wf = wf; }

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

      int u_ext_offset;

      Hermes::vector<MeshFunction<Scalar>*> ext;

      void set_current_stage_time(double time);

      double get_current_stage_time() const;

      double scaling_factor;

      WeakForm<Scalar>* wf;
      double stage_time;
      void setScalingFactor(double scalingFactor);
      void set_uExtOffset(int u_ext_offset);
      friend class WeakForm<Scalar>;
      friend class RungeKutta<Scalar>;
      friend class DiscreteProblem<Scalar>;
      friend class DiscreteProblemLinear<Scalar>;
    };

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

      virtual ~MatrixForm();

      unsigned int i;
      unsigned int j;
      unsigned int previous_iteration_space_index;

      SymFlag sym;

      typedef Scalar valueFunction(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *u, Func<double> *v,
        Geom<double> *e, Func<Scalar> **ext) const;

      virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *u, Func<double> *v,
        Geom<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,
        Geom<Hermes::Ord> *e, Func<Ord> **ext) const;
    
    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 MatrixFormVol* clone() const;
    };

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

      virtual ~MatrixFormSurf();

      virtual MatrixFormSurf* clone() const;
    };

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

      virtual ~MatrixFormDG();

      virtual MatrixFormDG* clone() const;
    };

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

      virtual ~VectorForm();

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

      virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[], Func<Hermes::Ord> *v, Geom<Hermes::Ord> *e,
        Func<Ord> **ext) const;
      unsigned int i;

    protected:
      friend class DiscreteProblem<Scalar>;
    };

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

      virtual ~VectorFormVol();

      virtual VectorFormVol* clone() const;
    };

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

      virtual ~VectorFormSurf();

      virtual VectorFormSurf* clone() const;
    };

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

      virtual ~VectorFormDG();

      virtual VectorFormDG* clone() const;
    };
  }
}
#endif