solution.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_SOLUTION_H
#define __H2D_SOLUTION_H

#include "../function/mesh_function.h"
#include "../space/space.h"
#include "../mesh/refmap.h"
#include "exceptions.h"

namespace Hermes
{
  namespace Hermes2D
  {
    class Quad2DCheb;

    enum SolutionType {
      HERMES_UNDEF = -1,
      HERMES_SLN = 0,
      HERMES_EXACT = 1
    };

    template<typename Scalar>
    class HERMES_API Solution : public MeshFunction<Scalar>, public Hermes2D::Mixins::XMLParsing
    {
    public:
      Solution();
      Solution(const Mesh *mesh);
      Solution (Space<Scalar>* s, Vector<Scalar>* coeff_vec);
      Solution (Space<Scalar>* s, Scalar* coeff_vec);
      virtual ~Solution();

      inline std::string getClassName() const { return "Solution"; }

      void assign(Solution* sln);
      inline Solution& operator = (Solution& sln) { assign(&sln); return *this; }

      virtual void copy(const Solution<Scalar>* sln);

      void set_dirichlet_lift(const Space<Scalar>* space, PrecalcShapeset* pss = NULL);

      virtual void save(const char* filename) const;

      void load(const char* filename, Space<Scalar>* space);

      Scalar get_ref_value(Element* e, double xi1, double xi2, int component = 0, int item = 0);

      Scalar get_ref_value_transformed(Element* e, double xi1, double xi2, int a, int b);

      virtual Func<Scalar>* get_pt_value(double x, double y);

      void multiply(Scalar coef);

      inline SolutionType get_type() const { return sln_type; };

      inline SpaceType get_space_type() const { return space_type; };

      static void vector_to_solutions(const Scalar* solution_vector, Hermes::vector<const Space<Scalar> *> spaces,
          Hermes::vector<Solution<Scalar>*> solutions,
          Hermes::vector<bool> add_dir_lift = Hermes::vector<bool>(),
          Hermes::vector<int> start_indices = Hermes::vector<int>());

      static void vector_to_solution(const Scalar* solution_vector, const Space<Scalar>* space, Solution<Scalar>* solution,
          bool add_dir_lift = true, int start_index = 0);

      static void vector_to_solutions(const Vector<Scalar>* vec, Hermes::vector<const Space<Scalar> *> spaces,
          Hermes::vector<Solution<Scalar>*> solutions,
          Hermes::vector<bool> add_dir_lift = Hermes::vector<bool>(),
          Hermes::vector<int> start_indices = Hermes::vector<int>());

      static void vector_to_solutions_common_dir_lift(const Vector<Scalar>* vec, Hermes::vector<const Space<Scalar> *> spaces,
          Hermes::vector<Solution<Scalar>*> solutions,
          bool add_dir_lift = false);

      static void vector_to_solutions_common_dir_lift(const Scalar* solution_vector, Hermes::vector<const Space<Scalar> *> spaces,
          Hermes::vector<Solution<Scalar>*> solutions,
          bool add_dir_lift = false);

      static void vector_to_solution(const Vector<Scalar>* vec, const Space<Scalar>* space, Solution<Scalar>* solution,
          bool add_dir_lift = true, int start_index = 0);

      static void vector_to_solutions(const Scalar* solution_vector, Hermes::vector<const Space<Scalar> *> spaces,
          Hermes::vector<Solution<Scalar>*> solutions, Hermes::vector<PrecalcShapeset *> pss,
          Hermes::vector<bool> add_dir_lift = Hermes::vector<bool>(),
          Hermes::vector<int> start_indices = Hermes::vector<int>());

      static void vector_to_solution(const Scalar* solution_vector, const Space<Scalar>* space, Solution<Scalar>* solution,
          PrecalcShapeset* pss, bool add_dir_lift = true, int start_index = 0);

      virtual void set_active_element(Element* e);

      virtual MeshFunction<Scalar>* clone() const;

      static void set_static_verbose_output(bool verbose);

    protected:
      static bool static_verbose_output;

      virtual int get_edge_fn_order(int edge) { return MeshFunction<Scalar>::get_edge_fn_order(edge); }

      void enable_transform(bool enable = true);

      virtual void init();

      virtual void free();

      virtual void set_coeff_vector(const Space<Scalar>* space, const Vector<Scalar>* vec, bool add_dir_lift, int start_index);

      virtual void set_coeff_vector(const Space<Scalar>* space, PrecalcShapeset* pss, const Scalar* coeffs, bool add_dir_lift, int start_index);

      virtual void set_coeff_vector(const Space<Scalar>* space, const Scalar* coeffs, bool add_dir_lift, int start_index);

      SolutionType sln_type;
      SpaceType space_type;

      bool transform;

      std::map<uint64_t, LightArray<struct Function<Scalar>::Node*>*>* tables[H2D_MAX_QUADRATURES][H2D_SOLUTION_ELEMENT_CACHE_SIZE];

      Element* elems[H2D_MAX_QUADRATURES][H2D_SOLUTION_ELEMENT_CACHE_SIZE];
      int cur_elem, oldest[H2D_SOLUTION_ELEMENT_CACHE_SIZE];

      Scalar* mono_coeffs;  
      int* elem_coeffs[H2D_MAX_SOLUTION_COMPONENTS];  

      int* elem_orders;
      int num_coeffs, num_elems;
      int num_dofs;

      void transform_values(int order, struct Function<Scalar>::Node* node, int newmask, int oldmask, int np);

      virtual void precalculate(int order, int mask);

      Scalar* dxdy_coeffs[H2D_MAX_SOLUTION_COMPONENTS][6];

      Scalar* dxdy_buffer;

      double** calc_mono_matrix(int o, int*& perm);

      void init_dxdy_buffer();

      void free_tables();

      Element* e_last; 

      friend class RefMap;
      template<typename T> friend class KellyTypeAdapt;
      template<typename T> friend class CalculationContinuity;
      template<typename T> friend class OGProjection;
      template<typename T> friend class OGProjectionNOX;
      template<typename T> friend class Adapt;
      template<typename T> friend class Func;
      template<typename T> friend class DiscontinuousFunc;
      template<typename T> friend class DiscreteProblem;
      template<typename T> friend class DiscreteProblemLinear;
      template<typename T> friend class NeighborSearch;
      template<typename T> friend HERMES_API Func<T>* init_fn(Solution<T>*fu, const int order);
      template<typename T> friend class RefinementSelectors::ProjBasedSelector;
      template<typename T> friend class RefinementSelectors::H1ProjBasedSelector;
      template<typename T> friend class RefinementSelectors::L2ProjBasedSelector;
      template<typename T> friend class RefinementSelectors::HcurlProjBasedSelector;
    };
  }
}

#endif