kelly_type_adapt.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 KELLY_TYPE_ADAPT_H
#define KELLY_TYPE_ADAPT_H

#include "adapt.h"
#include "../neighbor.h"
#include "../discrete_problem.h"

namespace Hermes
{
  namespace Hermes2D
  {
    class InterfaceEstimatorScalingFunction
    {
      private:
        virtual double value(double e_diam, const std::string& e_marker) const = 0;
        template<typename Scalar> friend class KellyTypeAdapt;
    };

    class ScaleByElementDiameter : public InterfaceEstimatorScalingFunction
    {
      private:
        virtual double value(double e_diam, const std::string& e_marker) const {
          return e_diam;
        }
        template<typename Scalar> friend class KellyTypeAdapt;
    };

    template<typename Scalar>
    class HERMES_API KellyTypeAdapt : public Adapt<Scalar>
    {
    public:
      class HERMES_API ErrorEstimatorForm : public Form<Scalar>
      {
      public:
        int i; 
        std::string area; 
        Hermes::vector<MeshFunction<Scalar>*> ext; 

        void setAsInterface();
        ErrorEstimatorForm(int i, std::string area = HERMES_ANY,
                           Hermes::vector<MeshFunction<Scalar>*> ext = Hermes::vector<MeshFunction<Scalar>*>())
          : i(i), area(area), ext(ext) {}

        virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[],
                             Func<Scalar> *u, Geom<double> *e,
                             Func<Scalar> **ext) const
        {
          throw Exceptions::MethodNotOverridenException("KellyTypeAdapt::ErrorEstimatorForm::value()");
          return 0.0;
        }

        virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[],
                                Func<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
                                Func<Ord> **ext) const
        {
          throw Exceptions::MethodNotOverridenException("KellyTypeAdapt::ErrorEstimatorForm::ord().");
          return Hermes::Ord();
        }

        KellyTypeAdapt *adapt;
      };

    protected:
      DiscreteProblem<Scalar> dp; // Only needed for gaining access to NeighborSearch methods.

