Hermes2D  3.0
kelly_type_adapt.h
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 /*
16 #ifndef KELLY_TYPE_ADAPT_H
17 #define KELLY_TYPE_ADAPT_H
18 
19 #include "adapt.h"
20 #include "../neighbor_search.h"
21 #include "../discrete_problem.h"
22 
23 namespace Hermes
24 {
25 namespace Hermes2D
26 {
28 class InterfaceEstimatorScalingFunction
29 {
30 private:
31 virtual double value(double e_diam, const std::string& e_marker) const = 0;
32 template<typename Scalar> friend class KellyTypeAdapt;
33 };
34 
36 class ScaleByElementDiameter : public InterfaceEstimatorScalingFunction
37 {
38 private:
39 virtual double value(double e_diam, const std::string& e_marker) const {
40 return e_diam;
41 }
42 template<typename Scalar> friend class KellyTypeAdapt;
43 };
44 
65 template<typename Scalar>
66 class HERMES_API KellyTypeAdapt : public Adapt<Scalar>
67 {
68 public:
90 class HERMES_API ErrorEstimatorForm : public Form<Scalar>
91 {
92 public:
94 int i;
96 std::string area;
98 std::vector<MeshFunctionSharedPtr<Scalar> > ext;
100 void setAsInterface();
102 ErrorEstimatorForm(int i, std::string area = HERMES_ANY,
103 std::vector<MeshFunctionSharedPtr<Scalar> > ext = std::vector<MeshFunctionSharedPtr<Scalar> >())
104 : i(i), area(area), ext(ext) {}
105 
107 virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[],
108 DiscontinuousFunc<Scalar> *u, Geom<double> *e,
109 Func<Scalar> **ext) const
110 {
111 throw Exceptions::MethodNotOverridenException("KellyTypeAdapt::ErrorEstimatorForm::value()");
112 return 0.0;
113 }
114 
116 virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[],
117 DiscontinuousFunc<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
118 Func<Ord> **ext) const
119 {
120 throw Exceptions::MethodNotOverridenException("KellyTypeAdapt::ErrorEstimatorForm::ord().");
121 return Hermes::Ord();
122 }
123 
125 KellyTypeAdapt *adapt;
126 };
127 
128 protected:
129 // Only needed for gaining access to NeighborSearch methods.
130 DiscreteProblem<Scalar> dp;
131 
135 double eval_volumetric_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
136 RefMap* rm);
137 double eval_boundary_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
138 RefMap* rm,
139 SurfPos* surf_pos);
140 double eval_interface_estimator(typename KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
141 RefMap *rm,
142 SurfPos* surf_pos,
143 LightArray<NeighborSearch<Scalar>*>& neighbor_searches,
144 int neighbor_index);
145 double eval_solution_norm(typename Adapt<Scalar>::MatrixFormVolError* form,
146 RefMap* rm,
147 MeshFunctionSharedPtr<Scalar> sln);
148 
152 std::vector<typename KellyTypeAdapt::ErrorEstimatorForm *> error_estimators_vol;
153 std::vector<typename KellyTypeAdapt::ErrorEstimatorForm *> error_estimators_surf;
154 
155 Mesh::ElementMarkersConversion element_markers_conversion;
156 Mesh::BoundaryMarkersConversion boundary_markers_conversion;
157 
160 std::vector<const InterfaceEstimatorScalingFunction*> interface_scaling_fns;
162 bool use_aposteriori_interface_scaling;
165 
170 double interface_scaling_const;
172 double volumetric_scaling_const;
174 double boundary_scaling_const;
175 
179 bool ignore_visited_segments;
180 
185 virtual double calc_err_internal(std::vector<MeshFunctionSharedPtr<Scalar> > slns,
186 std::vector<double>* component_errors,
187 unsigned int error_flags);
188 
189 public:
190 
220 KellyTypeAdapt(std::vector<SpaceSharedPtr<Scalar> > spaces,
221 bool ignore_visited_segments = true,
222 std::vector<const InterfaceEstimatorScalingFunction*>
223 interface_scaling_fns_ = std::vector<const InterfaceEstimatorScalingFunction*>(),
224 std::vector<NormType> norms_ = std::vector<NormType>());
225 
226 KellyTypeAdapt(SpaceSharedPtr<Scalar> space,
227 bool ignore_visited_segments = true,
228 const InterfaceEstimatorScalingFunction* interface_scaling_fn_ = nullptr,
229 NormType norm_ = HERMES_UNSET_NORM);
230 
232 virtual ~KellyTypeAdapt()
233 {
234 for (unsigned int i = 0; i < error_estimators_surf.size(); i++)
235 delete error_estimators_surf[i];
236 error_estimators_surf.clear();
237 
238 for (unsigned int i = 0; i < error_estimators_vol.size(); i++)
239 delete error_estimators_vol[i];
240 error_estimators_vol.clear();
241 
242 for (unsigned int i = 0; i < interface_scaling_fns.size(); i++)
243 delete interface_scaling_fns[i];
244 interface_scaling_fns.clear();
245 }
246 
247 Mesh::ElementMarkersConversion* get_element_markers_conversion()
248 {
249 return &this->element_markers_conversion;
250 }
251 Mesh::BoundaryMarkersConversion* get_boundary_markers_conversion()
252 {
253 return &this->boundary_markers_conversion;
254 }
255 
263 void add_error_estimator_vol(ErrorEstimatorForm* form);
264 
271 void add_error_estimator_surf(ErrorEstimatorForm* form);
272 
276 
277 double calc_err_est(MeshFunctionSharedPtr<Scalar>sln,
278 unsigned int error_flags = HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL)
279 {
280 if(this->num != 1)
281 throw Exceptions::Exception("Wrong number of solutions.");
282 std::vector<MeshFunctionSharedPtr<Scalar> > slns;
283 slns.push_back(sln);
284 return calc_err_est(slns, nullptr, error_flags);
285 }
286 
287 double calc_err_est(std::vector<MeshFunctionSharedPtr<Scalar> > slns,
288 std::vector<double>* component_errors = nullptr,
289 unsigned int error_flags = HERMES_TOTAL_ERROR_REL | HERMES_ELEMENT_ERROR_REL)
290 {
291 return calc_err_internal(slns, component_errors, error_flags);
292 }
293 
297 bool adapt(double thr, int strat = 0, int regularize = -1, double to_be_processed = 0.0);
298 
299 void disable_aposteriori_interface_scaling() { use_aposteriori_interface_scaling = false; }
300 
301 void set_volumetric_scaling_const(double C) { volumetric_scaling_const = C; }
302 void set_boundary_scaling_const(double C) { boundary_scaling_const = C; }
303 };
304 
319 template<typename Scalar>
320 class HERMES_API BasicKellyAdapt : public KellyTypeAdapt<Scalar>
321 {
322 public:
323 class HERMES_API ErrorEstimatorFormKelly : public KellyTypeAdapt<Scalar>::ErrorEstimatorForm
324 {
325 public:
327 ErrorEstimatorFormKelly(int i = 0, double const_by_laplacian = 1.0);
328 
329 virtual Scalar value(int n, double *wt, Func<Scalar> *u_ext[],
330 DiscontinuousFunc<Scalar> *u, Geom<double> *e,
331 Func<Scalar> **ext) const
332 {
333 Scalar result = 0.;
334 if(u->fn_central != nullptr)
335 for (int i = 0; i < n; i++)
336 result += wt[i] * Hermes::sqr( const_by_laplacian * ( e->nx[i] * u->dx[i] + e->ny[i] * u->dy[i]));
337 else
338 for (int i = 0; i < n; i++)
339 result += wt[i] * Hermes::sqr( const_by_laplacian * ( e->nx[i] * u->dx_neighbor[i] + e->ny[i] * u->dy_neighbor[i]));
340 
341 return result;
342 }
343 
344 virtual Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[],
345 DiscontinuousFunc<Hermes::Ord> *u, Geom<Hermes::Ord> *e,
346 Func<Ord> **ext) const
347 {
348 if(u->fn_central != nullptr)
349 return Hermes::sqr(u->dx[0] + u->dy[0]);
350 else
351 return Hermes::sqr(u->dx_neighbor[0] + u->dy_neighbor[0]);
352 }
353 
354 private:
355 double const_by_laplacian;
356 };
357 
362 BasicKellyAdapt(std::vector<SpaceSharedPtr<Scalar> > spaces_,
363 double const_by_laplacian = 1.0,
364 std::vector<NormType> norms_ = std::vector<NormType>())
365 : KellyTypeAdapt<Scalar>(spaces_, true, std::vector<const InterfaceEstimatorScalingFunction*>(), norms_)
366 {
367 set_scaling_consts(const_by_laplacian);
368 for (int i = 0; i < this->num; i++)
369 this->error_estimators_surf.push_back(new ErrorEstimatorFormKelly(i, const_by_laplacian));
370 }
371 
372 BasicKellyAdapt(SpaceSharedPtr<Scalar> space_, double const_by_laplacian = 1.0, NormType norm_ = HERMES_UNSET_NORM)
373 : KellyTypeAdapt<Scalar>(space_, true, nullptr, norm_)
374 {
375 set_scaling_consts(const_by_laplacian);
376 this->error_estimators_surf.push_back(new ErrorEstimatorFormKelly(0, const_by_laplacian));
377 }
378 
379 private:
380 void set_scaling_consts(double C)
381 {
382 this->interface_scaling_const = 1./(24.*C);
383 this->volumetric_scaling_const = this->interface_scaling_const;
384 this->boundary_scaling_const = this->interface_scaling_const;
385 }
386 };
387 }
388 }
389 #endif
390 */