Line data Source code
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 :
39 : #include "error_indicator_util.h"
40 : #include "gradient_evaluators.h"
41 : #include "common/common.h"
42 : #include "common/util/provider.h"
43 : #include "lib_grid/refinement/refiner_interface.h"
44 : #include "lib_disc/common/geometry_util.h"
45 : #include "lib_disc/reference_element/reference_element_util.h"
46 : #include "lib_disc/reference_element/reference_mapping_provider.h"
47 : #include "lib_disc/spatial_disc/disc_util/fvcr_geom.h"
48 : #include "integrate.h"
49 : #include "common/profiler/profiler.h"
50 :
51 : #ifdef UG_PARALLEL
52 : #include "lib_grid/parallelization/util/compol_attachment_reduce.h"
53 : #endif
54 :
55 : namespace ug{
56 :
57 :
58 : template <typename TFunction>
59 0 : void ComputeGradientLagrange1(TFunction& u, size_t fct,
60 : MultiGrid::AttachmentAccessor<
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 0 : = u.domain()->position_accessor();
71 :
72 : // some storage
73 : MathMatrix<dim, dim> JTInv;
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 0 : for(; iter != iterEnd; ++iter)
84 : {
85 : // get the element
86 : element_type* elem = *iter;
87 :
88 : // reference object type
89 0 : ReferenceObjectID roid = elem->reference_object_id();
90 :
91 : // get trial space
92 : const LocalShapeFunctionSet<dim>& lsfs =
93 0 : LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::LAGRANGE, dim, 1));
94 :
95 : // create a reference mapping
96 : DimReferenceMapping<dim, dim>& map
97 0 : = ReferenceMappingProvider::get<dim, dim>(roid);
98 :
99 : // get local Mid Point
100 0 : MathVector<dim> localIP = ReferenceElementCenter<dim>(roid);
101 :
102 : // number of shape functions
103 0 : const size_t numSH = lsfs.num_sh();
104 0 : vLocalGrad.resize(numSH);
105 0 : vGlobalGrad.resize(numSH);
106 :
107 : // evaluate reference gradient at local midpoint
108 0 : lsfs.grads(&vLocalGrad[0], localIP);
109 :
110 : // get corners of element
111 : CollectCornerCoordinates(vCorner, *elem, aaPos);
112 :
113 : // update mapping
114 0 : map.update(&vCorner[0]);
115 :
116 : // compute jacobian
117 0 : map.jacobian_transposed_inverse(JTInv, localIP);
118 :
119 : // compute size (volume) of element
120 0 : 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 0 : 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 0 : 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 0 : aaError[elem] = VecTwoNorm(MidGrad) * pow(elemSize, 2./dim);
145 : }
146 0 : }
147 :
148 : template <typename TFunction>
149 0 : void ComputeGradientCrouzeixRaviart(TFunction& u, size_t fct,
150 : MultiGrid::AttachmentAccessor<
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 0 : = u.domain()->position_accessor();
166 :
167 0 : grid_type& grid = *u.domain()->grid();
168 :
169 : // some storage
170 : MathMatrix<dim, dim> JTInv;
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 0 : for(; iter != iterEnd; ++iter)
181 : {
182 : // get the element
183 : element_type* elem = *iter;
184 :
185 : // get sides of element
186 0 : grid.associated_elements_sorted(sides, elem );
187 :
188 : // reference object type
189 0 : ReferenceObjectID roid = elem->reference_object_id();
190 :
191 : // get trial space
192 : const LocalShapeFunctionSet<dim>& lsfs =
193 0 : LocalFiniteElementProvider::get<dim>(roid, LFEID(LFEID::CROUZEIX_RAVIART, dim, 1));
194 :
195 : // create a reference mapping
196 : DimReferenceMapping<dim, dim>& map
197 0 : = ReferenceMappingProvider::get<dim, dim>(roid);
198 :
199 : // get local Mid Point
200 0 : MathVector<dim> localIP = ReferenceElementCenter<dim>(roid);
201 :
202 : // number of shape functions
203 0 : const size_t numSH = lsfs.num_sh();
204 0 : vLocalGrad.resize(numSH);
205 0 : vGlobalGrad.resize(numSH);
206 :
207 : // evaluate reference gradient at local midpoint
208 0 : lsfs.grads(&vLocalGrad[0], localIP);
209 :
210 : // get corners of element
211 : CollectCornerCoordinates(vCorner, *elem, aaPos);
212 :
213 : // update mapping
214 0 : map.update(&vCorner[0]);
215 :
216 : // compute jacobian
217 0 : map.jacobian_transposed_inverse(JTInv, localIP);
218 :
219 : // compute size (volume) of element
220 0 : 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 0 : 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 0 : aaError[elem] = VecTwoNorm(MidGrad) * pow(elemSize, 2./dim);
242 : }
243 0 : }
244 :
245 : template <int dim> struct face_type_traits
246 : {
247 : typedef void face_type0;
248 : typedef void face_type1;
249 : };
250 :
251 : template <> struct face_type_traits<1>
252 : {
253 : typedef ReferenceVertex face_type0;
254 : typedef ReferenceVertex face_type1;
255 : };
256 :
257 : template <> struct face_type_traits<2>
258 : {
259 : typedef ReferenceEdge face_type0;
260 : typedef ReferenceEdge face_type1;
261 : };
262 :
263 : template <> struct face_type_traits<3>
264 : {
265 : typedef ReferenceTriangle face_type0;
266 : typedef ReferenceQuadrilateral face_type1;
267 : };
268 :
269 : template <typename TFunction>
270 0 : void ComputeGradientPiecewiseConstant(TFunction& u, size_t fct,
271 : MultiGrid::AttachmentAccessor<
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 0 : = u.domain()->position_accessor();
290 :
291 0 : grid_type& grid = *u.domain()->grid();
292 :
293 : // some storage
294 : MathMatrix<dim, dim> JTInv;
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 0 : 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 0 : grid.associated_elements_sorted(sides, elem );
314 :
315 : // reference object type
316 0 : 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 0 : 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 0 : 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 0 : for (size_t i=0;i<numOfAsso;i++){
336 : std::vector<DoFIndex> ind;
337 : u.inner_dof_indices(assoElements[i], fct, ind);
338 0 : faceValue+=DoFRef(u, ind[0]);
339 : }
340 0 : faceValue/=(number)numOfAsso;
341 : MathVector<dim> normal;
342 : size_t numSideCo = rRefElem.num(dim-1,s,0);
343 :
344 0 : for (size_t i = 0; i < numSideCo; ++i)
345 0 : 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 0 : if ((int)numSideCo==dim)
350 0 : ElementNormal<face_type0,dim>(normal,&vSideCoPos[0]);
351 : else
352 0 : ElementNormal<face_type1,dim>(normal,&vSideCoPos[0]);
353 :
354 0 : for (int d=0;d<dim;d++){
355 0 : vGlobalGrad[d] += faceValue * normal[d];
356 : }
357 : }
358 : vGlobalGrad/=(number)elemSize;
359 :
360 : // write result in array storage
361 0 : aaError[elem] = VecTwoNorm(vGlobalGrad) * pow(elemSize, 2./dim);
362 : }
363 0 : }
364 :
365 :
366 : template <typename TDomain, typename TAlgebra>
367 0 : void MarkForAdaption_GradientIndicator(IRefiner& refiner,
368 : GridFunction<TDomain, TAlgebra>& u,
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 0 : SmartPtr<grid_type> pMG = u.domain()->grid();
386 :
387 : // attach error field
388 : typedef Attachment<number> ANumber;
389 : ANumber aError;
390 0 : pMG->template attach_to<element_type>(aError);
391 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaError(*pMG, aError);
392 :
393 : // Compute error on elements
394 : if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
395 0 : ComputeGradientLagrange1(u, fct, aaError);
396 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
397 0 : ComputeGradientCrouzeixRaviart(u, fct, aaError);
398 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
399 0 : ComputeGradientPiecewiseConstant(u,fct,aaError);
400 : else{
401 0 : UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
402 : }
403 :
404 : // Mark elements for refinement
405 0 : 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 0 : };
410 :
411 :
412 : template <typename TDomain, typename TAlgebra>
413 0 : void MarkForAdaption_AbsoluteGradientIndicator(IRefiner& refiner,
414 : GridFunction<TDomain, TAlgebra>& u,
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 0 : SmartPtr<grid_type> pMG = u.domain()->grid();
433 :
434 : // attach error field
435 : typedef Attachment<number> ANumber;
436 : ANumber aError;
437 0 : pMG->template attach_to<element_type>(aError);
438 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaError(*pMG, aError);
439 :
440 : // Compute error on elements
441 : if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
442 0 : ComputeGradientLagrange1(u, fct, aaError);
443 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
444 0 : ComputeGradientCrouzeixRaviart(u, fct, aaError);
445 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
446 0 : ComputeGradientPiecewiseConstant(u,fct,aaError);
447 : else{
448 0 : UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
449 : }
450 :
451 : // Mark elements for refinement
452 0 : 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 0 : };
457 :
458 :
459 : template <typename TFunction>
460 0 : void computeGradientJump(TFunction& u,
461 : MultiGrid::AttachmentAccessor<
462 : typename TFunction::element_type,
463 : ug::Attachment<number> >& aaGrad,
464 : MultiGrid::AttachmentAccessor<
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 0 : 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 0 : 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 0 : grid.associated_elements(neighElements,side);
486 0 : if (neighElements.size()!=2) continue;
487 0 : number localJump = std::abs(aaGrad[neighElements[0]]-aaGrad[neighElements[1]]);
488 0 : for (size_t i=0;i<2;i++)
489 0 : if (aaError[neighElements[i]]<localJump) aaError[neighElements[i]]=localJump;
490 : }
491 0 : }
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
496 : template <typename TDomain, typename TAlgebra>
497 0 : void MarkForAdaption_GradientJumpIndicator(IRefiner& refiner,
498 : GridFunction<TDomain, TAlgebra>& u,
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 0 : SmartPtr<grid_type> pMG = u.domain()->grid();
516 :
517 : // attach error field
518 : typedef Attachment<number> ANumber;
519 : ANumber aGrad;
520 0 : pMG->template attach_to<element_type>(aGrad);
521 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaGrad(*pMG, aGrad);
522 :
523 : ANumber aError;
524 0 : pMG->template attach_to<element_type>(aError);
525 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaError(*pMG, aError);
526 :
527 : // Compute error on elements
528 : if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
529 0 : ComputeGradientLagrange1(u, fct, aaGrad);
530 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
531 0 : ComputeGradientCrouzeixRaviart(u, fct, aaGrad);
532 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
533 0 : ComputeGradientPiecewiseConstant(u,fct,aaGrad);
534 : else{
535 0 : UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
536 : }
537 0 : computeGradientJump(u,aaGrad,aaError);
538 :
539 : // Mark elements for refinement
540 0 : 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 0 : };
545 :
546 : template <typename TDomain, typename TAlgebra>
547 0 : void MarkForAdaption_AbsoluteGradientJumpIndicator(IRefiner& refiner,
548 : GridFunction<TDomain, TAlgebra>& u,
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 0 : SmartPtr<grid_type> pMG = u.domain()->grid();
567 :
568 : // attach error field
569 : typedef Attachment<number> ANumber;
570 : ANumber aGrad;
571 0 : pMG->template attach_to<element_type>(aGrad);
572 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaGrad(*pMG, aGrad);
573 :
574 : ANumber aError;
575 0 : pMG->template attach_to<element_type>(aError);
576 : MultiGrid::AttachmentAccessor<element_type, ANumber> aaError(*pMG, aError);
577 :
578 : // Compute error on elements
579 : if (u.local_finite_element_id(fct) == LFEID(LFEID::LAGRANGE, dim, 1))
580 0 : ComputeGradientLagrange1(u, fct, aaGrad);
581 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::CROUZEIX_RAVIART, dim, 1))
582 0 : ComputeGradientCrouzeixRaviart(u, fct, aaGrad);
583 : else if (u.local_finite_element_id(fct) == LFEID(LFEID::PIECEWISE_CONSTANT, dim, 0))
584 0 : ComputeGradientPiecewiseConstant(u,fct,aaGrad);
585 : else{
586 0 : UG_THROW("Non-supported finite element type " << u.local_finite_element_id(fct) << "\n");
587 : }
588 0 : computeGradientJump(u,aaGrad,aaError);
589 :
590 : // Mark elements for refinement
591 0 : 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 0 : };
596 :
597 : /**
598 : * \param[in] maxL2Error errors below maxL2Error are considered fine.
599 : */
600 : template <typename TDomain, typename TAlgebra>
601 0 : void MarkForAdaption_L2ErrorExact(IRefiner& refiner,
602 : SmartPtr<GridFunction<TDomain, TAlgebra> > u,
603 : SmartPtr<UserData<number, TDomain::dim> > spExactSol,
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 0 : 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 0 : grid_t& mg = *u->domain()->grid();
626 :
627 : typename TFunction::domain_type::position_accessor_type&
628 0 : aaPos = u->domain()->position_accessor();
629 :
630 : // attach error field
631 : typedef Attachment<number> ANumber;
632 : ANumber aError;
633 0 : mg.template attach_to<elem_t>(aError);
634 : MultiGrid::AttachmentAccessor<elem_t, ANumber> aaError(mg, 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 0 : L2ErrorIntegrand<TFunction> integrand(spExactSol, *u, fct, time);
640 : number l2Error =
641 0 : Integrate<dim, dim>(u->template begin<elem_t>(), u->template end<elem_t>(),
642 : aaPos, integrand, quadOrder, "best", &aaError);
643 :
644 : #ifdef UG_PARALLEL
645 : pcl::ProcessCommunicator com;
646 : l2Error = com.allreduce(l2Error, PCL_RO_SUM);
647 : #endif
648 :
649 0 : l2Error = sqrt(l2Error);
650 :
651 : UG_LOG("maxError " << maxL2Error << ", l2Error " << l2Error << std::endl);
652 :
653 0 : 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 0 : number maxElemError = 0; // error in elements which may be refined
658 0 : number minElemError = numeric_limits<number>::max();
659 : number fixedError = 0; // error in elements which can't be refined any further
660 0 : for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
661 : ++numElemsTotal;
662 0 : if(mg.get_level(*iter) < maxLvl){
663 : ++numElemsActive;
664 0 : maxElemError = max(maxElemError, aaError[*iter]);
665 0 : minElemError = min(minElemError, aaError[*iter]);
666 : }
667 : else{
668 : fixedError += aaError[*iter];
669 : }
670 : }
671 :
672 : #ifdef UG_PARALLEL
673 : pcl::ProcessCommunicator com;
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 0 : 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 0 : 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 0 : };
701 :
702 :
703 : /** Exchages error-values and nbr-element-numbers between distributed sides and
704 : * adjusts those values in rim-shadows and rim-shadowing elements on the fly
705 : * (e.g. constrained and constraining elements).
706 : *
707 : * \param[in] u The grid-function on which error-indication is based
708 : * \param[in] aSideError ANumber attachment on side_t:
709 : * The total error accumulated in element-sides.
710 : * \param[in] aNumElems ANumber attachment on side_t:
711 : * The number of elements which have contributed to aSideError.
712 : */
713 : template <class side_t, class TFunction>
714 0 : void 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 0 : 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();
731 : pcl::InterfaceCommunicator<layout_t > icom;
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 0 : if(u.grid_level().type() == GridLevel::SURFACE){
761 : //todo: avoid iteration over the whole grid!
762 : //todo: reduce communication
763 0 : 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 0 : iter != g.template end<side_t>(); ++iter)
770 : {
771 : side_t* s = *iter;
772 0 : if(!surfView->surface_state(s).partially_contains(SurfaceView::SHADOW_RIM))
773 0 : continue;
774 :
775 : const size_t numChildren = g.template num_children<side_t>(s);
776 0 : if(numChildren > 0){
777 0 : number w = 1. / number(numChildren);
778 0 : 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 0 : aaSideError[s] += w * aaSideError[c];
781 0 : 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 0 : iter != g.template end<side_t>(); ++iter)
794 : {
795 : side_t* s = *iter;
796 0 : if(!surfView->surface_state(s).partially_contains(SurfaceView::SHADOW_RIM))
797 0 : continue;
798 :
799 : const size_t numChildren = g.template num_children<side_t>(s);
800 0 : if(numChildren > 0){
801 0 : number w = 1. / number(numChildren);
802 0 : 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 0 : aaSideError[c] = w * aaSideError[s];
805 0 : 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 0 : }
822 :
823 : /** Calculates the jump between normal components of element-gradients on element sides,
824 : * then calculates the integral over those jumps on each side and finally sums those
825 : * integrals into aaError.*/
826 : template <typename TGradientEvaluator, typename TFunction>
827 0 : void EvaluateGradientJump_SideIntegral(TFunction& u, size_t fct,
828 : MultiGrid::AttachmentAccessor<
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 0 : = 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 0 : grid_t& g = *u.domain()->grid();
849 : ANumber aSideError;
850 : ANumber aNumElems;
851 0 : g.template attach_to_dv<side_t>(aSideError, 0);
852 0 : 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 0 : TGradientEvaluator gradEvaluator(&u, fct);
858 : const_iterator iterEnd = u.template end<elem_t>();
859 0 : for(const_iterator iter = u.template begin<elem_t>();
860 0 : iter != iterEnd; ++iter)
861 : {
862 : // get the element
863 : elem_t* elem = *iter;
864 0 : 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 0 : for(size_t i = 0; i < elem->num_sides(); ++i){
869 0 : vector_t outerNorm = CalculateOuterNormal(elem, i, aaPos);
870 : number ng = VecDot(elemGrad, outerNorm);
871 0 : side_t* s = g.get_side(elem, i);
872 0 : aaSideError[s] += ng;
873 0 : ++aaNumElems[s];
874 : }
875 : }
876 :
877 0 : 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
882 : typename Grid::traits<side_t>::secure_container sides;
883 : Grid::edge_traits::secure_container edges;
884 0 : for(const_iterator iter = u.template begin<elem_t>();
885 0 : iter != iterEnd; ++iter)
886 : {
887 : elem_t* elem = *iter;
888 : number err = 0;
889 : g.associated_elements(sides, elem);
890 0 : for(size_t i = 0; i < sides.size(); ++i){
891 : side_t* s = sides[i];
892 0 : if(aaNumElems[s] > 1){
893 : g.associated_elements(edges, s);
894 0 : number a = CalculateVolume(s, aaPos);
895 : number hs = ElementDiameter(g, aaPos, s);
896 0 : err += hs * sq(a * aaSideError[s]);
897 : }
898 : }
899 :
900 0 : if(addErrSquareToAAError)
901 0 : aaError[elem] += 0.5 * err;
902 : else
903 0 : 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 0 : }
909 :
910 :
911 : /** calculates the length of the gradient in each element and then
912 : * fills aaError depending on the difference in length between neighbored elements.*/
913 : template <typename TGradientEvaluator, typename TFunction>
914 0 : void EvaluateGradientJump_Norm(TFunction& u, size_t fct,
915 : MultiGrid::AttachmentAccessor<
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 0 : = u.domain()->position_accessor();
931 :
932 : // some storage
933 : typename Grid::traits<side_t>::secure_container sides;
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 0 : grid_t& g = *u.domain()->grid();
938 : ANumber aSideError;
939 : ANumber aNumElems;
940 0 : g.template attach_to_dv<side_t>(aSideError, 0);
941 0 : 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 0 : TGradientEvaluator gradEvaluator(&u, fct);
947 : const_iterator iterEnd = u.template end<elem_t>();
948 0 : for(const_iterator iter = u.template begin<elem_t>();
949 0 : iter != iterEnd; ++iter)
950 : {
951 : // get the element
952 : elem_t* elem = *iter;
953 0 : vector_t elemGrad = gradEvaluator.evaluate(elem);
954 :
955 : number elemContrib = VecTwoNorm(elemGrad);
956 0 : aaError[elem] = elemContrib;
957 : g.associated_elements(sides, elem);
958 0 : for(size_t i = 0; i < sides.size(); ++i){
959 : side_t* s = sides[i];
960 0 : aaSideError[s] += elemContrib;
961 0 : ++aaNumElems[s];
962 : }
963 : }
964 :
965 0 : 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 0 : for(const_iterator iter = u.template begin<elem_t>();
970 0 : iter != iterEnd; ++iter)
971 : {
972 : elem_t* elem = *iter;
973 0 : const number elemErr = aaError[elem];
974 0 : number err = 0;
975 : g.associated_elements(sides, elem);
976 0 : for(size_t i = 0; i < sides.size(); ++i){
977 : side_t* s = sides[i];
978 0 : if(aaNumElems[s] > 0)
979 0 : 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 0 : 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 0 : }
990 :
991 : /** This gradient jump error indicator runs in parallel environments and
992 : * works with h-nodes.
993 : *
994 : * \param[in] refFrac given minElemError and maxElemError, all elements with
995 : * an error > minElemError+refFrac*(maxElemError-minElemError)
996 : * will be refined.
997 : * \param[in] jumpType defines the type of jump that shall be evaluated:
998 : * - norm: the difference between gradient norms on
999 : * neighboring elements is regarded.
1000 : * - sideInt: The integral over the normal component
1001 : * of the gradient on each side is regarded.
1002 : *
1003 : * \todo: Add coarsenFrac and apply coarsen-marks, too.
1004 : */
1005 : template <typename TDomain, typename TAlgebra>
1006 0 : void MarkForAdaption_GradientJump(IRefiner& refiner,
1007 : SmartPtr<GridFunction<TDomain, TAlgebra> > u,
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 0 : 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 0 : grid_t& mg = *u->domain()->grid();
1028 :
1029 : // attach error field
1030 : typedef Attachment<number> ANumber;
1031 : ANumber aError;
1032 0 : mg.template attach_to<elem_t>(aError);
1033 : MultiGrid::AttachmentAccessor<elem_t, ANumber> aaError(mg, 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 0 : if(jumpType == std::string("norm"))
1038 0 : EvaluateGradientJump_Norm<LagrangeP1Evaluator>(*u, fct, aaError);
1039 0 : else if(jumpType == std::string("sideInt"))
1040 0 : EvaluateGradientJump_SideIntegral<LagrangeP1Evaluator>(*u, fct, aaError);
1041 : else{
1042 0 : UG_THROW("Unsupported jumpType in MarkForAdaption_GradientJump: "
1043 : "Valid values are: norm, sideInt");
1044 : }
1045 :
1046 :
1047 0 : number maxElemError = 0; // error in elements which may be refined
1048 0 : number minElemError = numeric_limits<number>::max();
1049 0 : for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
1050 0 : if(mg.get_level(*iter) < maxLvl){
1051 0 : maxElemError = max(maxElemError, aaError[*iter]);
1052 0 : minElemError = min(minElemError, aaError[*iter]);
1053 : }
1054 : }
1055 :
1056 : #ifdef UG_PARALLEL
1057 : pcl::ProcessCommunicator com;
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 0 : number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1064 0 : 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 0 : };
1070 :
1071 :
1072 : /// Evaluates the residual error for P1 shape functions (with some simplifications)
1073 : template <typename TFunction>
1074 0 : void EvaluateResidualErrorP1(SmartPtr<TFunction> u,
1075 : SmartPtr<UserData<number, TFunction::dim> > f,
1076 : const char* cmp,
1077 : number time,
1078 : int quadOrder, std::string quadType,
1079 : MultiGrid::AttachmentAccessor<
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 0 : 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 0 : grid_t& mg = *u->domain()->grid();
1097 : typename TFunction::domain_type::position_accessor_type&
1098 0 : 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 0 : UserDataIntegrand<number, TFunction> integrand(f, &(*u), time);
1104 0 : 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 0 : for(ElemIter iter = u->template begin<elem_t>();
1110 0 : iter != u->template end<elem_t>(); ++iter)
1111 : {
1112 : elem_t* elem = *iter;
1113 0 : 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 0 : EvaluateGradientJump_SideIntegral<LagrangeP1Evaluator>(*u, fct, aaError, true);
1120 :
1121 : // finally take the square root of aaError
1122 0 : for(ElemIter iter = u->template begin<elem_t>();
1123 0 : iter != u->template end<elem_t>(); ++iter)
1124 : {
1125 0 : aaError[*iter] = sqrt(aaError[*iter]);
1126 : }
1127 0 : };
1128 :
1129 : /** A classical residual error for the poisson problem on linear shape functions
1130 : * works with h-nodes.
1131 : *
1132 : * Evaluates the element residual error sqrt{hT^2 * ||RT||^2 + 0.5*sum(hS*||RS||^2)}
1133 : * where the sum contains all sides S of an element T. RT denotes the residuum
1134 : * and RS the gradient-jump over sides of T.
1135 : *
1136 : * \param[in] refTol Threshold value. Only elements with a higher residual
1137 : * error than refTol are refined.
1138 : *
1139 : * \returns the fraction of the residual error in marked elements compared to
1140 : * the total residual error.
1141 : * \todo: Add coarsenFrac and apply coarsen-marks, too.
1142 : */
1143 : template <typename TDomain, typename TAlgebra>
1144 0 : number MarkForAdaption_ResidualErrorP1Absolute(IRefiner& refiner,
1145 : SmartPtr<GridFunction<TDomain, TAlgebra> > u,
1146 : SmartPtr<UserData<number, TDomain::dim> > f,
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 0 : grid_t& mg = *u->domain()->grid();
1165 : typedef Attachment<number> ANumber;
1166 : ANumber aError;
1167 0 : mg.template attach_to<elem_t>(aError);
1168 : MultiGrid::AttachmentAccessor<elem_t, ANumber> aaError(mg, aError);
1169 :
1170 : // Evaluate the residual error
1171 0 : EvaluateResidualErrorP1(u, f, cmp, time, quadOrder, quadType, aaError);
1172 :
1173 : // mark
1174 0 : 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 0 : for(ElemIter iter = u->template begin<elem_t>();
1181 0 : iter != u->template end<elem_t>(); ++iter)
1182 : {
1183 : elem_t* e = *iter;
1184 0 : errs[1] += sq(aaError[e]);
1185 0 : if(refiner.get_mark(e) & (RM_REFINE | RM_ANISOTROPIC))
1186 0 : errs[0] += sq(aaError[e]);
1187 : }
1188 :
1189 : #ifdef UG_PARALLEL
1190 : pcl::ProcessCommunicator com;
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 0 : if(gErrs[1] > 0)
1199 0 : frac = sqrt(gErrs[0] / gErrs[1]);
1200 :
1201 : // detach error field
1202 : mg.template detach_from<elem_t>(aError);
1203 0 : return frac;
1204 : };
1205 :
1206 : /** A classical residual error for the poisson problem on linear shape functions
1207 : * works with h-nodes.
1208 : *
1209 : * Evaluates the element residual error sqrt{hT^2 * ||RT||^2 + 0.5*sum(hS*||RS||^2)}
1210 : * where the sum contains all sides S of an element T. RT denotes the residuum
1211 : * and RS the gradient-jump over sides of T.
1212 : *
1213 : * \param[in] refFrac given minElemError and maxElemError, all elements with
1214 : * an error > minElemError+refFrac*(maxElemError-minElemError)
1215 : * will be refined.
1216 : *
1217 : * \todo: Add coarsenFrac and apply coarsen-marks, too.
1218 : */
1219 : template <typename TDomain, typename TAlgebra>
1220 0 : void MarkForAdaption_ResidualErrorP1Relative(IRefiner& refiner,
1221 : SmartPtr<GridFunction<TDomain, TAlgebra> > u,
1222 : SmartPtr<UserData<number, TDomain::dim> > f,
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 0 : grid_t& mg = *u->domain()->grid();
1239 : typedef Attachment<number> ANumber;
1240 : ANumber aError;
1241 0 : mg.template attach_to<elem_t>(aError);
1242 : MultiGrid::AttachmentAccessor<elem_t, ANumber> aaError(mg, aError);
1243 :
1244 : // Evaluate the residual error
1245 0 : EvaluateResidualErrorP1(u, f, cmp, time, quadOrder, quadType, aaError);
1246 :
1247 : // find min- and max-values
1248 0 : number maxElemError = 0; // error in elements which may be refined
1249 0 : number minElemError = numeric_limits<number>::max();
1250 0 : for(ElemIter iter = u->template begin<elem_t>();
1251 0 : iter != u->template end<elem_t>(); ++iter)
1252 : {
1253 0 : if(mg.get_level(*iter) < maxLvl){
1254 0 : number err = aaError[*iter] = sqrt(aaError[*iter]);
1255 0 : maxElemError = max(maxElemError, err);
1256 0 : minElemError = min(minElemError, err);
1257 : }
1258 : }
1259 :
1260 : #ifdef UG_PARALLEL
1261 : pcl::ProcessCommunicator com;
1262 : minElemError = com.allreduce(minElemError, PCL_RO_MIN);
1263 : maxElemError = com.allreduce(maxElemError, PCL_RO_MAX);
1264 : #endif
1265 :
1266 0 : number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1267 0 : 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 0 : };
1273 :
1274 :
1275 : template <typename TDomain, typename TAlgebra>
1276 0 : void MarkForAdaption_GradientAverage(IRefiner& refiner,
1277 : SmartPtr<GridFunction<TDomain, TAlgebra> > u,
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 0 : = u->domain()->position_accessor();
1295 :
1296 : // function id
1297 : const size_t fct = u->fct_id_by_name(cmp);
1298 0 : 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 0 : grid_t& mg = *u->domain()->grid();
1303 :
1304 : // attach error field
1305 : typedef Attachment<number> ANumber;
1306 : ANumber aError;
1307 0 : 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);
1314 : MultiGrid::AttachmentAccessor<Vertex, AGrad> aaGradVrt(mg, 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 0 : LagrangeP1Evaluator gradEvaluator(u.get(), fct);
1333 : Grid::vertex_traits::secure_container vrts;
1334 :
1335 : ElemIter iterElemEnd = u->template end<elem_t>();
1336 0 : for(ElemIter iter = u->template begin<elem_t>();
1337 0 : iter != iterElemEnd; ++iter)
1338 : {
1339 : // get the element
1340 : elem_t* elem = *iter;
1341 0 : vector_t elemGrad = gradEvaluator.evaluate(elem);
1342 :
1343 : mg.associated_elements(vrts, elem);
1344 0 : for(size_t i = 0; i < vrts.size(); ++i){
1345 : Vertex* v = vrts[i];
1346 : aaGradVrt[v] += elemGrad;
1347 0 : ++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();
1357 : pcl::InterfaceCommunicator<layout_t> icom;
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 0 : if(u->grid_level().type() == GridLevel::SURFACE){
1385 : //todo: avoid iteration over the whole grid!
1386 : //todo: reduce communication
1387 0 : 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 0 : iter != mg.template end<Vertex>(); ++iter)
1394 : {
1395 : Vertex* s = *iter;
1396 0 : if(!surfView->surface_state(s).partially_contains(SurfaceView::SHADOW_RIM))
1397 0 : continue;
1398 :
1399 : Vertex* c = mg.get_child_vertex(s);
1400 0 : if(c){
1401 : aaGradVrt[s] += aaGradVrt[c];
1402 0 : 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 0 : iter != mg.template end<Vertex>(); ++iter)
1415 : {
1416 : Vertex* s = *iter;
1417 0 : if(!surfView->surface_state(s).partially_contains(SurfaceView::SHADOW_RIM))
1418 0 : continue;
1419 :
1420 : Vertex* c = mg.get_child_vertex(s);
1421 0 : if(c){
1422 : aaGradVrt[c] = aaGradVrt[s];
1423 0 : 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 0 : iter != mg.template end<ConstrainedVertex>(); ++iter)
1433 : {
1434 : ConstrainedVertex* v = *iter;
1435 : GridObject* p = v->get_constraining_object();
1436 0 : if(p){
1437 : VecSet(aaGradVrt[v], 0);
1438 0 : aaNumContribsVrt[v] = 0;
1439 : mg.associated_elements(vrts, p);
1440 : int numConstr = 0;
1441 0 : for(size_t i = 0; i < vrts.size(); ++i){
1442 0 : if(!surfView->surface_state(vrts[i]).partially_contains(SurfaceView::SURFACE_RIM))
1443 0 : continue;
1444 : aaGradVrt[v] += aaGradVrt[vrts[i]];
1445 0 : aaNumContribsVrt[v] += aaNumContribsVrt[vrts[i]];
1446 0 : ++numConstr;
1447 : }
1448 0 : aaGradVrt[v] /= (number)numConstr;
1449 0 : 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 0 : for(ElemIter iter = u->template begin<elem_t>();
1468 0 : iter != iterElemEnd; ++iter)
1469 : {
1470 : // get the element
1471 : elem_t* elem = *iter;
1472 0 : vector_t elemGrad = gradEvaluator.evaluate(elem);
1473 :
1474 : vector_t vrtAvrgGrad;
1475 : VecSet(vrtAvrgGrad, 0);
1476 : mg.associated_elements(vrts, elem);
1477 0 : 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 0 : vrtAvrgGrad /= (number)vrts.size();
1484 0 : aaErrorElem[elem] = CalculateVolume(elem, aaPos)
1485 0 : * VecDistance(elemGrad, vrtAvrgGrad)
1486 0 : / (number)(dim+1);
1487 : }
1488 :
1489 : // mark for adaption
1490 0 : number maxElemError = 0; // error in elements which may be refined
1491 0 : number minElemError = numeric_limits<number>::max();
1492 0 : for(ElemIter iter = u->template begin<elem_t>(); iter != u->template end<elem_t>(); ++iter){
1493 0 : if(mg.get_level(*iter) < maxLvl){
1494 0 : maxElemError = max(maxElemError, aaErrorElem[*iter]);
1495 0 : minElemError = min(minElemError, aaErrorElem[*iter]);
1496 : }
1497 : }
1498 :
1499 : #ifdef UG_PARALLEL
1500 : pcl::ProcessCommunicator com;
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 0 : number refTol = minElemError + (maxElemError - minElemError) * refFrac;
1507 0 : 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 0 : };
1515 :
1516 : } // end namespace ug
1517 :
1518 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__ERROR_INDICATOR__ */
|