      double eval_volumetric_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
                                       RefMap* rm);
      double eval_boundary_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
                                     RefMap* rm,
                                     SurfPos* surf_pos);
      double eval_interface_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
                                      RefMap *rm,
                                      SurfPos* surf_pos,
                                      LightArray<NeighborSearch<Scalar>*>& neighbor_searches,
                                      int neighbor_index);
      double eval_solution_norm(typename Adapt<Scalar>::MatrixFormVolError* form,
                                RefMap* rm,
                                MeshFunction<Scalar>* sln);

      Hermes::vector<typename KellyTypeAdapt::ErrorEstimatorForm *> error_estimators_vol;
      Hermes::vector<typename KellyTypeAdapt::ErrorEstimatorForm *> error_estimators_surf;

      Mesh::ElementMarkersConversion element_markers_conversion;
      Mesh::BoundaryMarkersConversion boundary_markers_conversion;

      Hermes::vector<const InterfaceEstimatorScalingFunction*> interface_scaling_fns;
      bool use_aposteriori_interface_scaling; 



      double interface_scaling_const;   
      double volumetric_scaling_const;  
      double boundary_scaling_const;    

      bool ignore_visited_segments;

      virtual double calc_err_internal(Hermes::vector<Solution<Scalar>*> slns,
                                       Hermes::vector<double>* component_errors,
                                       unsigned int error_flags);

    public:

      KellyTypeAdapt(Hermes::vector<Space<Scalar>*> spaces,
                     bool ignore_visited_segments = true,
                     Hermes::vector<const InterfaceEstimatorScalingFunction*>
                       interface_scaling_fns_ = Hermes::vector<const InterfaceEstimatorScalingFunction*>(),
                     Hermes::vector<ProjNormType> norms_ = Hermes::vector<ProjNormType>());

      KellyTypeAdapt(Space<Scalar>* space,
                     bool ignore_visited_segments = true,
                     const InterfaceEstimatorScalingFunction* interface_scaling_fn_ = NULL,
                     ProjNormType norm_ = HERMES_UNSET_NORM);

      virtual ~KellyTypeAdapt()
      {
        for (unsigned int i = 0; i < error_estimators_surf.size(); i++)
          delete error_estimators_surf[i];
        error_estimators_surf.clear();

        for (unsigned int i = 0; i < error_estimators_vol.size(); i++)
          delete error_estimators_vol[i];
        error_estimators_vol.clear();

        for (unsigned int i = 0; i < interface_scaling_fns.size(); i++)
          delete interface_scaling_fns[i];
        interface_scaling_fns.clear();
      }

      Mesh::ElementMarkersConversion* get_element_markers_conversion()
      {
        return &this->element_markers_conversion;
      }
      Mesh::BoundaryMarkersConversion* get_boundary_markers_conversion()
      {
        return &this->boundary_markers_conversion;
      }

      void add_error_estimator_vol(ErrorEstimatorForm* form);

      void add_error_estimator_surf(ErrorEstimatorForm* form);


      double calc_err_est(Solution<Scalar>*sln,
                          unsigned int error_flags = HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL)
      {
        if(this->num != 1)
          throw Exceptions::Exception("Wrong number of solutions.");
        Hermes::vector<Solution<Scalar>*> slns;
        slns.push_back(sln);
        return calc_err_est(slns, NULL, error_flags);
      }

      double calc_err_est(Hermes::vector<Solution<Scalar>*> slns,
                          Hermes::vector<double>* component_errors = NULL,
                          unsigned int error_flags = HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL)
      {
        return calc_err_internal(slns, component_errors, error_flags);
      }

      bool adapt(double thr, int strat = 0, int regularize = -1, double to_be_processed = 0.0);

      void disable_aposteriori_interface_scaling() { use_aposteriori_interface_scaling = false; }

      void set_volumetric_scaling_const(double C) { volumetric_scaling_const = C; }
      void set_boundary_scaling_const(double C) { boundary_scaling_const = C; }
    };

    template<typename Scalar>
    class HERMES_API BasicKellyAdapt : public KellyTypeAdapt<Scalar>
    {
    public:
      class HERMES_API ErrorEstimatorFormKelly : public KellyTypeAdapt<Scalar>::ErrorEstimatorForm
      {
      public:
        ErrorEstimatorFormKelly(int i = 0, double const_by_laplacian = 1.0);

        virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[],
                             Func<Scalar> *u, Geom<double> *e,
                             Func<Scalar> **ext) const
        {
          Scalar result = 0.;
          for (int i = 0; i < n; i++)
            result += wt[i] * Hermes::sqr( const_by_laplacian * ( e->nx[i] * (u->get_dx_central(i) - u->get_dx_neighbor(i)) +
                                                                  e->ny[i] * (u->get_dy_central(i) - u->get_dy_neighbor(i)) ) );
          return result;
        }

        virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[],
                                Func<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
                                Func<Ord> **ext) const
        {
          return Hermes::sqr( (u->get_dx_central(0) - u->get_dx_neighbor(0)) +
                              (u->get_dy_central(0) - u->get_dy_neighbor(0)) );
        }

      private:
        double const_by_laplacian;
      };

      BasicKellyAdapt(Hermes::vector<Space<Scalar>*> spaces_,
                      double const_by_laplacian = 1.0,
                      Hermes::vector<ProjNormType> norms_ = Hermes::vector<ProjNormType>())
        : KellyTypeAdapt<Scalar>(spaces_, true, Hermes::vector<const InterfaceEstimatorScalingFunction*>(), norms_)
      {
        set_scaling_consts(const_by_laplacian);
        for (int i = 0; i < this->num; i++)
          this->error_estimators_surf.push_back(new ErrorEstimatorFormKelly(i, const_by_laplacian));
      }

      BasicKellyAdapt(Space<Scalar>* space_, double const_by_laplacian = 1.0, ProjNormType norm_ = HERMES_UNSET_NORM)
        : KellyTypeAdapt<Scalar>(space_, true, NULL, norm_)
      {
        set_scaling_consts(const_by_laplacian);
        this->error_estimators_surf.push_back(new ErrorEstimatorFormKelly(0, const_by_laplacian));
      }

    private:
      void set_scaling_consts(double C)
      {
        this->interface_scaling_const = 1./(24.*C);
        this->volumetric_scaling_const = this->interface_scaling_const;
        this->boundary_scaling_const = this->interface_scaling_const;
      }
    };
  }
}
#endif