ug4
Loading...
Searching...
No Matches
error_indicator.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt
3 * Authors: Andreas Vogel, Sebastian Reiter
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
34#ifndef __H__UG__LIB_DISC__FUNCTION_SPACE__ERROR_INDICATOR__
35#define __H__UG__LIB_DISC__FUNCTION_SPACE__ERROR_INDICATOR__
36
37#include <vector>
38
40#include "gradient_evaluators.h"
41#include "common/common.h"
48#include "integrate.h"
50
51#ifdef UG_PARALLEL
53#endif
54
55namespace ug{
56
57
58template <typename TFunction>
59void ComputeGradientLagrange1(TFunction& u, size_t fct,
61 typename TFunction::element_type,
62 ug::Attachment<number> >& aaError)
63{
64 static const int dim = TFunction::dim;
65 typedef typename TFunction::const_element_iterator const_iterator;
66 typedef typename TFunction::element_type element_type;
67
68// get position accessor
69 typename TFunction::domain_type::position_accessor_type& aaPos
70 = u.domain()->position_accessor();
71
72// some storage
74 std::vector<MathVector<dim> > vLocalGrad;
75 std::vector<MathVector<dim> > vGlobalGrad;
76 std::vector<MathVector<dim> > vCorner;
77
78// get iterator over elements
79 const_iterator iter = u.template begin<element_type>();
80 const_iterator iterEnd = u.template end<element_type>();
81
82// loop elements
83 for(; iter != iterEnd; ++iter)
84 {
85 // get the element
86 element_type* elem = *iter;
87
88 // reference object type
89 ReferenceObjectID roid = elem->reference_object_id();
90
91 // get trial space
92 const LocalShapeFunctionSet<dim>& lsfs =
93 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
94
95 // create a reference mapping
97 = ReferenceMappingProvider::get<dim, dim>(roid);
98
99 // get local Mid Point
100 MathVector<dim> localIP = ReferenceElementCenter<dim>(roid);
101
102 // number of shape functions
103 const size_t numSH = lsfs.num_sh();
104 vLocalGrad.resize(numSH);
105 vGlobalGrad.resize(numSH);
106
107 // evaluate reference gradient at local midpoint
108 lsfs.grads(&vLocalGrad[0], localIP);
109
110 // get corners of element
111 CollectCornerCoordinates(vCorner, *elem, aaPos);
112
113 // update mapping
114 map.update(&vCorner[0]);
115
116 // compute jacobian
117 map.jacobian_transposed_inverse(JTInv, localIP);
118
119 // compute size (volume) of element
120 const number elemSize = ElementSize<dim>(roid, &vCorner[0]);
121
122 // compute gradient at mid point by summing contributions of all shape fct
123 MathVector<dim> MidGrad; VecSet(MidGrad, 0.0);
124 for(size_t sh = 0 ; sh < numSH; ++sh)
125 {
126 // get global Gradient
127 MatVecMult(vGlobalGrad[sh], JTInv, vLocalGrad[sh]);
128
129 // get vertex
130 Vertex* vert = elem->vertex(sh);
131
132 // get of of vertex
133 std::vector<DoFIndex> ind;
134 u.inner_dof_indices(vert, fct, ind);
135
136 // scale global gradient
137 vGlobalGrad[sh] *= DoFRef(u, ind[0]);
138
139 // sum up
140 MidGrad += vGlobalGrad[sh];
141 }
142
143 // write result in array storage
144 aaError[elem] = VecTwoNorm(MidGrad) * pow(elemSize, 2./dim);
145 }
146}
147
148template <typename TFunction>
149void ComputeGradientCrouzeixRaviart(TFunction& u, size_t fct,
151 typename TFunction::element_type,
152 ug::Attachment<number> >& aaError)
153{
154 static const int dim = TFunction::dim;
155 typedef typename TFunction::const_element_iterator const_iterator;
156 typedef typename TFunction::domain_type domain_type;
157 typedef typename domain_type::grid_type grid_type;
158 typedef typename TFunction::element_type element_type;
159 typedef typename element_type::side side_type;
160
161 typename grid_type::template traits<side_type>::secure_container sides;
162
163// get position accessor
164 typename TFunction::domain_type::position_accessor_type& aaPos
165 = u.domain()->position_accessor();
166
167 grid_type& grid = *u.domain()->grid();
168
169// some storage
171 std::vector<MathVector<dim> > vLocalGrad;
172 std::vector<MathVector<dim> > vGlobalGrad;
173 std::vector<MathVector<dim> > vCorner;
174
175// get iterator over elements
176 const_iterator iter = u.template begin<element_type>();
177 const_iterator iterEnd = u.template end<element_type>();
178
179// loop elements
180 for(; iter != iterEnd; ++iter)
181 {
182 // get the element
183 element_type* elem = *iter;
184
185 // get sides of element
186 grid.associated_elements_sorted(sides, elem );
187
188 // reference object type
189 ReferenceObjectID roid = elem->reference_object_id();
190
191 // get trial space
192 const LocalShapeFunctionSet<dim>& lsfs =
193 LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART, dim, 1));
194
195 // create a reference mapping
197 = ReferenceMappingProvider::get<dim, dim>(roid);
198
199 // get local Mid Point
200 MathVector<dim> localIP = ReferenceElementCenter<dim>(roid);
201
202 // number of shape functions
203 const size_t numSH = lsfs.num_sh();
204 vLocalGrad.resize(numSH);
205 vGlobalGrad.resize(numSH);
206
207 // evaluate reference gradient at local midpoint
208 lsfs.grads(&vLocalGrad[0], localIP);
209
210 // get corners of element
211 CollectCornerCoordinates(vCorner, *elem, aaPos);
212
213 // update mapping
214 map.update(&vCorner[0]);
215
216 // compute jacobian
217 map.jacobian_transposed_inverse(JTInv, localIP);
218
219 // compute size (volume) of element
220 const number elemSize = ElementSize<dim>(roid, &vCorner[0]);
221
222 // compute gradient at mid point by summing contributions of all shape fct
223 MathVector<dim> MidGrad; VecSet(MidGrad, 0.0);
224 for(size_t sh = 0 ; sh < numSH; ++sh)
225 {
226 // get global Gradient
227 MatVecMult(vGlobalGrad[sh], JTInv, vLocalGrad[sh]);
228
229 // get of of vertex
230 std::vector<DoFIndex> ind;
231 u.inner_dof_indices(sides[sh], fct, ind);
232
233 // scale global gradient
234 vGlobalGrad[sh] *= DoFRef(u, ind[0]);
235
236 // sum up
237 MidGrad += vGlobalGrad[sh];
238 }
239
240 // write result in array storage
241 aaError[elem] = VecTwoNorm(MidGrad) * pow(elemSize, 2./dim);
242 }
243}
244
245template <int dim> struct face_type_traits
246{
247 typedef void face_type0;
248 typedef void face_type1;
249};
250
251template <> struct face_type_traits<1>
252{
255};
256
257template <> struct face_type_traits<2>
258{
261};
262
263template <> struct face_type_traits<3>
264{
267};
268
269template <typename TFunction>
270void ComputeGradientPiecewiseConstant(TFunction& u, size_t fct,
272 typename TFunction::element_type,
273 ug::Attachment<number> >& aaError)
274{
275 static const int dim = TFunction::dim;
276 typedef typename TFunction::const_element_iterator const_iterator;
277 typedef typename TFunction::domain_type domain_type;
278 typedef typename domain_type::grid_type grid_type;
279 typedef typename TFunction::element_type element_type;
280 typedef typename element_type::side side_type;
281
282 typename grid_type::template traits<side_type>::secure_container sides;
283
284 typedef typename face_type_traits<dim>::face_type0 face_type0;
285 typedef typename face_type_traits<dim>::face_type1 face_type1;
286
287// get position accessor
288 typename TFunction::domain_type::position_accessor_type& aaPos
289 = u.domain()->position_accessor();
290
291 grid_type& grid = *u.domain()->grid();
292
293// some storage
295
296 std::vector<MathVector<dim> > vCorner;
297 std::vector<MathVector<dim> > sideCorners;
298 std::vector<MathVector<dim> > vSideCoPos;
299
300// get iterator over elements
301 const_iterator iter = u.template begin<element_type>();
302 const_iterator iterEnd = u.template end<element_type>();
303
304// loop elements
305 for(; iter != iterEnd; ++iter)
306 {
307 // get the element
308 element_type* elem = *iter;
309
310 MathVector<dim> vGlobalGrad=0;
311
312 // get sides of element
313 grid.associated_elements_sorted(sides, elem );
314
315 // reference object type
316 ReferenceObjectID roid = elem->reference_object_id();
317
318 const DimReferenceElement<dim>& rRefElem
319 = ReferenceElementProvider::get<dim>(roid);
320
321 // get corners of element
322 CollectCornerCoordinates(vCorner, *elem, aaPos);
323
324 // compute size (volume) of element
325 const number elemSize = ElementSize<dim>(roid, &vCorner[0]);
326
327 typename grid_type::template traits<element_type>::secure_container assoElements;
328
329 // assemble element-wise finite volume gradient
330 for (size_t s=0;s<sides.size();s++){
331 grid.associated_elements(assoElements,sides[s]);
332 // face value is average of associated elements
333 number faceValue = 0;
334 size_t numOfAsso = assoElements.size();
335 for (size_t i=0;i<numOfAsso;i++){
336 std::vector<DoFIndex> ind;
337 u.inner_dof_indices(assoElements[i], fct, ind);
338 faceValue+=DoFRef(u, ind[0]);
339 }
340 faceValue/=(number)numOfAsso;
341 MathVector<dim> normal;
342 size_t numSideCo = rRefElem.num(dim-1,s,0);
343
344 for (size_t i = 0; i < numSideCo; ++i)
345 vSideCoPos.push_back(vCorner[rRefElem.id(dim-1, s, 0, i)]);
346
347 // faces have dim corners in 1d, 2d
348 // in 3d they have dim corners (triangle) or dim+1 corners (quadrilateral)
349 if ((int)numSideCo==dim)
350 ElementNormal<face_type0,dim>(normal,&vSideCoPos[0]);
351 else
352 ElementNormal<face_type1,dim>(normal,&vSideCoPos[0]);
353
354 for (int d=0;d<dim;d++){
355 vGlobalGrad[d] += faceValue * normal[d];
356 }
357 }
358 vGlobalGrad/=(number)elemSize;
359
360 // write result in array storage
361 aaError[elem] = VecTwoNorm(vGlobalGrad) * pow(elemSize, 2./dim);
362 }
363}
364
365
366template <typename TDomain, typename TAlgebra>
369 const char* fctName,
370 number TOL,
371 number refineFrac, number coarseFrac,
372 int maxLevel)
373{
374 PROFILE_FUNC();
375// types
376 typedef GridFunction<TDomain, TAlgebra> TFunction;
377 typedef typename TFunction::domain_type::grid_type grid_type;
378 typedef typename TFunction::element_type element_type;
379 const int dim = TFunction::dim;
380
381// function id
382 const size_t fct = u.fct_id_by_name(fctName);
383
384// get multigrid
385 SmartPtr<grid_type> pMG = u.domain()->grid();
386
387// attach error field
389 ANumber aError;
390 pMG->template attach_to<element_type>(aError);
392
393// Compute error on elements
394 if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
395 ComputeGradientLagrange1(u, fct, aaError);
396 else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
397 ComputeGradientCrouzeixRaviart(u, fct, aaError);
398 else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
400 else{
401 UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
402 }
403
404// Mark elements for refinement
405 MarkElements<element_type> (aaError, refiner, u.dof_distribution(), TOL, refineFrac, coarseFrac, maxLevel);
406
407// detach error field
408 pMG->template detach_from<element_type>(aError);
409};
410
411
412template <typename TDomain, typename TAlgebra>
415 const char* fctName,
416 number refTol,
417 number coarsenTol,
418 int minLvl,
419 int maxLevel)
420{
421 PROFILE_FUNC();
422// types
423 typedef GridFunction<TDomain, TAlgebra> TFunction;
424 typedef typename TFunction::domain_type::grid_type grid_type;
425 typedef typename TFunction::element_type element_type;
426 const int dim = TFunction::dim;
427
428// function id
429 const size_t fct = u.fct_id_by_name(fctName);
430
431// get multigrid
432 SmartPtr<grid_type> pMG = u.domain()->grid();
433
434// attach error field
436 ANumber aError;
437 pMG->template attach_to<element_type>(aError);
439
440// Compute error on elements
441 if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
442 ComputeGradientLagrange1(u, fct, aaError);
443 else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
444 ComputeGradientCrouzeixRaviart(u, fct, aaError);
445 else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
447 else{
448 UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
449 }
450
451// Mark elements for refinement
452 MarkElementsAbsolute<element_type> (aaError, refiner, u.dof_distribution (), refTol, coarsenTol, minLvl, maxLevel);
453
454// detach error field
455 pMG->template detach_from<element_type>(aError);
456};
457
458
459template <typename TFunction>
460void computeGradientJump(TFunction& u,
462 typename TFunction::element_type,
463 ug::Attachment<number> >& aaGrad,
465 typename TFunction::element_type,
466 ug::Attachment<number> >& aaError)
467{
468 typedef typename TFunction::domain_type domain_type;
469 typedef typename domain_type::grid_type grid_type;
470 typedef typename TFunction::element_type element_type;
471 typedef typename element_type::side side_type;
472 typedef typename TFunction::template traits<side_type>::const_iterator side_iterator;
473
474 grid_type& grid = *u.domain()->grid();
475
476 // get iterator over elements
477 side_iterator iter = u.template begin<side_type>();
478 side_iterator iterEnd = u.template end<side_type>();
479 // loop elements
480 for(; iter != iterEnd; ++iter)
481 {
482 // get the element
483 side_type* side = *iter;
484 typename grid_type::template traits<element_type>::secure_container neighElements;
485 grid.associated_elements(neighElements,side);
486 if (neighElements.size()!=2) continue;
487 number localJump = std::abs(aaGrad[neighElements[0]]-aaGrad[neighElements[1]]);
488 for (size_t i=0;i<2;i++)
489 if (aaError[neighElements[i]]<localJump) aaError[neighElements[i]]=localJump;
490 }
491}
492
493// indicator is based on elementwise computation of jump between elementwise gradients
494// the value in an element is the maximum jump between the element gradient and gradient
495// in elements with common face
496template <typename TDomain, typename TAlgebra>
499 const char* fctName,
500 number TOL,
501 number refineFrac, number coarseFrac,
502 int maxLevel)
503{
504 PROFILE_FUNC();
505// types
506 typedef GridFunction<TDomain, TAlgebra> TFunction;
507 typedef typename TFunction::domain_type::grid_type grid_type;
508 typedef typename TFunction::element_type element_type;
509 const int dim = TFunction::dim;
510
511// function id
512 const size_t fct = u.fct_id_by_name(fctName);
513
514// get multigrid
515 SmartPtr<grid_type> pMG = u.domain()->grid();
516
517// attach error field
519 ANumber aGrad;
520 pMG->template attach_to<element_type>(aGrad);
522
523 ANumber aError;
524 pMG->template attach_to<element_type>(aError);
526
527// Compute error on elements
528 if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
529 ComputeGradientLagrange1(u, fct, aaGrad);
530 else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
531 ComputeGradientCrouzeixRaviart(u, fct, aaGrad);
532 else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
534 else{
535 UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
536 }
537 computeGradientJump(u,aaGrad,aaError);
538
539// Mark elements for refinement
540 MarkElements<element_type> (aaError, refiner, u.dof_distribution(), TOL, refineFrac, coarseFrac, maxLevel);
541
542// detach error field
543 pMG->template detach_from<element_type>(aError);
544};
545
546template <typename TDomain, typename TAlgebra>
549 const char* fctName,
550 number refTol,
551 number coarsenTol,
552 int minLvl,
553 int maxLevel)
554{
555 PROFILE_FUNC();
556// types
557 typedef GridFunction<TDomain, TAlgebra> TFunction;
558 typedef typename TFunction::domain_type::grid_type grid_type;
559 typedef typename TFunction::element_type element_type;
560 const int dim = TFunction::dim;
561
562// function id
563 const size_t fct = u.fct_id_by_name(fctName);
564
565// get multigrid
566 SmartPtr<grid_type> pMG = u.domain()->grid();
567
568// attach error field
570 ANumber aGrad;
571 pMG->template attach_to<element_type>(aGrad);
573
574 ANumber aError;
575 pMG->template attach_to<element_type>(aError);
577
578// Compute error on elements
579 if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
580 ComputeGradientLagrange1(u, fct, aaGrad);
581 else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
582 ComputeGradientCrouzeixRaviart(u, fct, aaGrad);
583 else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
585 else{
586 UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
587 }
588 computeGradientJump(u,aaGrad,aaError);
589
590// Mark elements for refinement
591 MarkElementsAbsolute<element_type> (aaError, refiner, u.dof_distribution(), refTol, coarsenTol, minLvl, maxLevel);
592
593// detach error field
594 pMG->template detach_from<element_type>(aError);
595};
596
600template <typename TDomain, typename TAlgebra>
604 const char* cmp,
605 number minL2Error,
606 number maxL2Error,
607 number refFrac,
608 int minLvl, int maxLvl,
609 number time, int quadOrder)
610{
611 PROFILE_FUNC();
612 using namespace std;
613// types
614 typedef GridFunction<TDomain, TAlgebra> TFunction;
615 typedef typename TFunction::domain_type::grid_type grid_t;
616 typedef typename TFunction::element_type elem_t;
617 const int dim = TFunction::dim;
618
619// function id
620 const size_t fct = u->fct_id_by_name(cmp);
621 UG_COND_THROW(fct >= u->num_fct(),
622 "Function space does not contain a function with name " << cmp);
623
624// get multigrid and position accessor
625 grid_t& mg = *u->domain()->grid();
626
627 typename TFunction::domain_type::position_accessor_type&
628 aaPos = u->domain()->position_accessor();
629
630// attach error field
632 ANumber aError;
633 mg.template attach_to<elem_t>(aError);
635
636// create integrand and perform integration
637 /*SmartPtr<IIntegrand<number, dim> > spIntegrand
638 = make_sp(new L2ErrorIntegrand<TFunction>(spExactSol, u, fct, time));*/
639 L2ErrorIntegrand<TFunction> integrand(spExactSol, *u, fct, time);
640 number l2Error =
641 Integrate<dim, dim>(u->template begin<elem_t>(), u->template end<elem_t>(),
642 aaPos, integrand, quadOrder, "best", &aaError);
643
644 #ifdef UG_PARALLEL
646 l2Error = com.allreduce(l2Error, PCL_RO_SUM);
647 #endif
648
649 l2Error = sqrt(l2Error);
650
651 UG_LOG("maxError " << maxL2Error << ", l2Error " << l2Error << std::endl);
652
653 if(l2Error > maxL2Error){
654 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
655 size_t numElemsActive = 0; // number of elements which may be refined
656 size_t numElemsTotal = 0; // total number of elements
657 number maxElemError = 0; // error in elements which may be refined
658 number minElemError = numeric_limits<number>::max();
659 number fixedError = 0; // error in elements which can't be refined any further
660 for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
661 ++numElemsTotal;
662 if(mg.get_level(*iter) < maxLvl){
663 ++numElemsActive;
664 maxElemError = max(maxElemError, aaError[*iter]);
665 minElemError = min(minElemError, aaError[*iter]);
666 }
667 else{
668 fixedError += aaError[*iter];
669 }
670 }
671
672 #ifdef UG_PARALLEL
674 minElemError = com.allreduce(minElemError, PCL_RO_MIN);
675 maxElemError = com.allreduce(maxElemError, PCL_RO_MAX);
676 fixedError = com.allreduce(fixedError, PCL_RO_MAX);
677 #endif
678
679 // note that aaError contains error-squares
680 maxElemError = minElemError + (maxElemError - minElemError) * sq(refFrac);
681 number refThreshold = maxElemError;
682 if(maxL2Error > 0){
683 // here we try to reduce the number of elements marked in the last steps.
684 // note that each element stores error-squares...
685 //refThreshold = max(maxElemError, (sq(maxL2Error) - fixedError) / (number)numElemsActive);
686 //refThreshold = max(maxElemError, sq(maxL2Error)/ (number)numElemsTotal);
687 }
688
689 MarkElementsAbsolute(aaError, refiner, u->dof_distribution(), refThreshold, -1,
690 minLvl, maxLvl);
691 }
692
693// coarsening
694 // number coarsenThreshold = -1;
695 // if(minL2Error > 0)
696 // coarsenThreshold = minL2Error * minL2Error / (number)numElemsActive;
697
698// detach error field
699 mg.template detach_from<elem_t>(aError);
700};
701
702
713template <class side_t, class TFunction>
714void ExchangeAndAdjustSideErrors(TFunction& u, ANumber aSideError, ANumber aNumElems)
715{
716 //typedef typename TFunction::template traits<side_t>::const_iterator side_iter_t;
717 typedef typename TFunction::domain_type::grid_type grid_t;
718
719 grid_t& g = *u.domain()->grid();
720 Grid::AttachmentAccessor<side_t, ANumber> aaSideError(g, aSideError);
721 Grid::AttachmentAccessor<side_t, ANumber> aaNumElems(g, aNumElems);
722
723// in a parallel environment we now have to sum the gradients over parallel
724// interfaces
725// since there may be constrained sides which locally do not have a constraining
726// side we do this before adding the constrained values to their constraining objects.
727 #ifdef UG_PARALLEL
728 typedef typename GridLayoutMap::Types<side_t>::Layout layout_t;
729 DistributedGridManager& dgm = *g.distributed_grid_manager();
730 GridLayoutMap& glm = dgm.grid_layout_map();
732
733 // sum all copies at the h-master attachment
734 ComPol_AttachmentReduce<layout_t, ANumber> compolSumErr(g, aSideError, PCL_RO_SUM);
735 ComPol_AttachmentReduce<layout_t, ANumber> compolSumNum(g, aNumElems, PCL_RO_SUM);
736 icom.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolSumErr);
737 icom.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolSumNum);
738 icom.communicate();
739
740 // copy the sum from the master to all of its slave-copies
741 ComPol_CopyAttachment<layout_t, ANumber> compolCopyErr(g, aSideError);
742 ComPol_CopyAttachment<layout_t, ANumber> compolCopyNum(g, aNumElems);
743 icom.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolCopyErr);
744 icom.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolCopyNum);
745 icom.communicate();
746
747 // since we're copying from vmasters to vslaves later on, we have to make
748 // sure, that also all v-masters contain the correct values.
749 //todo: communication to vmasters may not be necessary here...
750 // it is performed to make sure that all surface-rim-children
751 // contain their true value.
752 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyErr);
753 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyNum);
754
755 icom.communicate();
756 #endif
757
758// if we're currently operating on a surface function, we have to adjust
759// errors between shadowed and shadowing sides (constraining and constrained sides).
760 if(u.grid_level().type() == GridLevel::SURFACE){
761 //todo: avoid iteration over the whole grid!
762 //todo: reduce communication
763 const SurfaceView* surfView = u.approx_space()->surface_view().get();
764
765 // for(side_iter_t iter = u.template begin<side_t>(SurfaceView::SHADOW_RIM);
766 // iter != u.template end<side_t>(SurfaceView::SHADOW_RIM); ++iter)
767 typedef typename Grid::traits<side_t>::iterator grid_side_iter_t;
768 for(grid_side_iter_t iter = g.template begin<side_t>();
769 iter != g.template end<side_t>(); ++iter)
770 {
771 side_t* s = *iter;
773 continue;
774
775 const size_t numChildren = g.template num_children<side_t>(s);
776 if(numChildren > 0){
777 number w = 1. / number(numChildren);
778 for(size_t i_child = 0; i_child < numChildren; ++i_child){
779 side_t* c = g.template get_child<side_t>(s, i_child);
780 aaSideError[s] += w * aaSideError[c];
781 aaNumElems[s] += w * aaNumElems[c];
782 }
783 }
784 }
785 // those elems which have local children now contain their true value (vslaves...)
786 #ifdef UG_PARALLEL
787 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyErr);
788 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyNum);
789 icom.communicate();
790 #endif
791
792 for(grid_side_iter_t iter = g.template begin<side_t>();
793 iter != g.template end<side_t>(); ++iter)
794 {
795 side_t* s = *iter;
797 continue;
798
799 const size_t numChildren = g.template num_children<side_t>(s);
800 if(numChildren > 0){
801 number w = 1. / number(numChildren);
802 for(size_t i_child = 0; i_child < numChildren; ++i_child){
803 side_t* c = g.template get_child<side_t>(s, i_child);
804 aaSideError[c] = w * aaSideError[s];
805 aaNumElems[c] = w * aaNumElems[s];
806 }
807 }
808 }
809
810 // those elems which have a local parent now contain their true value (vmasters...)
811 #ifdef UG_PARALLEL
812 // copy from v-masters to v-slaves, since there may be constrained
813 // sides which locally do not have a constraining element. Note that
814 // constrained V-Masters always have a local constraining element and thus
815 // contain the correct value.
816 icom.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, compolCopyErr);
817 icom.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, compolCopyNum);
818 icom.communicate();
819 #endif
820 }
821}
822
826template <typename TGradientEvaluator, typename TFunction>
827void EvaluateGradientJump_SideIntegral(TFunction& u, size_t fct,
829 typename TFunction::element_type,
830 ug::Attachment<number> >& aaError,
831 bool addErrSquareToAAError = false)
832{
833 using std::max;
834
835 static const int dim = TFunction::dim;
836 typedef typename TFunction::domain_type::grid_type grid_t;
837 typedef typename TFunction::const_element_iterator const_iterator;
838 typedef typename TFunction::element_type elem_t;
839 typedef typename elem_t::side side_t;
840 typedef MathVector<dim> vector_t;
841
842// get position accessor
843 typename TFunction::domain_type::position_accessor_type& aaPos
844 = u.domain()->position_accessor();
845
846// we need an attachment on the sides to gather gradient-jumps in parallel.
847// Note that a serial version could be created without this additional attachment.
848 grid_t& g = *u.domain()->grid();
849 ANumber aSideError;
850 ANumber aNumElems;
851 g.template attach_to_dv<side_t>(aSideError, 0);
852 g.template attach_to_dv<side_t>(aNumElems, 0);
853 Grid::AttachmentAccessor<side_t, ANumber> aaSideError(g, aSideError);
854 Grid::AttachmentAccessor<side_t, ANumber> aaNumElems(g, aNumElems);
855
856// loop elements
857 TGradientEvaluator gradEvaluator(&u, fct);
858 const_iterator iterEnd = u.template end<elem_t>();
859 for(const_iterator iter = u.template begin<elem_t>();
860 iter != iterEnd; ++iter)
861 {
862 // get the element
863 elem_t* elem = *iter;
864 vector_t elemGrad = gradEvaluator.evaluate(elem);
865
866 // iterate over all sides and add the normal component of the vector
867 // to the sides normGrad attachment
868 for(size_t i = 0; i < elem->num_sides(); ++i){
869 vector_t outerNorm = CalculateOuterNormal(elem, i, aaPos);
870 number ng = VecDot(elemGrad, outerNorm);
871 side_t* s = g.get_side(elem, i);
872 aaSideError[s] += ng;
873 ++aaNumElems[s];
874 }
875 }
876
877 ExchangeAndAdjustSideErrors<side_t>(u, aSideError, aNumElems);
878
879// now let's iterate over all elements again, this time integrating the jump
880// over the side and summing everything in aaError. Note that each element receives
881// 50% of a sides error
884 for(const_iterator iter = u.template begin<elem_t>();
885 iter != iterEnd; ++iter)
886 {
887 elem_t* elem = *iter;
888 number err = 0;
889 g.associated_elements(sides, elem);
890 for(size_t i = 0; i < sides.size(); ++i){
891 side_t* s = sides[i];
892 if(aaNumElems[s] > 1){
893 g.associated_elements(edges, s);
894 number a = CalculateVolume(s, aaPos);
895 number hs = ElementDiameter(g, aaPos, s);
896 err += hs * sq(a * aaSideError[s]);
897 }
898 }
899
900 if(addErrSquareToAAError)
901 aaError[elem] += 0.5 * err;
902 else
903 aaError[elem] = sqrt(0.5 * err);
904 }
905
906 g.template detach_from<side_t>(aSideError);
907 g.template detach_from<side_t>(aNumElems);
908}
909
910
913template <typename TGradientEvaluator, typename TFunction>
914void EvaluateGradientJump_Norm(TFunction& u, size_t fct,
916 typename TFunction::element_type,
917 ug::Attachment<number> >& aaError)
918{
919 using std::max;
920
921 static const int dim = TFunction::dim;
922 typedef typename TFunction::domain_type::grid_type grid_t;
923 typedef typename TFunction::const_element_iterator const_iterator;
924 typedef typename TFunction::element_type elem_t;
925 typedef typename elem_t::side side_t;
926 typedef MathVector<dim> vector_t;
927
928// get position accessor
929 typename TFunction::domain_type::position_accessor_type& aaPos
930 = u.domain()->position_accessor();
931
932// some storage
934
935// we need attachments on the sides to gather gradient-jumps in parallel.
936// Note that a serial version could be created without this additional attachment.
937 grid_t& g = *u.domain()->grid();
938 ANumber aSideError;
939 ANumber aNumElems;
940 g.template attach_to_dv<side_t>(aSideError, 0);
941 g.template attach_to_dv<side_t>(aNumElems, 0);
942 Grid::AttachmentAccessor<side_t, ANumber> aaSideError(g, aSideError);
943 Grid::AttachmentAccessor<side_t, ANumber> aaNumElems(g, aNumElems);
944
945// loop elements and evaluate gradient
946 TGradientEvaluator gradEvaluator(&u, fct);
947 const_iterator iterEnd = u.template end<elem_t>();
948 for(const_iterator iter = u.template begin<elem_t>();
949 iter != iterEnd; ++iter)
950 {
951 // get the element
952 elem_t* elem = *iter;
953 vector_t elemGrad = gradEvaluator.evaluate(elem);
954
955 number elemContrib = VecTwoNorm(elemGrad);
956 aaError[elem] = elemContrib;
957 g.associated_elements(sides, elem);
958 for(size_t i = 0; i < sides.size(); ++i){
959 side_t* s = sides[i];
960 aaSideError[s] += elemContrib;
961 ++aaNumElems[s];
962 }
963 }
964
965 ExchangeAndAdjustSideErrors<side_t>(u, aSideError, aNumElems);
966
967// finally store the highest difference between an element-error and
968// averaged errors in associated sides in each element-error.
969 for(const_iterator iter = u.template begin<elem_t>();
970 iter != iterEnd; ++iter)
971 {
972 elem_t* elem = *iter;
973 const number elemErr = aaError[elem];
974 number err = 0;
975 g.associated_elements(sides, elem);
976 for(size_t i = 0; i < sides.size(); ++i){
977 side_t* s = sides[i];
978 if(aaNumElems[s] > 0)
979 err = max(err, fabs(elemErr - aaSideError[s] / aaNumElems[s]));
980 }
981 //aaError[elem] = 2. * err * pow(CalculateVolume(elem, aaPos), number(dim-1)/number(dim));
982 aaError[elem] = 2. * err * pow(CalculateVolume(elem, aaPos), 2./number(dim));
983 //aaError[elem] = 2. * err * pow(CalculateVolume(elem, aaPos), 2./ (1. + 0.5*number(dim)));
984 //aaError[elem] = 2. * err * CalculateVolume(elem, aaPos);
985 }
986
987 g.template detach_from<side_t>(aSideError);
988 g.template detach_from<side_t>(aNumElems);
989}
990
1005template <typename TDomain, typename TAlgebra>
1008 const char* cmp,
1009 number refFrac,
1010 int minLvl, int maxLvl,
1011 std::string jumpType)
1012{
1013 PROFILE_FUNC();
1014 using namespace std;
1015// types
1016 typedef GridFunction<TDomain, TAlgebra> TFunction;
1017 typedef typename TFunction::domain_type::grid_type grid_t;
1018 typedef typename TFunction::element_type elem_t;
1019 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
1020
1021// function id
1022 const size_t fct = u->fct_id_by_name(cmp);
1023 UG_COND_THROW(fct >= u->num_fct(),
1024 "Function space does not contain a function with name " << cmp);
1025
1026// get multigrid and position accessor
1027 grid_t& mg = *u->domain()->grid();
1028
1029// attach error field
1031 ANumber aError;
1032 mg.template attach_to<elem_t>(aError);
1034
1035//todo: the type of the used gradient evaluator should depend on the used grid function.
1036 typedef GradientEvaluator_LagrangeP1<TFunction> LagrangeP1Evaluator;
1037 if(jumpType == std::string("norm"))
1038 EvaluateGradientJump_Norm<LagrangeP1Evaluator>(*u, fct, aaError);
1039 else if(jumpType == std::string("sideInt"))
1040 EvaluateGradientJump_SideIntegral<LagrangeP1Evaluator>(*u, fct, aaError);
1041 else{
1042 UG_THROW("Unsupported jumpType in MarkForAdaption_GradientJump: "
1043 "Valid values are: norm, sideInt");
1044 }
1045
1046
1047 number maxElemError = 0; // error in elements which may be refined
1048 number minElemError = numeric_limits<number>::max();
1049 for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
1050 if(mg.get_level(*iter) < maxLvl){
1051 maxElemError = max(maxElemError, aaError[*iter]);
1052 minElemError = min(minElemError, aaError[*iter]);
1053 }
1054 }
1055
1056 #ifdef UG_PARALLEL
1058 minElemError = com.allreduce(minElemError, PCL_RO_MIN);
1059 maxElemError = com.allreduce(maxElemError, PCL_RO_MAX);
1060 #endif
1061
1062// note that aaError contains error-squares
1063 number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1064 MarkElementsAbsolute(aaError, refiner, u->dof_distribution(), refTol, -1,
1065 minLvl, maxLvl);
1066
1067// detach error field
1068 mg.template detach_from<elem_t>(aError);
1069};
1070
1071
1073template <typename TFunction>
1076 const char* cmp,
1077 number time,
1078 int quadOrder, std::string quadType,
1080 typename TFunction::element_type,
1081 ug::Attachment<number> >& aaError)
1082{
1083 using namespace std;
1084// types
1085 typedef typename TFunction::domain_type::grid_type grid_t;
1086 typedef typename TFunction::element_type elem_t;
1087 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
1088 const int dim = TFunction::dim;
1089
1090// function id
1091 const size_t fct = u->fct_id_by_name(cmp);
1092 UG_COND_THROW(fct >= u->num_fct(),
1093 "Function space does not contain a function with name " << cmp);
1094
1095// get multigrid and position accessor
1096 grid_t& mg = *u->domain()->grid();
1097 typename TFunction::domain_type::position_accessor_type&
1098 aaPos = u->domain()->position_accessor();
1099
1100// evaluate L2-Norm of f and store element contributions in aaError
1101 /*SmartPtr<IIntegrand<number, dim> > spIntegrand
1102 = make_sp(new UserDataIntegrand<number, TFunction>(f, u, time));*/
1103 UserDataIntegrand<number, TFunction> integrand(f, &(*u), time);
1104 Integrate<dim, dim>(u->template begin<elem_t>(), u->template end<elem_t>(),
1105 aaPos, integrand, quadOrder, quadType, &aaError);
1106
1107// we have to square contributions in aaError and multiply them with the
1108// square of the element-diameter
1109 for(ElemIter iter = u->template begin<elem_t>();
1110 iter != u->template end<elem_t>(); ++iter)
1111 {
1112 elem_t* elem = *iter;
1113 aaError[elem] *= aaError[elem] * ElementDiameterSq(mg, aaPos, elem);
1114 }
1115
1116// now evaluate contributions of the integral over gradient jumps on the element sides
1117// their squares will be added to aaError.
1118 typedef GradientEvaluator_LagrangeP1<TFunction> LagrangeP1Evaluator;
1119 EvaluateGradientJump_SideIntegral<LagrangeP1Evaluator>(*u, fct, aaError, true);
1120
1121// finally take the square root of aaError
1122 for(ElemIter iter = u->template begin<elem_t>();
1123 iter != u->template end<elem_t>(); ++iter)
1124 {
1125 aaError[*iter] = sqrt(aaError[*iter]);
1126 }
1127};
1128
1143template <typename TDomain, typename TAlgebra>
1147 const char* cmp,
1148 number time,
1149 number refTol,
1150 number coarsenTol,
1151 int maxLvl,
1152 int quadOrder, std::string quadType,
1153 bool refTopLvlOnly = false)
1154{
1155 PROFILE_FUNC();
1156 using namespace std;
1157// types
1158 typedef GridFunction<TDomain, TAlgebra> TFunction;
1159 typedef typename TFunction::domain_type::grid_type grid_t;
1160 typedef typename TFunction::element_type elem_t;
1161 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
1162
1163// attach error field
1164 grid_t& mg = *u->domain()->grid();
1166 ANumber aError;
1167 mg.template attach_to<elem_t>(aError);
1169
1170// Evaluate the residual error
1171 EvaluateResidualErrorP1(u, f, cmp, time, quadOrder, quadType, aaError);
1172
1173// mark
1174 MarkElementsAbsolute(aaError, refiner, u->dof_distribution(), refTol, coarsenTol,
1175 0, maxLvl, refTopLvlOnly);
1176
1177// evaluate fraction of error in marked elements compared to the total error
1178 number errs[2] = {0, 0}; //0: marked error, 1: total error
1179
1180 for(ElemIter iter = u->template begin<elem_t>();
1181 iter != u->template end<elem_t>(); ++iter)
1182 {
1183 elem_t* e = *iter;
1184 errs[1] += sq(aaError[e]);
1185 if(refiner.get_mark(e) & (RM_REFINE | RM_ANISOTROPIC))
1186 errs[0] += sq(aaError[e]);
1187 }
1188
1189 #ifdef UG_PARALLEL
1191 number gErrs[2];
1192 com.allreduce(errs, gErrs, 2, PCL_RO_SUM);
1193 #else
1194 number* gErrs = errs;
1195 #endif
1196
1197 number frac = 1;
1198 if(gErrs[1] > 0)
1199 frac = sqrt(gErrs[0] / gErrs[1]);
1200
1201// detach error field
1202 mg.template detach_from<elem_t>(aError);
1203 return frac;
1204};
1205
1219template <typename TDomain, typename TAlgebra>
1223 const char* cmp,
1224 number time,
1225 number refFrac,
1226 int minLvl, int maxLvl,
1227 int quadOrder, std::string quadType)
1228{
1229 PROFILE_FUNC();
1230 using namespace std;
1231// types
1232 typedef GridFunction<TDomain, TAlgebra> TFunction;
1233 typedef typename TFunction::domain_type::grid_type grid_t;
1234 typedef typename TFunction::element_type elem_t;
1235 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
1236
1237// attach error field
1238 grid_t& mg = *u->domain()->grid();
1240 ANumber aError;
1241 mg.template attach_to<elem_t>(aError);
1243
1244// Evaluate the residual error
1245 EvaluateResidualErrorP1(u, f, cmp, time, quadOrder, quadType, aaError);
1246
1247// find min- and max-values
1248 number maxElemError = 0; // error in elements which may be refined
1249 number minElemError = numeric_limits<number>::max();
1250 for(ElemIter iter = u->template begin<elem_t>();
1251 iter != u->template end<elem_t>(); ++iter)
1252 {
1253 if(mg.get_level(*iter) < maxLvl){
1254 number err = aaError[*iter] = sqrt(aaError[*iter]);
1255 maxElemError = max(maxElemError, err);
1256 minElemError = min(minElemError, err);
1257 }
1258 }
1259
1260 #ifdef UG_PARALLEL
1262 minElemError = com.allreduce(minElemError, PCL_RO_MIN);
1263 maxElemError = com.allreduce(maxElemError, PCL_RO_MAX);
1264 #endif
1265
1266 number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1267 MarkElementsAbsolute(aaError, refiner, u->dof_distribution(), refTol, -1,
1268 minLvl, maxLvl);
1269
1270// detach error field
1271 mg.template detach_from<elem_t>(aError);
1272};
1273
1274
1275template <typename TDomain, typename TAlgebra>
1278 const char* cmp,
1279 number refFrac,
1280 int minLvl, int maxLvl)
1281{
1282 PROFILE_FUNC();
1283 using namespace std;
1284// types
1285 typedef GridFunction<TDomain, TAlgebra> TFunction;
1286 static const int dim = TFunction::dim;
1287 typedef typename TFunction::domain_type::grid_type grid_t;
1288 typedef typename TFunction::element_type elem_t;
1289 typedef typename TFunction::template traits<elem_t>::const_iterator ElemIter;
1290 typedef MathVector<dim> vector_t;
1291
1292// get position accessor
1293 typename TFunction::domain_type::position_accessor_type& aaPos
1294 = u->domain()->position_accessor();
1295
1296// function id
1297 const size_t fct = u->fct_id_by_name(cmp);
1298 UG_COND_THROW(fct >= u->num_fct(),
1299 "Function space does not contain a function with name " << cmp);
1300
1301// get multigrid and position accessor
1302 grid_t& mg = *u->domain()->grid();
1303
1304// attach error field
1306 ANumber aError;
1307 mg.template attach_to<elem_t>(aError);
1308 MultiGrid::AttachmentAccessor<elem_t, ANumber> aaErrorElem(mg, aError);
1309
1310// attach gradient to vertices and initialize it with zero
1311 typedef Attachment<vector_t> AGrad;
1312 AGrad aGrad;
1313 mg.attach_to_vertices(aGrad);
1315 vector_t zeroVec;
1316 VecSet(zeroVec, 0);
1317 SetAttachmentValues(aaGradVrt, mg.template begin<Vertex>(),
1318 mg.template end<Vertex>(), zeroVec);
1319
1320// attach contributor-counter to vertices. Required for parallel execution!
1321 ANumber aNumContribs;
1322 mg.attach_to_vertices(aNumContribs);
1323 MultiGrid::AttachmentAccessor<Vertex, ANumber> aaNumContribsVrt(mg, aNumContribs);
1324 SetAttachmentValues(aaNumContribsVrt, mg.template begin<Vertex>(),
1325 mg.template end<Vertex>(), 0);
1326
1327// average the gradients in the grids vertices
1328//todo: the type of the used gradient evaluator should depend on the used grid function.
1329 typedef GradientEvaluator_LagrangeP1<TFunction> LagrangeP1Evaluator;
1330
1331// loop elements and evaluate gradient
1332 LagrangeP1Evaluator gradEvaluator(u.get(), fct);
1334
1335 ElemIter iterElemEnd = u->template end<elem_t>();
1336 for(ElemIter iter = u->template begin<elem_t>();
1337 iter != iterElemEnd; ++iter)
1338 {
1339 // get the element
1340 elem_t* elem = *iter;
1341 vector_t elemGrad = gradEvaluator.evaluate(elem);
1342
1343 mg.associated_elements(vrts, elem);
1344 for(size_t i = 0; i < vrts.size(); ++i){
1345 Vertex* v = vrts[i];
1346 aaGradVrt[v] += elemGrad;
1347 ++aaNumContribsVrt[v];
1348 }
1349 }
1350
1351// in a parallel environment we now have to sum the gradients over parallel
1352// interfaces
1353 #ifdef UG_PARALLEL
1354 typedef typename GridLayoutMap::Types<Vertex>::Layout layout_t;
1355 DistributedGridManager& dgm = *mg.distributed_grid_manager();
1356 GridLayoutMap& glm = dgm.grid_layout_map();
1358
1359 // sum all copies at the h-master attachment
1360 ComPol_AttachmentReduce<layout_t, AGrad> compolSumGrad(mg, aGrad, PCL_RO_SUM);
1361 ComPol_AttachmentReduce<layout_t, ANumber> compolSumNum(mg, aNumContribs, PCL_RO_SUM);
1362 icom.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolSumGrad);
1363 icom.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolSumNum);
1364 icom.communicate();
1365
1366 // copy the sum from the master to all of its slave-copies
1367 ComPol_CopyAttachment<layout_t, AGrad> compolCopyGrad(mg, aGrad);
1368 ComPol_CopyAttachment<layout_t, ANumber> compolCopyNum(mg, aNumContribs);
1369 icom.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolCopyGrad);
1370 icom.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolCopyNum);
1371 icom.communicate();
1372
1373 //todo: communication to vmasters may not be necessary here...
1374 // it is performed to make sure that all surface-rim-children
1375 // contain their true value.
1376 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyGrad);
1377 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyNum);
1378
1379 icom.communicate();
1380 #endif
1381
1382// if we're currently operating on a surface function, we have to adjust
1383// errors between shadowed and shadowing vertices
1384 if(u->grid_level().type() == GridLevel::SURFACE){
1385 //todo: avoid iteration over the whole grid!
1386 //todo: reduce communication
1387 const SurfaceView* surfView = u->approx_space()->surface_view().get();
1388
1389 // for(side_iter_t iter = u.template begin<side_t>(SurfaceView::SHADOW_RIM);
1390 // iter != u.template end<side_t>(SurfaceView::SHADOW_RIM); ++iter)
1391 typedef Grid::traits<Vertex>::iterator grid_side_iter_t;
1392 for(grid_side_iter_t iter = mg.template begin<Vertex>();
1393 iter != mg.template end<Vertex>(); ++iter)
1394 {
1395 Vertex* s = *iter;
1397 continue;
1398
1399 Vertex* c = mg.get_child_vertex(s);
1400 if(c){
1401 aaGradVrt[s] += aaGradVrt[c];
1402 aaNumContribsVrt[s] += aaNumContribsVrt[c];
1403 }
1404 }
1405 // those elems which have local children now contain their true value (vslaves...)
1406 #ifdef UG_PARALLEL
1407 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyGrad);
1408 icom.exchange_data(glm, INT_V_SLAVE, INT_V_MASTER, compolCopyNum);
1409 icom.communicate();
1410 #endif
1411
1412 // copy values up to children
1413 for(grid_side_iter_t iter = mg.template begin<Vertex>();
1414 iter != mg.template end<Vertex>(); ++iter)
1415 {
1416 Vertex* s = *iter;
1418 continue;
1419
1420 Vertex* c = mg.get_child_vertex(s);
1421 if(c){
1422 aaGradVrt[c] = aaGradVrt[s];
1423 aaNumContribsVrt[c] = aaNumContribsVrt[s];
1424 }
1425 }
1426
1427 // we'll also average values in h-nodes now. If a parent is locally available,
1428 // it's associated values are correct at this point. Communication from vmaster
1429 // to vslaves is performed afterwards.
1430 typedef MultiGrid::traits<ConstrainedVertex>::iterator constr_vrt_iter;
1431 for(constr_vrt_iter iter = mg.template begin<ConstrainedVertex>();
1432 iter != mg.template end<ConstrainedVertex>(); ++iter)
1433 {
1434 ConstrainedVertex* v = *iter;
1436 if(p){
1437 VecSet(aaGradVrt[v], 0);
1438 aaNumContribsVrt[v] = 0;
1439 mg.associated_elements(vrts, p);
1440 int numConstr = 0;
1441 for(size_t i = 0; i < vrts.size(); ++i){
1443 continue;
1444 aaGradVrt[v] += aaGradVrt[vrts[i]];
1445 aaNumContribsVrt[v] += aaNumContribsVrt[vrts[i]];
1446 ++numConstr;
1447 }
1448 aaGradVrt[v] /= (number)numConstr;
1449 aaNumContribsVrt[v] /= (number)numConstr;
1450 }
1451 }
1452
1453 // those elems which have a local parent now contain their true value (vmasters...)
1454 #ifdef UG_PARALLEL
1455 // copy from v-masters to v-slaves, since there may be constrained
1456 // sides which locally do not have a constraining element. Note that
1457 // constrained V-Masters always have a local constraining element and thus
1458 // contain the correct value.
1459 icom.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, compolCopyGrad);
1460 icom.exchange_data(glm, INT_V_MASTER, INT_V_SLAVE, compolCopyNum);
1461 icom.communicate();
1462 #endif
1463 }
1464
1465// evaluate integral over difference of element gradients and averaged
1466// vertex gradients
1467 for(ElemIter iter = u->template begin<elem_t>();
1468 iter != iterElemEnd; ++iter)
1469 {
1470 // get the element
1471 elem_t* elem = *iter;
1472 vector_t elemGrad = gradEvaluator.evaluate(elem);
1473
1474 vector_t vrtAvrgGrad;
1475 VecSet(vrtAvrgGrad, 0);
1476 mg.associated_elements(vrts, elem);
1477 for(size_t i = 0; i < vrts.size(); ++i){
1478 Vertex* v = vrts[i];
1479 vector_t vg = aaGradVrt[v];
1480 vg /= aaNumContribsVrt[v];
1481 vrtAvrgGrad += vg;
1482 }
1483 vrtAvrgGrad /= (number)vrts.size();
1484 aaErrorElem[elem] = CalculateVolume(elem, aaPos)
1485 * VecDistance(elemGrad, vrtAvrgGrad)
1486 / (number)(dim+1);
1487 }
1488
1489// mark for adaption
1490 number maxElemError = 0; // error in elements which may be refined
1491 number minElemError = numeric_limits<number>::max();
1492 for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
1493 if(mg.get_level(*iter) < maxLvl){
1494 maxElemError = max(maxElemError, aaErrorElem[*iter]);
1495 minElemError = min(minElemError, aaErrorElem[*iter]);
1496 }
1497 }
1498
1499 #ifdef UG_PARALLEL
1501 minElemError = com.allreduce(minElemError, PCL_RO_MIN);
1502 maxElemError = com.allreduce(maxElemError, PCL_RO_MAX);
1503 #endif
1504
1505// note that aaErrorElem contains error-squares
1506 number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1507 MarkElementsAbsolute(aaErrorElem, refiner, u->dof_distribution(), refTol, -1,
1508 minLvl, maxLvl);
1509
1510// detach error field
1511 mg.detach_from_vertices(aNumContribs);
1512 mg.detach_from_vertices(aGrad);
1513 mg.template detach_from<elem_t>(aError);
1514};
1515
1516} // end namespace ug
1517
1518#endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__ERROR_INDICATOR__ */
parameterString p
parameterString s
Definition smart_pointer.h:108
Performs communication between interfaces on different processes.
Definition pcl_interface_communicator.h:68
bool communicate(int tag=749345)
sends and receives the collected data.
Definition pcl_interface_communicator_impl.hpp:409
void exchange_data(const TLayoutMap &layoutMap, const typename TLayoutMap::Key &keyFrom, const typename TLayoutMap::Key &keyTo, ICommunicationPolicy< TLayout > &commPol)
internally calls send_data and receive_data with the specified layouts.
Definition pcl_interface_communicator_impl.hpp:213
Definition pcl_process_communicator.h:70
void allreduce(const void *sendBuf, void *recBuf, int count, DataType type, ReduceOperation op) const
performs MPI_Allreduce on the processes of the communicator.
Definition pcl_process_communicator.cpp:318
A generic specialization of IAttachment.
Definition attachment_pipe.h:263
Performs reduce operations on the specified attachment.
Definition compol_attachment_reduce.h:274
copies values from a specified attachment to a stream and back.
Definition compol_copy_attachment.h:54
A vertex appearing on edges or faces.
Definition grid_objects_0d.h:110
GridObject * get_constraining_object()
Definition grid_objects_0d.h:150
dimension dependent base class for reference elements
Definition reference_element.h:183
virtual base class for reference mappings
Definition reference_mapping_provider.h:53
virtual number jacobian_transposed_inverse(MathMatrix< worldDim, dim > &JTInv, const MathVector< dim > &locPos) const =0
returns transposed of the inverse of the jacobian and returns sqrt of gram determinante
virtual void update(const MathVector< worldDim > *vCornerCoord)=0
refresh mapping for new set of corners
manages the layouts and interfaces which are associated with a distributed grid.
Definition distributed_grid.h:88
GridLayoutMap & grid_layout_map()
Definition distributed_grid.h:103
const LFEID & local_finite_element_id(size_t fct) const
Definition dof_distribution_info.h:225
size_t fct_id_by_name(const char *name) const
returns fct id by name
Definition dof_distribution_info.h:188
bool partially_contains(TStorageType flag) const
Definition flags.h:72
Definition gradient_evaluators.h:48
the generic attachment-accessor for access to grids attachment pipes.
Definition grid.h:182
represents numerical solutions on a grid using an algebraic vector
Definition grid_function.h:121
SmartPtr< TDomain > domain()
returns domain
Definition grid_function.h:342
SmartPtr< DoFDistribution > dof_distribution()
Definition grid_function.h:236
lets you access layouts by type and key
Definition parallel_grid_layout.h:152
@ SURFACE
Definition grid_level.h:48
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition grid_base_objects.h:157
The refiner interface allows to mark elements for refinement and to call refine.
Definition refiner_interface.h:67
virtual RefinementMark get_mark(Vertex *v) const
Returns the mark of a given element. Default returns RM_REFINE.
Definition refiner_interface.h:198
Definition integrate.h:1073
Identifier for Local Finite Elements.
Definition local_finite_element_id.h:98
@ CROUZEIX_RAVIART
Definition local_finite_element_id.h:105
@ PIECEWISE_CONSTANT
Definition local_finite_element_id.h:106
@ LAGRANGE
Definition local_finite_element_id.h:104
virtual size_t num_sh() const
Definition local_dof_set.cpp:46
virtual base class for local shape function sets
Definition local_shape_function_set.h:70
virtual void grads(grad_type *vGrad, const MathVector< dim > &x) const =0
returns all gradients evaluated at a point
A class for fixed size, dense matrices.
Definition math_matrix.h:63
a mathematical Vector with N entries.
Definition math_vector.h:97
Definition reference_element.h:345
size_t num(int dim) const
returns the number of geometric objects of dim
Definition reference_element.h:95
int id(int dim_i, size_t i, int dim_j, size_t j) const
id of object j in dimension dim_j of obj i in dimension dim_i
Definition reference_element.h:127
Definition reference_element.h:445
Definition reference_element.h:394
Definition reference_element.h:303
Represents the surface view of a multi-grid hierarchy.
Definition surface_view.h:61
@ SHADOW_RIM
Definition surface_view.h:98
@ SURFACE_RIM
Definition surface_view.h:94
SurfaceState surface_state(TElem *elem, const GridLevel &gl) const
returns the surface states, when considered as part of grid level
Type based UserData.
Definition user_data.h:143
For arbitrary UserData , this class defines the integrand .
Definition integrate.h:468
Base-class for all vertex-types.
Definition grid_base_objects.h:231
number ElementDiameter(const TElem &elem, TDomain &domain)
returns the maximal distance between to element vertices
number ElementDiameterSq(const TElem &elem, TDomain &domain)
returns the maximal squared distance between to element vertices
void CollectCornerCoordinates(std::vector< typename TAAPos::ValueType > &vCornerCoordsOut, const TElem &elem, const TAAPos &aaPos, bool clearContainer=true)
returns the corner coordinates of a geometric object
Definition domain_util_impl.h:75
void SetAttachmentValues(TAttachmentAccessor &aaVal, TIter elemsBegin, TIter elemsEnd, const TVal &val)
sets attachment-values for elements between elemsBegin and elemsEnd.
Definition attachment_util_impl.hpp:44
TAAPos::ValueType CalculateOuterNormal(Vertex *v, int sideIndex, TAAPos aaPos)
Calculates the outer normal of the i-th side of the given grid object.
Definition normal_calculation_impl.h:44
@ RM_REFINE
DEPRECATED. Use RM_FULL instead.
Definition refiner_interface.h:55
@ RM_ANISOTROPIC
DEPRECATED. Use RM_CLOSURE instead.
Definition refiner_interface.h:52
number CalculateVolume(Volume *elem, TAAPos aaPos)
Calculates the volume of the given element.
Definition volume_calculation_impl.hpp:43
#define PCL_RO_SUM
Definition pcl_methods.h:63
#define PCL_RO_MAX
Definition pcl_methods.h:61
#define PCL_RO_MIN
Definition pcl_methods.h:62
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
double number
Definition types.h:124
void MatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication.
Definition math_matrix_vector_functions_common_impl.hpp:49
TNumber sq(TNumber val)
returns the square of a value (val*val)
Definition math_util_impl.hpp:91
void VecSet(vector_t &vInOut, typename vector_t::value_type s)
Set each vector component to scalar (componentwise)
Definition math_vector_functions_common_impl.hpp:539
vector_t::value_type VecDistance(const vector_t &v1, const vector_t &v2)
returns the distance of two vector_ts.
Definition math_vector_functions_common_impl.hpp:375
vector_t::value_type VecTwoNorm(const vector_t &v)
Definition math_vector_functions_common_impl.hpp:577
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
returns the dot-product of two vector_ts
Definition math_vector_functions_common_impl.hpp:385
parameterNumber frac
Definition smart_pointer.h:814
the ug namespace
void EvaluateResidualErrorP1(SmartPtr< TFunction > u, SmartPtr< UserData< number, TFunction::dim > > f, const char *cmp, number time, int quadOrder, std::string quadType, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Evaluates the residual error for P1 shape functions (with some simplifications)
Definition error_indicator.h:1074
void MarkForAdaption_AbsoluteGradientJumpIndicator(IRefiner &refiner, GridFunction< TDomain, TAlgebra > &u, const char *fctName, number refTol, number coarsenTol, int minLvl, int maxLevel)
Definition error_indicator.h:547
void MarkForAdaption_AbsoluteGradientIndicator(IRefiner &refiner, GridFunction< TDomain, TAlgebra > &u, const char *fctName, number refTol, number coarsenTol, int minLvl, int maxLevel)
Definition error_indicator.h:413
void MarkForAdaption_GradientAverage(IRefiner &refiner, SmartPtr< GridFunction< TDomain, TAlgebra > > u, const char *cmp, number refFrac, int minLvl, int maxLvl)
Definition error_indicator.h:1276
ReferenceObjectID
these ids are used to identify the shape of a geometric object.
Definition grid_base_objects.h:74
void MarkForAdaption_GradientJump(IRefiner &refiner, SmartPtr< GridFunction< TDomain, TAlgebra > > u, const char *cmp, number refFrac, int minLvl, int maxLvl, std::string jumpType)
Definition error_indicator.h:1006
void MarkElementsAbsolute(MultiGrid::AttachmentAccessor< TElem, ug::Attachment< number > > &aaError, IRefiner &refiner, ConstSmartPtr< DoFDistribution > dd, number refTol, number coarsenTol, int minLevel, int maxLevel, bool refTopLvlOnly=false)
marks elements according to an attached error value field
Definition error_indicator_util.h:497
void EvaluateGradientJump_Norm(TFunction &u, size_t fct, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Definition error_indicator.h:914
void ComputeGradientCrouzeixRaviart(TFunction &u, size_t fct, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Definition error_indicator.h:149
void ExchangeAndAdjustSideErrors(TFunction &u, ANumber aSideError, ANumber aNumElems)
Definition error_indicator.h:714
void MarkForAdaption_ResidualErrorP1Relative(IRefiner &refiner, SmartPtr< GridFunction< TDomain, TAlgebra > > u, SmartPtr< UserData< number, TDomain::dim > > f, const char *cmp, number time, number refFrac, int minLvl, int maxLvl, int quadOrder, std::string quadType)
Definition error_indicator.h:1220
@ INT_V_MASTER
vertical master node
Definition parallel_grid_layout.h:106
@ INT_H_MASTER
horizontal master node
Definition parallel_grid_layout.h:104
@ INT_H_SLAVE
horizontal slave node
Definition parallel_grid_layout.h:105
@ INT_V_SLAVE
vertical slave node
Definition parallel_grid_layout.h:107
void MarkForAdaption_GradientJumpIndicator(IRefiner &refiner, GridFunction< TDomain, TAlgebra > &u, const char *fctName, number TOL, number refineFrac, number coarseFrac, int maxLevel)
Definition error_indicator.h:497
Attachment< number > ANumber
Definition common_attachments.h:57
void EvaluateGradientJump_SideIntegral(TFunction &u, size_t fct, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError, bool addErrSquareToAAError=false)
Definition error_indicator.h:827
void computeGradientJump(TFunction &u, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaGrad, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Definition error_indicator.h:460
void ComputeGradientPiecewiseConstant(TFunction &u, size_t fct, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Definition error_indicator.h:270
void MarkForAdaption_GradientIndicator(IRefiner &refiner, GridFunction< TDomain, TAlgebra > &u, const char *fctName, number TOL, number refineFrac, number coarseFrac, int maxLevel)
Definition error_indicator.h:367
number MarkForAdaption_ResidualErrorP1Absolute(IRefiner &refiner, SmartPtr< GridFunction< TDomain, TAlgebra > > u, SmartPtr< UserData< number, TDomain::dim > > f, const char *cmp, number time, number refTol, number coarsenTol, int maxLvl, int quadOrder, std::string quadType, bool refTopLvlOnly=false)
Definition error_indicator.h:1144
void MarkForAdaption_L2ErrorExact(IRefiner &refiner, SmartPtr< GridFunction< TDomain, TAlgebra > > u, SmartPtr< UserData< number, TDomain::dim > > spExactSol, const char *cmp, number minL2Error, number maxL2Error, number refFrac, int minLvl, int maxLvl, number time, int quadOrder)
Definition error_indicator.h:601
void ComputeGradientLagrange1(TFunction &u, size_t fct, MultiGrid::AttachmentAccessor< typename TFunction::element_type, ug::Attachment< number > > &aaError)
Definition error_indicator.h:59
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
#define PROFILE_FUNC()
Definition profiler.h:257
The traits class holds some important types for each element-type.
Definition grid.h:136
geometry_traits< TElem >::iterator iterator
Definition grid.h:143
PointerConstArray< TElem * > secure_container
Definition grid.h:146
pcl::MultiLevelLayout< Interface > Layout
Definition parallel_grid_layout.h:161
ReferenceVertex face_type1
Definition error_indicator.h:254
ReferenceVertex face_type0
Definition error_indicator.h:253
ReferenceEdge face_type0
Definition error_indicator.h:259
ReferenceEdge face_type1
Definition error_indicator.h:260
ReferenceTriangle face_type0
Definition error_indicator.h:265
ReferenceQuadrilateral face_type1
Definition error_indicator.h:266
Definition active_set_impl.h:45
void face_type0
Definition error_indicator.h:247
void face_type1
Definition error_indicator.h:248