error_calculator.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_ERROR_CALCULATOR_H
#define __H2D_ERROR_CALCULATOR_H

#include "../weakform/weakform.h"
#include "norm_form.h"
#include "error_thread_calculator.h"

namespace Hermes
{
  namespace Hermes2D
  {
    template<typename Scalar> class ErrorThreadCalculator;

    enum CalculatedErrorType
    {
      AbsoluteError,
      RelativeErrorToElementNorm,
      RelativeErrorToGlobalNorm
    };


    template<typename Scalar>
    class HERMES_API ErrorCalculator :
      public Hermes::Mixins::TimeMeasurable,
      public Hermes::Mixins::Loggable,
      public Hermes::Hermes2D::Mixins::Parallel,
      public Hermes::Mixins::StateQueryable
    {
    public:
      ErrorCalculator(CalculatedErrorType errorType);

      void calculate_errors(std::vector<MeshFunctionSharedPtr<Scalar> > coarse_solutions, std::vector<MeshFunctionSharedPtr<Scalar> > fine_solutions, bool sort_and_store = true);

      void calculate_errors(MeshFunctionSharedPtr<Scalar> coarse_solution, MeshFunctionSharedPtr<Scalar> fine_solution, bool sort_and_store = true);

      virtual ~ErrorCalculator();

      void add_error_form(NormFormVol<Scalar>* form);
      void add_error_form(NormFormSurf<Scalar>* form);
      void add_error_form(NormFormDG<Scalar>* form);


      double get_element_error_squared(int component, int element_id) const;
      double get_element_norm_squared(int component, int element_id) const;
      double get_error_squared(int component) const;
      double get_norm_squared(int component) const;
      double get_total_error_squared() const;
      double get_total_norm_squared() const;

      int get_component_count() const { return this->component_count; }

      struct ElementReference {
        ElementReference() : element_id(-1), comp(-1), error(nullptr), norm(nullptr) {};
        ElementReference(int comp, int element_id, double* error, double* norm) : element_id(element_id), comp(comp), error(error), norm(norm) {};

        int element_id;
        int comp;
        double* error;
        double* norm;
      };

      const ElementReference& get_element_reference(unsigned int id) const;

      MeshFunctionSharedPtr<double> get_errorMeshFunction(int component = 0);

      void free();

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

      bool data_prepared_for_querying() const;

      void init_data_storage();

      void postprocess_error();

      std::vector<MeshFunctionSharedPtr<Scalar> > coarse_solutions;
      std::vector<MeshFunctionSharedPtr<Scalar> > fine_solutions;

      CalculatedErrorType errorType;

      ElementReference* element_references;

      int component_count;

      int num_act_elems;

      double* errors[H2D_MAX_COMPONENTS];
      double* norms[H2D_MAX_COMPONENTS];
      double component_errors[H2D_MAX_COMPONENTS];
      double component_norms[H2D_MAX_COMPONENTS];
      int element_count[H2D_MAX_COMPONENTS];
      double  errors_squared_sum;
      double  norms_squared_sum;

      MeshFunctionSharedPtr<double> errorMeshFunction[H2D_MAX_COMPONENTS];

      std::vector<NormFormVol<Scalar> *> mfvol;
      std::vector<NormFormSurf<Scalar> *> mfsurf;
      std::vector<NormFormDG<Scalar> *> mfDG;

      bool elements_stored;

      static int compareElementReference(const void * a, const void * b)
      {
        ElementReference* ref_a = (ElementReference*)(a);
        ElementReference* ref_b = (ElementReference*)(b);
        if (*((*ref_a).error) > *((*ref_b).error))
          return -1;
        else
          return 1;
      };

      friend class Adapt < Scalar > ;
      friend class ErrorThreadCalculator < Scalar > ;
    };

    template<typename Scalar, NormType normType>
    class HERMES_API DefaultErrorCalculator : public ErrorCalculator < Scalar >
    {
    public:
      DefaultErrorCalculator(CalculatedErrorType errorType, int component_count);
      virtual ~DefaultErrorCalculator();
    };

    template<typename Scalar, NormType normType>
    class HERMES_API DefaultNormCalculator : public ErrorCalculator < Scalar >
    {
    public:
      DefaultNormCalculator(int component_count);
      virtual ~DefaultNormCalculator();

      double calculate_norms(std::vector<MeshFunctionSharedPtr<Scalar> >& solutions);
      double calculate_norm(MeshFunctionSharedPtr<Scalar> solution);
    };
  }
}
#endif