space.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_SPACE_H
#define __H2D_SPACE_H

#include "../shapeset/shapeset.h"
#include "asmlist.h"
#include "../boundary_conditions/essential_boundary_conditions.h"
#include "../mesh/traverse.h"
#include "../quadrature/quad_all.h"
#include "algebra/dense_matrix_operations.h"

using namespace Hermes::Algebra::DenseMatrixOperations;

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

      SpaceSharedPtr(const SpaceSharedPtr<Scalar>& other);

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

    namespace Mixins
    {
      template<typename Scalar>
      class HERMES_API SettableSpaces
      {
      public:
        virtual void set_spaces(std::vector<SpaceSharedPtr<Scalar> > spaces) = 0;
        virtual void set_space(SpaceSharedPtr<Scalar> space);
        virtual std::vector<SpaceSharedPtr<Scalar> > get_spaces();
        virtual SpaceSharedPtr<Scalar> get_space(int n);
      };
    }

    template<typename Scalar> class Adapt;
    template<typename Scalar> class DiscreteProblem;
    template<typename Scalar> class DiscreteProblemDGAssembler;
    template<typename Scalar> class DiscreteProblemThreadAssembler;
    template<typename Scalar> class DiscreteProblemIntegrationOrderCalculator;
    namespace Views
    {
      template<typename Scalar> class BaseView;
      template<typename Scalar> class VectorBaseView;
      class Orderizer;
      class OrderView;
    };
    class Shapeset;

    template<typename Scalar> class L2Space;
    template<typename Scalar> class H1Space;
    template<typename Scalar> class HcurlSpace;
    template<typename Scalar> class HdivSpace;

    template<typename Scalar>
    class HERMES_API Space : public Hermes::Mixins::Loggable, public Hermes::Mixins::StateQueryable, public Hermes::Hermes2D::Mixins::XMLParsing
    {
    public:
      Space();
      Space(MeshSharedPtr mesh, Shapeset* shapeset, EssentialBCs<Scalar>* essential_bcs);

      void init();

      virtual void init(Shapeset* shapeset, int p_init, bool assign_dofs_init = true) = 0;

      virtual ~Space();

      void free();
      void free_bc_data();

#pragma region Getters
      inline int get_element_order(int id) const {
        return edata[id].order;
      }

      int get_num_dofs() const;

      static int get_num_dofs(std::vector<SpaceSharedPtr<Scalar> > spaces);

      static int get_num_dofs(SpaceSharedPtr<Scalar> space);

      MeshSharedPtr get_mesh() const;

      virtual SpaceType get_type() const = 0;

      int get_vertex_functions_count();
      int get_edge_functions_count();
      int get_bubble_functions_count();

      int get_seq() const;

      EssentialBCs<Scalar>* get_essential_bcs() const;

      virtual void get_element_assembly_list(Element* e, AsmList<Scalar>* al) const;

      virtual int get_edge_order(Element* e, int edge) const;

      int get_max_dof() const;

      void get_boundary_assembly_list(Element* e, int surf_num, AsmList<Scalar>* al) const;

      Shapeset* get_shapeset() const;
#pragma endregion

#pragma region Setters
      virtual void set_shapeset(Shapeset* shapeset, bool clone = false) = 0;

      void set_mesh(MeshSharedPtr mesh);

      void set_mesh_seq(int seq);

      void set_essential_bcs(EssentialBCs<Scalar>* essential_bcs);
#pragma endregion

#pragma region Order setting
      virtual void set_element_order(int id, int order, int order_v = -1);

      void set_uniform_order(int order, std::string marker = HERMES_ANY);

      void adjust_element_order(int order_change, int min_order);

      void adjust_element_order(int horizontal_order_change, int vertical_order_change, unsigned int horizontal_min_order, unsigned int vertical_min_order);

      void set_uniform_order_internal(int order, int marker);

      virtual int assign_dofs(int first_dof = 0);

      static int assign_dofs(std::vector<SpaceSharedPtr<Scalar> > spaces);
#pragma endregion

#pragma region Mesh handling
      void unrefine_all_mesh_elements(bool keep_initial_refinements = true);

      static void unrefine_all_mesh_elements(std::vector<SpaceSharedPtr<Scalar> > spaces, bool keep_initial_refinements = true);
#pragma endregion

#pragma region Boundary conditions
      virtual Scalar* get_bc_projection(SurfPos* surf_pos, int order, EssentialBoundaryCondition<Scalar> *bc) = 0;

      void update_essential_bc_values();

      static void update_essential_bc_values(std::vector<SpaceSharedPtr<Scalar> > spaces, double time);

      static void update_essential_bc_values(SpaceSharedPtr<Scalar> space, double time);
#pragma endregion

#pragma region Save & Load
      void save(const char *filename) const;
#ifdef WITH_BSON
      void save_bson(const char* filename) const;
#endif

      static SpaceSharedPtr<Scalar> load(const char *filename, MeshSharedPtr mesh, bool validate = false, EssentialBCs<Scalar>* essential_bcs = nullptr, Shapeset* shapeset = nullptr);
      void load(const char *filename);

#ifdef WITH_BSON
      static SpaceSharedPtr<Scalar> load_bson(const char *filename, MeshSharedPtr mesh, EssentialBCs<Scalar>* essential_bcs = nullptr, Shapeset* shapeset = nullptr);
      void load_bson(const char *filename);
#endif
#pragma endregion

      virtual void copy(SpaceSharedPtr<Scalar> space, MeshSharedPtr new_mesh);

      class HERMES_API ReferenceSpaceCreator
      {
      public:
        ReferenceSpaceCreator(unsigned int order_increase = 1);
        ReferenceSpaceCreator(SpaceSharedPtr<Scalar> coarse_space, MeshSharedPtr ref_mesh, unsigned int order_increase = 1);

        virtual void handle_orders(SpaceSharedPtr<Scalar> ref_space);

        virtual SpaceSharedPtr<Scalar> create_ref_space(bool assign_dofs = true);
        virtual SpaceSharedPtr<Scalar> create_ref_space(SpaceSharedPtr<Scalar> coarse_space, MeshSharedPtr ref_mesh, bool assign_dofs = true);

      private:
        SpaceSharedPtr<Scalar> init_construction_l2();
        SpaceSharedPtr<Scalar> init_construction_h1();
        SpaceSharedPtr<Scalar> init_construction_hcurl();
        SpaceSharedPtr<Scalar> init_construction_hdiv();

        virtual void finish_construction(SpaceSharedPtr<Scalar> ref_space);

        SpaceSharedPtr<Scalar> coarse_space;
        MeshSharedPtr ref_mesh;
        unsigned int order_increase;
      };

    protected:

      virtual void set_element_order_internal(int id, int order, int order_v = -1);

      void distribute_orders(MeshSharedPtr mesh, int* parents);

      bool is_up_to_date() const;

      static Node* get_mid_edge_vertex_node(Element* e, unsigned char i, unsigned char j);

      virtual bool isOkay() const;
      inline std::string getClassName() const { return "Space"; }

      int ndof;

      static const int H2D_UNASSIGNED_DOF = -2;
      static const int H2D_CONSTRAINED_DOF = -1;

      Shapeset* shapeset;
      bool own_shapeset;

      EssentialBCs<Scalar>* essential_bcs;

      MeshSharedPtr mesh;

      int vertex_functions_count, edge_functions_count, bubble_functions_count;

      int first_dof, next_dof;

      unsigned int seq;
      unsigned int seq_assigned;
      int mesh_seq;

      struct BaseComponent
      {
        int dof;
        Scalar coef;
      };

      union NodeData
      {
        struct // regular node
        {
          int dof;
          union {
            Scalar* edge_bc_proj;
            Scalar* vertex_bc_coef;
          };
          short n;
        };
        struct // constrained vertex node
        {
          BaseComponent* baselist;
          int ncomponents;
        };
        struct // constrained edge node
        {
          Node* base;
          int part;
        };
        NodeData() : dof(0), edge_bc_proj(nullptr) {}
      };

      class ElementData
      {
      public:
        ElementData() : changed_in_last_adaptation(true) {};
        short order;
        int bdof;
        short n;
        bool changed_in_last_adaptation;
      };

      NodeData* ndata;
      int nsize;
      ElementData* edata;
      int esize;

      virtual int get_edge_order_internal(Node* en) const;

      void unrefine_all_mesh_elements_internal(bool keep_initial_refinements, bool only_unrefine_space_data);

      virtual void resize_tables();

      void update_orders_recurrent(Element* e, int order);

      virtual void reset_dof_assignment();
      virtual void assign_vertex_dofs() = 0;
      virtual void assign_edge_dofs() = 0;
      virtual void assign_bubble_dofs() = 0;

      virtual void get_vertex_assembly_list(Element* e, int iv, AsmList<Scalar>* al) const = 0;
      virtual void get_boundary_assembly_list_internal(Element* e, int surf_num, AsmList<Scalar>* al) const = 0;
      virtual void get_bubble_assembly_list(Element* e, AsmList<Scalar>* al) const;

      double** proj_mat;
      double*  chol_p;

      std::vector<Scalar*> bc_data_projections;
      std::vector<typename Space<Scalar>::BaseComponent*> bc_data_base_components;

      void precalculate_projection_matrix(int nv, double**& mat, double*& p);
      void update_edge_bc(Element* e, SurfPos* surf_pos);

      virtual void update_constraints();

      virtual void post_assign();

      static SpaceSharedPtr<Scalar> init_empty_space(SpaceType spaceType, MeshSharedPtr mesh, Shapeset* shapeset);

      struct EdgeInfo
      {
        Node* node;
        int part;
        int ori;
        double lo, hi;
      };

      template<typename T> friend class OGProjection;
      template<typename T> friend class NewtonSolver;
      template<typename T> friend class PicardSolver;
      template<typename T> friend class LinearSolver;
      template<typename T> friend class OGProjectionNOX;
      template<typename T> friend class Solution;
      template<typename T> friend class RungeKutta;
      template<typename T> friend class ExactSolution;
      template<typename T> friend class NeighborSearch;
      template<typename T> friend class ExactSolutionScalar;
      template<typename T> friend class ExactSolutionVector;
      template<typename T> friend class Views::BaseView;
      friend class Views::Orderizer;
      friend class Views::OrderView;
      template<typename T> friend class Views::VectorBaseView;
      friend class Adapt < Scalar > ;
      friend class DiscreteProblem < Scalar > ;
      friend class DiscreteProblemDGAssembler < Scalar > ;
      friend class DiscreteProblemThreadAssembler < Scalar > ;
      friend class DiscreteProblemIntegrationOrderCalculator < Scalar > ;
    };
  }
}
#endif