Line data Source code
1 : /*
2 : * Copyright (c) 2014-2017: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Naegel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__FUNCTION_SPACE__METRIC_SPACES_H_
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__METRIC_SPACES_H_
35 :
36 : // C++ includes
37 : #include <vector>
38 : #include <cmath>
39 :
40 : // UG includes
41 : #include "common/common.h"
42 : #include "common/util/smart_pointer.h"
43 :
44 : #include "lib_algebra/lib_algebra.h"
45 : #include "lib_algebra/operator/debug_writer.h"
46 :
47 : #include "lib_disc/function_spaces/integrate.h"
48 : #include "lib_disc/function_spaces/grid_function.h"
49 : #include "lib_disc/function_spaces/grid_function_user_data.h"
50 : #include "lib_disc/function_spaces/integrate.h"
51 :
52 : #include "lib_disc/spatial_disc/user_data/linker/darcy_velocity_linker.h" // required for Darcy velocity linker
53 :
54 : namespace ug {
55 : /// Abstract base class for (algebraic) vectors
56 : template<typename TVector>
57 : class IBanachSpace
58 : {
59 :
60 : public:
61 : virtual ~IBanachSpace() {}
62 :
63 : /// euclidean norm (default)
64 : virtual double norm(TVector &x)
65 : { return x.norm(); }
66 :
67 : virtual double distance(TVector& x, TVector& y)
68 : {
69 : SmartPtr<TVector> delta = x.clone();
70 : *delta -= y;
71 : return norm(*delta);
72 : }
73 :
74 : };
75 :
76 :
77 :
78 : /// Abstract base class for grid functions
79 : template<typename TGridFunction>
80 : class IGridFunctionSpace :
81 : public IBanachSpace<typename TGridFunction::vector_type>
82 : {
83 :
84 : public:
85 : typedef typename TGridFunction::vector_type vector_type;
86 : typedef TGridFunction grid_function_type;
87 :
88 : /// DTOR
89 : virtual ~IGridFunctionSpace() {}
90 :
91 : /// norm (for grid functions)
92 : virtual double norm(TGridFunction& x) = 0;
93 : virtual double norm2(TGridFunction& x) = 0;
94 :
95 : // { return x.norm(); }
96 :
97 : /// distance (for grid functions)
98 : virtual double distance(TGridFunction& x, TGridFunction& y) = 0;
99 : virtual double distance2(TGridFunction& x, TGridFunction& y) = 0;
100 : /*{
101 : SmartPtr<TGridFunction> delta = x.clone();
102 : *delta -= y;
103 : return norm(*delta);
104 : }*/
105 :
106 : /// OVERRIDE norm (for vectors)
107 0 : virtual double norm(vector_type &x)
108 : {
109 0 : TGridFunction* gfX=dynamic_cast< TGridFunction*>(&x);
110 : UG_ASSERT(gfX!=NULL, "Huhh: GridFunction required!");
111 0 : return norm(*gfX);
112 : }
113 :
114 : /// OVERRIDE distance (for vectors)
115 0 : virtual double distance(vector_type &x, vector_type &y)
116 0 : { return distance(static_cast<TGridFunction &>(x), static_cast<TGridFunction &>(y)); }
117 :
118 0 : virtual double scaling() const
119 0 : { return 1.0; }
120 :
121 : virtual std::string config_string() const
122 : { return std::string("IGridFunctionSpace"); }
123 : };
124 :
125 : template<typename TGridFunction>
126 : class AlgebraicSpace : public IGridFunctionSpace<TGridFunction>
127 : {
128 : typedef typename TGridFunction::vector_type vector_type;
129 : virtual ~AlgebraicSpace() {}
130 :
131 : virtual double norm(TGridFunction& x)
132 : {return x.norm();}
133 :
134 : virtual double norm2(TGridFunction& x)
135 : { double n = this->norm(x); return n*n; }
136 :
137 : virtual double distance(TGridFunction& x, TGridFunction& y)
138 : { SmartPtr<TGridFunction> delta = x.clone(); *delta -= y; return delta->norm(); }
139 :
140 : virtual double distance2(TGridFunction& x, TGridFunction& y)
141 : { double d = this->distance(x,y);; return d*d;}
142 : };
143 :
144 :
145 :
146 : /** Auxiliary class for providing weights*/
147 : template <typename W>
148 0 : class IObjectWithWeights
149 : {
150 : public:
151 : typedef W weight_type;
152 :
153 : IObjectWithWeights()
154 : : m_spWeight(SPNULL) {}
155 :
156 : IObjectWithWeights(ConstSmartPtr<weight_type> spW)
157 : : m_spWeight(spW) {}
158 :
159 : /// for weighted norms
160 0 : void set_weight(ConstSmartPtr<weight_type> spWeight)
161 0 : { m_spWeight = spWeight; }
162 :
163 0 : ConstSmartPtr<weight_type> get_weight()
164 0 : { return m_spWeight; }
165 :
166 : protected:
167 : ConstSmartPtr<weight_type> m_spWeight;
168 : };
169 :
170 :
171 : /** Auxiliary class for time dependence - SHOULD be replaced by product space*/
172 : /*
173 : class ITimeData {
174 :
175 : public:
176 : virtual ~ITimeData() {};
177 :
178 : /// characteristic time
179 : virtual void update_time_data(number dt) = 0;
180 :
181 : virtual bool is_time_dependent() const
182 : { return false; }
183 :
184 : };
185 :
186 : */
187 :
188 : //! Estimate the error (based on the difference between two grid functions)
189 : template <typename TGridFunction>
190 : class IComponentSpace :
191 : public IGridFunctionSpace<TGridFunction>
192 : {
193 : protected:
194 : std::string m_fctNames;
195 : const char* m_ssNames;
196 : int m_quadorder;
197 : //bool m_opposite;
198 : //number m_thickness;
199 : //std::string m_type;
200 :
201 : public:
202 : typedef IGridFunctionSpace<TGridFunction> base_type;
203 : static const int dim=TGridFunction::dim;
204 :
205 0 : IComponentSpace(const char *fctNames)
206 0 : : m_fctNames(fctNames), m_ssNames(NULL), m_quadorder(3){}
207 :
208 0 : IComponentSpace(const char *fctNames, int order)
209 0 : : m_fctNames(fctNames), m_ssNames(NULL), m_quadorder(order) {}
210 :
211 0 : IComponentSpace(const char *fctNames, const char* ssNames, int order)
212 0 : : m_fctNames(fctNames), m_ssNames(ssNames), m_quadorder(order){}
213 :
214 0 : virtual ~IComponentSpace() {};
215 :
216 : // per convention, norm must return sqrt of norm2
217 0 : virtual double norm(TGridFunction& uFine)
218 0 : { return sqrt(norm2(uFine)); }
219 :
220 : // per convention, distance must return sqrt of distance2
221 0 : virtual double distance(TGridFunction& uFine, TGridFunction& uCoarse)
222 0 : { return sqrt(distance2(uFine, uCoarse)); }
223 :
224 : virtual double norm2(TGridFunction& uFine) = 0;
225 : virtual double distance2(TGridFunction& uFine, TGridFunction& uCoarse) = 0;
226 :
227 : std::string function_name(){return m_fctNames;}
228 :
229 : //void set_opposite(bool opposite){m_opposite=opposite;}
230 :
231 : public:
232 : /// print config string
233 0 : virtual std::string config_string() const
234 : {
235 0 : std::stringstream ss;
236 :
237 0 : if (this->m_ssNames)
238 0 : ss << this->m_fctNames << ", " << this->m_ssNames << ", " << this->m_quadorder
239 0 : << ", type=" << std::endl;
240 :
241 : else
242 0 : ss << this->m_fctNames << ", (no name), " << this->m_quadorder
243 0 : << ", type=" << std::endl;
244 0 : return ss.str();
245 0 : }
246 :
247 : };
248 :
249 :
250 : /// Computes simple vector (2-)norm from grid function DoFs which belong to either of the
251 : /// subsets and either of the functions supplied in the constructor.
252 : template <typename TGridFunction>
253 : class GridFunctionComponentSpace
254 : : public IComponentSpace<TGridFunction>
255 : {
256 : public:
257 : GridFunctionComponentSpace(const char* fctNames)
258 0 : : IComponentSpace<TGridFunction>(fctNames) {}
259 :
260 : GridFunctionComponentSpace(const char* fctNames, const char* ssNames)
261 0 : : IComponentSpace<TGridFunction>(fctNames, ssNames, 1) {}
262 :
263 0 : virtual ~GridFunctionComponentSpace() {};
264 :
265 0 : virtual double norm2(TGridFunction& uFine)
266 : {
267 0 : ConstSmartPtr<DoFDistribution> dd = uFine.dof_distribution();
268 :
269 : // find function indices
270 0 : FunctionGroup fg(dd->function_pattern());
271 0 : try {fg.add(TokenizeTrimString(this->m_fctNames));}
272 0 : UG_CATCH_THROW("Could not convert function names to function indices.");
273 :
274 : // find subset indices
275 0 : SubsetGroup sg(dd->subset_handler());
276 0 : if (this->m_ssNames)
277 : {
278 0 : try {sg.add(TokenizeTrimString(this->m_ssNames));}
279 0 : UG_CATCH_THROW("Could not convert subset names to subset indices.");
280 : }
281 : else
282 0 : sg.add_all();
283 :
284 :
285 : // loop subsets
286 0 : number sum = 0.0;
287 : const size_t sgSz = sg.size();
288 0 : for (size_t i = 0; i < sgSz; ++i)
289 : {
290 : int si = sg[i];
291 0 : if (dd->max_dofs(VERTEX, si)) add_norm_values<Vertex>(sum, uFine, dd, sg[i], fg);
292 0 : if (dd->max_dofs(EDGE, si)) add_norm_values<Edge>(sum, uFine, dd, sg[i], fg);
293 0 : if (dd->max_dofs(FACE, si)) add_norm_values<Face>(sum, uFine, dd, sg[i], fg);
294 0 : if (dd->max_dofs(VOLUME, si)) add_norm_values<Volume>(sum, uFine, dd, sg[i], fg);
295 : }
296 :
297 : #ifdef UG_PARALLEL
298 : if (pcl::NumProcs() > 1)
299 : {
300 : pcl::ProcessCommunicator pc;
301 : sum = pc.allreduce(sum, PCL_RO_SUM);
302 : }
303 : #endif
304 :
305 0 : return sum;
306 : }
307 :
308 0 : virtual double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
309 : {
310 0 : ConstSmartPtr<DoFDistribution> dd = uFine.dof_distribution();
311 0 : UG_COND_THROW(dd != uCoarse.dof_distribution(),
312 : "GridFunctionComponentSpace::distance2: GF1 DoF distro is not the same as for GF2.\n"
313 : "This case is not implemented.");
314 :
315 : // find function indices
316 0 : FunctionGroup fg(dd->function_pattern());
317 0 : try {fg.add(TokenizeTrimString(this->m_fctNames));}
318 0 : UG_CATCH_THROW("Could not convert function names to function indices.");
319 :
320 : // find subset indices
321 0 : SubsetGroup sg(dd->subset_handler());
322 0 : if (this->m_ssNames)
323 : {
324 0 : try {sg.add(TokenizeTrimString(this->m_ssNames));}
325 0 : UG_CATCH_THROW("Could not convert subset names to subset indices.");
326 : }
327 : else
328 0 : sg.add_all();
329 :
330 :
331 : // loop subsets
332 0 : number sum = 0.0;
333 : const size_t sgSz = sg.size();
334 0 : for (size_t i = 0; i < sgSz; ++i)
335 : {
336 : int si = sg[i];
337 0 : if (dd->max_dofs(VERTEX, si)) add_distance_values<Vertex>(sum, uFine, uCoarse, dd, sg[i], fg);
338 0 : if (dd->max_dofs(EDGE, si)) add_distance_values<Edge>(sum, uFine, uCoarse, dd, sg[i], fg);
339 0 : if (dd->max_dofs(FACE, si)) add_distance_values<Face>(sum, uFine, uCoarse, dd, sg[i], fg);
340 0 : if (dd->max_dofs(VOLUME, si)) add_distance_values<Volume>(sum, uFine, uCoarse, dd, sg[i], fg);
341 : }
342 :
343 : #ifdef UG_PARALLEL
344 : if (pcl::NumProcs() > 1)
345 : {
346 : pcl::ProcessCommunicator pc;
347 : sum = pc.allreduce(sum, PCL_RO_SUM);
348 : }
349 : #endif
350 :
351 0 : return sum;
352 : }
353 :
354 : private:
355 : template <typename TBaseElem>
356 0 : void add_distance_values
357 : (
358 : number& sum,
359 : const TGridFunction& uFine,
360 : const TGridFunction& uCoarse,
361 : ConstSmartPtr<DoFDistribution> dd,
362 : int si,
363 : const FunctionGroup& fg
364 : ) const
365 : {
366 0 : const SurfaceView& sv = *dd->surface_view();
367 0 : const MultiGrid& mg = *dd->multi_grid();
368 :
369 0 : const number nFct = fg.size();
370 :
371 : // iterate all elements (including SHADOW_RIM_COPY!)
372 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
373 0 : iter = dd->template begin<TBaseElem>(si, SurfaceView::ALL);
374 0 : iterEnd = dd->template end<TBaseElem>(si, SurfaceView::ALL);
375 : std::vector<DoFIndex> vInd;
376 0 : for (; iter != iterEnd; ++iter)
377 : {
378 : TBaseElem* elem = *iter;
379 0 : if (sv.is_contained(elem, dd->grid_level(), SurfaceView::SHADOW_RIM_COPY))
380 : {
381 0 : if (mg.num_children<TBaseElem>(elem) > 0)
382 : {
383 : TBaseElem* child = mg.get_child<TBaseElem>(elem, 0);
384 0 : if (sv.is_contained(child, dd->grid_level(), SurfaceView::SURFACE_RIM))
385 0 : continue;
386 : }
387 : }
388 :
389 0 : for (size_t fi = 0; fi < nFct; ++fi)
390 : {
391 0 : dd->inner_dof_indices(elem, fg[fi], vInd);
392 :
393 : const size_t nDof = vInd.size();
394 0 : for (size_t dof = 0; dof < nDof; ++dof)
395 : {
396 0 : const number tmp = DoFRef(uFine, vInd[dof]) - DoFRef(uCoarse, vInd[dof]);
397 0 : sum += tmp*tmp;
398 : }
399 : }
400 : }
401 : // note: no duplicate additions possible
402 0 : }
403 :
404 : template <typename TBaseElem>
405 0 : void add_norm_values
406 : (
407 : number& sum,
408 : const TGridFunction& uFine,
409 : ConstSmartPtr<DoFDistribution> dd,
410 : int si,
411 : const FunctionGroup& fg
412 : ) const
413 : {
414 0 : const SurfaceView& sv = *dd->surface_view();
415 0 : const MultiGrid& mg = *dd->multi_grid();
416 :
417 0 : const number nFct = fg.size();
418 :
419 : // iterate all elements (including SHADOW_RIM_COPY!)
420 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
421 0 : iter = dd->template begin<TBaseElem>(si, SurfaceView::ALL);
422 0 : iterEnd = dd->template end<TBaseElem>(si, SurfaceView::ALL);
423 : std::vector<DoFIndex> vInd;
424 0 : for (; iter != iterEnd; ++iter)
425 : {
426 : TBaseElem* elem = *iter;
427 0 : if (sv.is_contained(elem, dd->grid_level(), SurfaceView::SHADOW_RIM_COPY))
428 : {
429 0 : if (mg.num_children<TBaseElem>(elem) > 0)
430 : {
431 : TBaseElem* child = mg.get_child<TBaseElem>(elem, 0);
432 0 : if (sv.is_contained(child, dd->grid_level(), SurfaceView::SURFACE_RIM))
433 0 : continue;
434 : }
435 : }
436 :
437 0 : for (size_t fi = 0; fi < nFct; ++fi)
438 : {
439 0 : dd->inner_dof_indices(elem, fg[fi], vInd);
440 :
441 : const size_t nDof = vInd.size();
442 0 : for (size_t dof = 0; dof < nDof; ++dof)
443 : {
444 0 : const number tmp = DoFRef(uFine, vInd[dof]);
445 0 : sum += tmp*tmp;
446 : }
447 : }
448 : }
449 : // note: no duplicate additions possible
450 0 : }
451 : };
452 :
453 :
454 :
455 : //! Wrapper class for time dependence.
456 : template <typename TGridFunction>
457 : class TimeDependentSpace :
458 : public IGridFunctionSpace<TGridFunction>
459 : {
460 : public:
461 : typedef IGridFunctionSpace<TGridFunction> base_type;
462 : typedef IComponentSpace<TGridFunction> comp_space_type;
463 :
464 : /// time dependent CTOR
465 : TimeDependentSpace(SmartPtr<comp_space_type> spSpace, number tScale)
466 : : m_spSpatialSpace(spSpace), m_tScale(tScale) {};
467 :
468 : /// DTOR
469 : virtual ~TimeDependentSpace() {};
470 :
471 : using base_type::norm;
472 : using base_type::distance;
473 : using base_type::scaling;
474 :
475 : /// scaling (OVERRIDE)
476 : double scaling() const
477 : { return (m_spSpatialSpace->scaling()*m_tScale); }
478 :
479 : /// characteristic time
480 : void update_time_data(number tScale)
481 0 : { m_tScale = tScale; }
482 :
483 : /// print config string
484 : std::string config_string() const
485 : {
486 : std::stringstream ss;
487 : ss << "TimeDependentSpace for " << std::endl;
488 : ss << m_spSpatialSpace->config_string() << std::endl;
489 : return ss.str();
490 : }
491 :
492 : protected:
493 : SmartPtr<comp_space_type> m_spSpatialSpace;
494 : number m_tScale;
495 : };
496 :
497 :
498 : /** Evaluates difference between two grid functions in L2 norm */
499 : template <typename TGridFunction>
500 : class L2ComponentSpace :
501 : public IComponentSpace<TGridFunction>,
502 : public IObjectWithWeights<typename L2DistIntegrand<TGridFunction>::weight_type >
503 : {
504 : public:
505 : typedef IComponentSpace<TGridFunction> base_type;
506 : typedef typename L2DistIntegrand<TGridFunction>::weight_type weight_type;
507 : typedef IObjectWithWeights<weight_type> weighted_obj_type;
508 :
509 : /// CTOR
510 0 : L2ComponentSpace(const char *fctNames)
511 0 : : base_type(fctNames), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))) {};
512 :
513 0 : L2ComponentSpace(const char *fctNames, int order)
514 0 : : base_type(fctNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))) {};
515 :
516 0 : L2ComponentSpace(const char *fctNames, int order, double weight, const char* ssNames=0)
517 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(weight))) {};
518 :
519 0 : L2ComponentSpace(const char *fctNames, int order, ConstSmartPtr<weight_type> spWeight, const char* ssNames=0)
520 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight) {};
521 :
522 : /// DTOR
523 0 : ~L2ComponentSpace() {};
524 :
525 : using IComponentSpace<TGridFunction>::norm;
526 : using IComponentSpace<TGridFunction>::distance;
527 :
528 : /// for weighted norms
529 : using weighted_obj_type::set_weight;
530 : using weighted_obj_type::get_weight;
531 : using weighted_obj_type::m_spWeight;
532 :
533 : /// \copydoc IComponentSpace<TGridFunction>::norm
534 0 : double norm2(TGridFunction& uFine)
535 0 : { return L2Norm2(uFine, base_type::m_fctNames.c_str(), base_type::m_quadorder, base_type::m_ssNames, weighted_obj_type::m_spWeight); }
536 :
537 : /// \copydoc IComponentSpace<TGridFunction>::distance
538 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
539 0 : { return L2Distance2(uFine, base_type::m_fctNames.c_str(), uCoarse, base_type::m_fctNames.c_str(),
540 0 : base_type::m_quadorder, base_type::m_ssNames, weighted_obj_type::m_spWeight);}
541 :
542 :
543 :
544 : };
545 :
546 : /** Evaluates difference between two grid functions in L2 norm */
547 : template <typename TGridFunction>
548 : class L2QuotientSpace :
549 : public IComponentSpace<TGridFunction>,
550 : public IObjectWithWeights<typename L2DistIntegrand<TGridFunction>::weight_type >
551 : {
552 : public:
553 : typedef IComponentSpace<TGridFunction> base_type;
554 : typedef typename L2DistIntegrand<TGridFunction>::weight_type weight_type;
555 : typedef IObjectWithWeights<weight_type> weighted_obj_type;
556 :
557 : /// CTOR
558 0 : L2QuotientSpace(const char *fctNames)
559 0 : : base_type(fctNames), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))) {};
560 :
561 0 : L2QuotientSpace(const char *fctNames, int order)
562 0 : : base_type(fctNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))) {};
563 :
564 0 : L2QuotientSpace(const char *fctNames, int order, double weight, const char* ssNames=0)
565 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(weight))) {};
566 :
567 0 : L2QuotientSpace(const char *fctNames, int order, ConstSmartPtr<weight_type> spWeight, const char* ssNames=0)
568 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight) {};
569 :
570 : /// DTOR
571 0 : ~L2QuotientSpace() {};
572 :
573 : using IComponentSpace<TGridFunction>::norm;
574 : using IComponentSpace<TGridFunction>::distance;
575 :
576 : /// for weighted norms
577 : using weighted_obj_type::set_weight;
578 : using weighted_obj_type::get_weight;
579 : using weighted_obj_type::m_spWeight;
580 :
581 : /// \copydoc IComponentSpace<TGridFunction>::norm
582 0 : double norm2(TGridFunction& u)
583 : {
584 : typedef ConstUserNumber<TGridFunction::dim> MyConstUserData;
585 : typedef SmartPtr<UserData<number, TGridFunction::dim> > SPUserData;
586 :
587 0 : SPUserData spConst= make_sp(new MyConstUserData(1.0));
588 0 : number Meas = Integral(spConst, u, base_type::m_ssNames, 0.0, 1, "best");
589 0 : number uAvg = Integral(u, base_type::m_fctNames.c_str(), base_type::m_ssNames, base_type::m_quadorder);
590 :
591 0 : std::cerr << "Average:=" << uAvg <<"/" << Meas << " = " << uAvg/Meas << std::endl;
592 0 : SPUserData spAvg = make_sp(new MyConstUserData(uAvg/Meas));
593 0 : double qnorm = L2Error(spAvg, u, base_type::m_fctNames.c_str(), 0.0, base_type::m_quadorder, base_type::m_ssNames);
594 :
595 0 : return qnorm*qnorm;
596 : }
597 :
598 : /// \copydoc IComponentSpace<TGridFunction>::distance
599 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
600 : {
601 : typedef ConstUserNumber<TGridFunction::dim> MyConstUserData;
602 : typedef SmartPtr<UserData<number, TGridFunction::dim> > SPUserData;
603 :
604 0 : SPUserData spConst= make_sp(new MyConstUserData(1.0));
605 0 : number Meas = Integral(spConst, uFine, base_type::m_ssNames, 0.0, 1, "best");
606 0 : number avgFine = Integral(uFine, base_type::m_fctNames.c_str(), base_type::m_ssNames, base_type::m_quadorder);
607 0 : number avgCoarse = Integral(uCoarse, base_type::m_fctNames.c_str(), base_type::m_ssNames, base_type::m_quadorder);
608 :
609 0 : std::cerr << "Average:=(" << avgFine << "-" << avgCoarse<<")/" << Meas << " = " << (avgFine-avgCoarse)/Meas << std::endl;
610 0 : return L2Distance2(uFine, base_type::m_fctNames.c_str(),
611 : uCoarse, base_type::m_fctNames.c_str(),
612 : base_type::m_quadorder, base_type::m_ssNames, weighted_obj_type::m_spWeight,
613 0 : (avgFine-avgCoarse)/Meas);
614 : }
615 :
616 :
617 :
618 : };
619 :
620 :
621 : /** Evaluates distance between two grid functions in H1 semi-norm */
622 : template <typename TGridFunction>
623 : class H1SemiComponentSpace
624 : : public IComponentSpace<TGridFunction>,
625 : public IObjectWithWeights<typename H1SemiDistIntegrand<TGridFunction>::weight_type >
626 : {
627 : public:
628 : typedef IComponentSpace<TGridFunction> base_type;
629 : typedef typename H1SemiDistIntegrand<TGridFunction>::weight_type weight_type;
630 : typedef IObjectWithWeights<weight_type> weighted_obj_type;
631 :
632 :
633 0 : H1SemiComponentSpace(const char *fctNames)
634 0 : : base_type(fctNames), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(1.0))) {};
635 :
636 0 : H1SemiComponentSpace(const char *fctNames, int order)
637 0 : : base_type(fctNames, order), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(1.0))) {};
638 :
639 0 : H1SemiComponentSpace(const char *fctNames, int order, number weight, const char* ssNames=0)
640 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(weight))) {};
641 :
642 0 : H1SemiComponentSpace(const char *fctNames, int order, ConstSmartPtr<weight_type> spWeight, const char* ssNames=0)
643 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight) {};
644 :
645 0 : H1SemiComponentSpace(const char *fctNames, int order, const char* ssNames, ConstSmartPtr<weight_type> spWeight)
646 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight) {};
647 :
648 : /// DTOR
649 0 : ~H1SemiComponentSpace() {};
650 :
651 :
652 : /// \copydoc IComponentSpace<TGridFunction>::norm
653 : /// \copydoc IComponentSpace<TGridFunction>::distance
654 : using IComponentSpace<TGridFunction>::norm;
655 : using IComponentSpace<TGridFunction>::distance;
656 :
657 : /// \copydoc IComponentSpace<TGridFunction>::norm2
658 0 : double norm2(TGridFunction& uFine)
659 0 : { return H1SemiNorm2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), base_type::m_quadorder, NULL, weighted_obj_type::m_spWeight); }
660 :
661 : /// \copydoc IComponentSpace<TGridFunction>::distance2
662 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
663 0 : { return H1SemiDistance2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), uCoarse, base_type::m_fctNames.c_str(), base_type::m_quadorder, m_spWeight); }
664 :
665 : /// for weighted norms
666 : using weighted_obj_type::set_weight;
667 : using weighted_obj_type::get_weight;
668 : using weighted_obj_type::m_spWeight;
669 :
670 : };
671 :
672 :
673 :
674 : /*
675 : //! Evaluates distance w.r.t kinetic energy
676 : template <typename TGridFunction>
677 : class KineticEnergyComponentSpace
678 : : public IComponentSpace<TGridFunction>,
679 : public IObjectWithWeights<typename L2Integrand<TGridFunction>::weight_type >
680 : {
681 : public:
682 : typedef IComponentSpace<TGridFunction> base_type;
683 : typedef typename L2Integrand<TGridFunction>::weight_type weight_type;
684 : typedef IObjectWithWeights<weight_type> weighted_obj_type;
685 : typedef DarcyVelocityLinker<TGridFunction::dim> velocity_type;
686 :
687 :
688 : KineticEnergyComponentSpace(SmartPtr<velocity_type> spVelocity, const char *fctNames)
689 : : base_type(fctNames), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))), m_spVelocity(spVelocity) {};
690 :
691 : KineticEnergyComponentSpace(SmartPtr<velocity_type> spVelocity, const char *fctNames, int order)
692 : : base_type(fctNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(1.0))), m_spVelocity(spVelocity) {};
693 :
694 : KineticEnergyComponentSpace(SmartPtr<velocity_type> spVelocity, const char *fctNames, int order, number weight, const char* ssNames=0)
695 : : base_type(fctNames, ssNames, order), weighted_obj_type(make_sp(new ConstUserNumber<TGridFunction::dim>(weight))), m_spVelocity(spVelocity) {};
696 :
697 : KineticEnergyComponentSpace(SmartPtr<velocity_type> spVelocity, const char *fctNames, int order, ConstSmartPtr<weight_type> spWeight, const char* ssNames=0)
698 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight), m_spVelocity(spVelocity) {};
699 :
700 : KineticEnergyComponentSpace(SmartPtr<velocity_type> spVelocity, const char *fctNames, int order, const char* ssNames, ConstSmartPtr<weight_type> spWeight)
701 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight), m_spVelocity(spVelocity) {};
702 :
703 : /// DTOR
704 : ~KineticEnergyComponentSpace() {};
705 :
706 :
707 : using IComponentSpace<TGridFunction>::m_quadorder;
708 : using IComponentSpace<TGridFunction>::norm;
709 : using IComponentSpace<TGridFunction>::distance;
710 :
711 : /// \copydoc IComponentSpace<TGridFunction>::norm2
712 : double norm2(TGridFunction& u)
713 : {
714 : const char* subsets = NULL;
715 : UserDataIntegrandSq<MathVector<TGridFunction::dim>, TGridFunction> integrand2(m_spVelocity, &u, 0.0);
716 : return IntegrateSubsets(integrand2, u, subsets, m_quadorder);
717 : }
718 :
719 : /// \copydoc IComponentSpace<TGridFunction>::distance2
720 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
721 : {
722 : // UG_THROW("Not implemented!");
723 : const char* subsets = NULL;
724 : UserDataDistIntegrandSq<MathVector<TGridFunction::dim>, TGridFunction> integrand2(m_spVelocity, uFine, 0, uCoarse, 0);
725 : return IntegrateSubsets(integrand2, uFine, subsets, m_quadorder);
726 : }
727 :
728 : /// for weighted norms
729 : using weighted_obj_type::set_weight;
730 : using weighted_obj_type::get_weight;
731 : using weighted_obj_type::m_spWeight;
732 :
733 : void set_velocity(SmartPtr<velocity_type> spVelocity)
734 : { m_spVelocity = spVelocity;}
735 :
736 : protected:
737 : SmartPtr<velocity_type> m_spVelocity;
738 :
739 : };
740 : */
741 :
742 : //!
743 : /** Evaluates distance between two grid functions in H1 semi-norm */
744 : template <typename TGridFunction>
745 : class H1EnergyComponentSpace
746 : : public IComponentSpace<TGridFunction>,
747 : public IObjectWithWeights<typename H1EnergyDistIntegrand<TGridFunction>::weight_type >
748 : {
749 : public:
750 : typedef IComponentSpace<TGridFunction> base_type;
751 : typedef typename H1SemiDistIntegrand<TGridFunction>::weight_type weight_type;
752 : typedef IObjectWithWeights<weight_type> weighted_obj_type;
753 : // typedef DarcyVelocityLinker<TGridFunction::dim> velocity_type;
754 : static const int dim = TGridFunction::dim;
755 : typedef UserData<MathVector<dim>, dim> velocity_type;
756 :
757 0 : H1EnergyComponentSpace(const char *fctNames)
758 0 : : base_type(fctNames), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(1.0))), m_spVelocity(SPNULL) {};
759 :
760 0 : H1EnergyComponentSpace(const char *fctNames, int order)
761 0 : : base_type(fctNames, order), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(1.0))), m_spVelocity(SPNULL) {};
762 :
763 0 : H1EnergyComponentSpace(const char *fctNames, int order, number weight, const char* ssNames=0)
764 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(make_sp(new ConstUserMatrix<TGridFunction::dim>(weight))), m_spVelocity(SPNULL) {};
765 :
766 0 : H1EnergyComponentSpace(const char *fctNames, int order, ConstSmartPtr<weight_type> spWeight, const char* ssNames=0)
767 0 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight), m_spVelocity(SPNULL) {};
768 :
769 : /*H1EnergyComponentSpace(const char *fctNames, int order, const char* ssNames, ConstSmartPtr<weight_type> spWeight)
770 : : base_type(fctNames, ssNames, order), weighted_obj_type(spWeight), m_spVelocity(SPNULL) {};*/
771 :
772 : /// DTOR
773 0 : ~H1EnergyComponentSpace() {};
774 :
775 :
776 : /// \copydoc IComponentSpace<TGridFunction>::norm
777 : /// \copydoc IComponentSpace<TGridFunction>::distance
778 : using IComponentSpace<TGridFunction>::norm;
779 : using IComponentSpace<TGridFunction>::distance;
780 :
781 : /// \copydoc IComponentSpace<TGridFunction>::norm2
782 0 : double norm2(TGridFunction& uFine)
783 : {
784 0 : if (m_spVelocity.valid()) {
785 : //const char* subsets = NULL; // [q^2]
786 0 : UserDataIntegrandSq<MathVector<TGridFunction::dim>, TGridFunction> integrand2(m_spVelocity, &uFine, 0.0);
787 0 : return IntegrateSubsets(integrand2, uFine, base_type::m_ssNames, base_type::m_quadorder);
788 : } else {
789 0 : return H1EnergyNorm2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), base_type::m_quadorder, base_type::m_ssNames, weighted_obj_type::m_spWeight);
790 : }
791 : }
792 :
793 : /// \copydoc IComponentSpace<TGridFunction>::distance2
794 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
795 0 : { return H1EnergyDistance2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), uCoarse, base_type::m_fctNames.c_str(), base_type::m_quadorder,base_type::m_ssNames, weighted_obj_type::m_spWeight); }
796 :
797 : /// for weighted norms
798 : using weighted_obj_type::set_weight;
799 : using weighted_obj_type::get_weight;
800 :
801 0 : void set_velocity(SmartPtr<velocity_type> spVelocity)
802 0 : { m_spVelocity = spVelocity;}
803 :
804 : protected:
805 : SmartPtr<velocity_type> m_spVelocity;
806 :
807 : };
808 :
809 : /** Defines a H1 space for single component */
810 : template <typename TGridFunction>
811 : class H1ComponentSpace :
812 : public IComponentSpace<TGridFunction>
813 : {
814 : public:
815 : typedef IComponentSpace<TGridFunction> base_type;
816 :
817 0 : H1ComponentSpace(const char *fctNames) : base_type(fctNames) {};
818 0 : H1ComponentSpace(const char *fctNames, int order) : base_type(fctNames, order) {};
819 0 : H1ComponentSpace(const char *fctNames, const char* ssNames, int order) : base_type(fctNames, ssNames, order) {};
820 0 : ~H1ComponentSpace() {};
821 :
822 : using IComponentSpace<TGridFunction>::norm;
823 : using IComponentSpace<TGridFunction>::distance;
824 :
825 : /// \copydoc IComponentSpace<TGridFunction>::norm
826 0 : double norm2(TGridFunction& uFine)
827 0 : { return H1Norm2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), base_type::m_quadorder, base_type::m_ssNames); }
828 :
829 : /// \copydoc IComponentSpace<TGridFunction>::norm
830 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
831 0 : { return H1Distance2<TGridFunction>(uFine, base_type::m_fctNames.c_str(), uCoarse, base_type::m_fctNames.c_str(), base_type::m_quadorder, base_type::m_ssNames); }
832 :
833 : };
834 :
835 :
836 :
837 : //! Defines a composite space, (i.e., additive composition from other spaces)
838 : /*! Employs a l2-type extension: \| u \|^2 := \sum_I \sigma_I \| u \|_I^2
839 : * TODO: Merge with ApproximationSpace?
840 : * */
841 : template <typename TGridFunction>
842 : class CompositeSpace : public IGridFunctionSpace<TGridFunction>
843 : {
844 : public:
845 : typedef IGridFunctionSpace<TGridFunction> base_type;
846 : typedef IComponentSpace<TGridFunction> obj_type;
847 : typedef TimeDependentSpace<TGridFunction> time_dependent_obj_type;
848 : typedef std::pair<SmartPtr<obj_type>, number> weighted_obj_type;
849 :
850 0 : CompositeSpace() {};
851 : // virtual ~CompositeSpace() {};
852 :
853 : using base_type::norm;
854 : using base_type::distance;
855 :
856 : /// \copydoc IComponentSpace<TGridFunction>::norm
857 0 : double norm(TGridFunction& uFine)
858 0 : { return(sqrt(norm2(uFine))); }
859 :
860 : /// \copydoc IComponentSpace<TGridFunction>::norm2
861 0 : double norm2(TGridFunction& uFine)
862 : {
863 : number unorm2 = 0.0;
864 0 : for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin();
865 0 : it!= m_spWeightedSubspaces.end(); ++it)
866 : {
867 0 : double snorm2 = it->first->norm2(uFine);
868 0 : unorm2 += it->second * snorm2; // scaling
869 0 : UG_LOG("composite-norm2:\t" << snorm2 << "\t*\t" << it->second
870 : << "\t=\t" << it->second * snorm2 << std::endl);
871 : }
872 : UG_LOG("composite-norm2-final:\t" << unorm2 << std::endl);
873 0 : return unorm2;
874 : }
875 :
876 : /// \copydoc IComponentSpace<TGridFunction>::distance2
877 0 : double distance2(TGridFunction& uFine, TGridFunction& uCoarse)
878 : {
879 : number unorm2 = 0.0;
880 0 : for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin();
881 0 : it!= m_spWeightedSubspaces.end(); ++it)
882 : {
883 0 : double sdist2 = it->first->distance2(uFine, uCoarse);
884 0 : unorm2 += it->second * sdist2; // scaling
885 0 : UG_LOG("composite-dist2:\t" << sdist2 << "\t*\t" << it->second
886 : << "\t=\t" << it->second * sdist2 << std::endl);
887 : }
888 : UG_LOG("composite-dist2-final:\t" << unorm2 << std::endl);
889 0 : return unorm2;
890 : }
891 :
892 : /// \copydoc IComponentSpace<TGridFunction>::distance2
893 0 : double distance(TGridFunction& uFine, TGridFunction& uCoarse)
894 0 : { return sqrt(distance2(uFine, uCoarse)); }
895 :
896 : /// add space to composite (with weight 1.0)
897 0 : void add(SmartPtr<obj_type> spSubSpace)
898 0 : { m_spWeightedSubspaces.push_back(std::make_pair(spSubSpace, 1.0)); }
899 :
900 : /// add space to composite (with variable weight)
901 0 : void add(SmartPtr<obj_type> spSubSpace, number sigma)
902 0 : { m_spWeightedSubspaces.push_back(std::make_pair(spSubSpace, sigma)); }
903 :
904 : /// print config string
905 0 : std::string config_string() const
906 : {
907 0 : std::stringstream ss;
908 : ss << "CompositeSpace:" << std::endl;
909 :
910 0 : for (typename std::vector<weighted_obj_type>::const_iterator it = m_spWeightedSubspaces.begin();
911 0 : it!= m_spWeightedSubspaces.end(); ++it)
912 0 : { ss << it->first->config_string(); }
913 :
914 0 : return ss.str();
915 0 : }
916 :
917 :
918 : //! Forward update to all members
919 0 : void update_time_data(number t)
920 : {
921 0 : for (typename std::vector<weighted_obj_type>::iterator it = m_spWeightedSubspaces.begin();
922 0 : it!= m_spWeightedSubspaces.end(); ++it)
923 : {
924 0 : SmartPtr<time_dependent_obj_type> spSpaceT = it->first.template cast_dynamic<time_dependent_obj_type> ();
925 0 : if (spSpaceT.valid()) spSpaceT->update_time_data(t);
926 : }
927 0 : }
928 :
929 : //! Check, if any object is time-dependent.
930 0 : bool is_time_dependent() const
931 : {
932 0 : for (typename std::vector<weighted_obj_type>::const_iterator it = m_spWeightedSubspaces.begin();
933 0 : it!= m_spWeightedSubspaces.end(); ++it)
934 : {
935 0 : SmartPtr<time_dependent_obj_type> spSpaceT = it->first.template cast_dynamic<time_dependent_obj_type>();
936 0 : if (spSpaceT.valid()) return true;
937 : }
938 : return false;
939 : }
940 :
941 : const std::vector<weighted_obj_type> &get_subspaces() const
942 : { return m_spWeightedSubspaces; }
943 :
944 : protected:
945 : std::vector<weighted_obj_type> m_spWeightedSubspaces;
946 :
947 : };
948 :
949 : /*
950 : template <typename TGridFunction>
951 : class TimeDependentCompositeSpace
952 : : public CompositeSpace<TGridFunction>, public ITimeDependentSpace<TGridFunction>
953 : {
954 : protected:
955 : using CompositeSpace<TGridFunction>::m_spSubspaces;
956 :
957 : public:
958 : typedef typename IComponentSpace<TGridFunction> obj_type;
959 :
960 :
961 : /// add spaces to composite
962 : void add(SmartPtr<obj_type> spSubSpace)
963 : { m_spSubspaces.push_back(make_sp(new ITimeDependentSpace(spSubSpace))); }
964 :
965 : //! Forward update to all members
966 : void update_time_data(number t)
967 : {
968 : for (typename std::vector<SmartPtr<obj_type> >::iterator it = m_spSubspaces.begin();
969 : it!= m_spSubspaces.end(); ++it)
970 : { if ((*it)->is_time_dependent()) (*it)->update_time_data(t); }
971 : }
972 :
973 : //! Check, if all objects are time-dependent.
974 : bool is_time_dependent() const
975 : {
976 : for (typename std::vector<SmartPtr<obj_type> >::const_iterator it = m_spSubspaces.begin();
977 : it!= m_spSubspaces.end(); ++it)
978 : { if ((*it)->is_time_dependent()) return true; }
979 : return false;
980 : }
981 : };
982 :
983 : */
984 :
985 : } // namespace ug
986 : #endif
|
