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_search.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;
std::vector<MeshFunctionSharedPtr<Scalar> > ext;
void setAsInterface();
ErrorEstimatorForm(int i, std::string area = HERMES_ANY,
std::vector<MeshFunctionSharedPtr<Scalar> > ext = std::vector<MeshFunctionSharedPtr<Scalar> >())
: i(i), area(area), ext(ext) {}

virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[],
DiscontinuousFunc<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[],
DiscontinuousFunc<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
Func<Ord> **ext) const
{
throw Exceptions::MethodNotOverridenException("KellyTypeAdapt::ErrorEstimatorForm::ord().");
return Hermes::Ord();
}

KellyTypeAdapt *adapt;
};

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

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,
MeshFunctionSharedPtr<Scalar> sln);

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

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

std::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(std::vector<MeshFunctionSharedPtr<Scalar> > slns,
std::vector<double>* component_errors,
unsigned int error_flags);

public:

KellyTypeAdapt(std::vector<SpaceSharedPtr<Scalar> > spaces,
bool ignore_visited_segments = true,
std::vector<const InterfaceEstimatorScalingFunction*>
interface_scaling_fns_ = std::vector<const InterfaceEstimatorScalingFunction*>(),
std::vector<NormType> norms_ = std::vector<NormType>());

KellyTypeAdapt(SpaceSharedPtr<Scalar> space,
bool ignore_visited_segments = true,
const InterfaceEstimatorScalingFunction* interface_scaling_fn_ = nullptr,
NormType 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(MeshFunctionSharedPtr<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.");
std::vector<MeshFunctionSharedPtr<Scalar> > slns;
slns.push_back(sln);
return calc_err_est(slns, nullptr, error_flags);
}

double calc_err_est(std::vector<MeshFunctionSharedPtr<Scalar> > slns,
std::vector<double>* component_errors = nullptr,
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[],
DiscontinuousFunc<Scalar> *u, Geom<double> *e,
Func<Scalar> **ext) const
{
Scalar result = 0.;
if(u->fn_central != nullptr)
for (int i = 0; i < n; i++)
result += wt[i] * Hermes::sqr( const_by_laplacian * ( e->nx[i] * u->dx[i] + e->ny[i] * u->dy[i]));
else
for (int i = 0; i < n; i++)
result += wt[i] * Hermes::sqr( const_by_laplacian * ( e->nx[i] * u->dx_neighbor[i] + e->ny[i] * u->dy_neighbor[i]));

return result;
}

virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[],
DiscontinuousFunc<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
Func<Ord> **ext) const
{
if(u->fn_central != nullptr)
return Hermes::sqr(u->dx[0] + u->dy[0]);
else
return Hermes::sqr(u->dx_neighbor[0] + u->dy_neighbor[0]);
}

private:
double const_by_laplacian;
};

BasicKellyAdapt(std::vector<SpaceSharedPtr<Scalar> > spaces_,
double const_by_laplacian = 1.0,
std::vector<NormType> norms_ = std::vector<NormType>())
: KellyTypeAdapt<Scalar>(spaces_, true, std::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(SpaceSharedPtr<Scalar> space_, double const_by_laplacian = 1.0, NormType norm_ = HERMES_UNSET_NORM)
: KellyTypeAdapt<Scalar>(space_, true, nullptr, 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
*/