Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
AdaptiveMeshRefinement_def.hpp
1#ifndef AdaptiveMeshRefinement_def_hpp
2#define AdaptiveMeshRefinement_def_hpp
3
4#ifndef MESH_TIMER_START
5#define MESH_TIMER_START(A,S) Teuchos::RCP<Teuchos::TimeMonitor> A = Teuchos::rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(std::string("Mesh Refinement") + std::string(S))));
6#endif
7
8#ifndef MESH_TIMER_STOP
9#define MESH_TIMER_STOP(A) A.reset();
10#endif
11
12#include <chrono>
13#include <iomanip>
21
22using Teuchos::reduceAll;
23using Teuchos::REDUCE_SUM;
24using Teuchos::REDUCE_MAX;
25using Teuchos::REDUCE_MIN;
26using Teuchos::outArg;
27
28namespace FEDD {
29template <class SC, class LO, class GO, class NO>
30AdaptiveMeshRefinement<SC,LO,GO,NO>::AdaptiveMeshRefinement():
31inputMeshP1_(),
32inputMeshP12_(),
33outputMesh_(),
34errorElementsMv_(),
35errorEstimationMv_(0),
36domainsP1_(0)
37{
38
39}
40
41
42void dummyFuncSol(double* x, double* res){
43
44 res[0] = 0.;
45
46 return;
47}
54template <class SC, class LO, class GO, class NO>
55AdaptiveMeshRefinement<SC,LO,GO,NO>::AdaptiveMeshRefinement(ParameterListPtr_Type parameterListAll):
56inputMeshP1_(),
57inputMeshP12_(),
58outputMesh_(),
59errorElementsMv_(),
60errorEstimationMv_(0),
61domainsP1_(0)
62{
63 parameterListAll_ = parameterListAll;
64 this->dim_ = parameterListAll->sublist("Parameter").get("Dimension",2);;
65 hasProblemType_ = false;
66
67 FEType1_ = "P1";
68 FEType2_ = parameterListAll->sublist("Parameter").get("Discretization","P1");
69
70 exportWithParaview_ = false;
71
72 tol_= parameterListAll->sublist("Mesh Refinement").get("Toleranz",0.001);
73 theta_ = parameterListAll->sublist("Mesh Refinement").get("Theta",0.35);
74 markingStrategy_ = parameterListAll->sublist("Mesh Refinement").get("RefinementType","Uniform");
75 maxIter_ = parameterListAll->sublist("Mesh Refinement").get("MaxIter",3);
76
77 writeRefinementInfo_ = parameterListAll->sublist("Mesh Refinement").get("Write Refinement Info",true);
78 writeMeshQuality_ = parameterListAll->sublist("Mesh Refinement").get("Write Mesh Quality",true);
79
80 coarseningCycle_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening Cycle",0);
81 coarseningM_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening m",1);
82 coarseningN_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening n" ,1);
83
84 refinementMode_ = parameterListAll->sublist("Mesh Refinement").get("Refinement Mode" ,"Regular");
85
86 exactSolFunc_ = dummyFuncSol;
87 exactSolPFunc_ = dummyFuncSol;
88
89
90}
91
98template <class SC, class LO, class GO, class NO>
99AdaptiveMeshRefinement<SC,LO,GO,NO>::AdaptiveMeshRefinement(std::string problemType, ParameterListPtr_Type parameterListAll ):
100inputMeshP1_(),
101inputMeshP12_(),
102outputMesh_(),
103errorElementsMv_(),
104errorEstimationMv_(0),
105domainsP1_(0)
106{
107 parameterListAll_ = parameterListAll;
108 this->dim_ = parameterListAll->sublist("Parameter").get("Dimension",2);;
109 problemType_ = problemType;
110
111 FEType1_ = "P1";
112 FEType2_ = parameterListAll->sublist("Parameter").get("Discretization","P1");
113
114 exportWithParaview_ = false;
115
116 tol_= parameterListAll->sublist("Mesh Refinement").get("Toleranz",0.001);
117 theta_ = parameterListAll->sublist("Mesh Refinement").get("Theta",0.35);
118 markingStrategy_ = parameterListAll->sublist("Mesh Refinement").get("RefinementType","Uniform");
119 maxIter_ = parameterListAll->sublist("Mesh Refinement").get("MaxIter",3);
120
121 writeRefinementInfo_ = parameterListAll->sublist("Mesh Refinement").get("Write Refinement Info",true);
122 writeMeshQuality_ = parameterListAll->sublist("Mesh Refinement").get("Write Mesh Quality",true);
123
124 coarseningCycle_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening Cycle",0);
125 coarseningM_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening m",1);
126 coarseningN_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening n" ,1);
127
128 refinementMode_ = parameterListAll->sublist("Mesh Refinement").get("Refinement Mode" ,"Regular");
129
130 exactSolFunc_ = dummyFuncSol;
131 exactSolPFunc_ = dummyFuncSol;
132
133
134}
135
142template <class SC, class LO, class GO, class NO>
143AdaptiveMeshRefinement<SC,LO,GO,NO>::AdaptiveMeshRefinement(std::string problemType, ParameterListPtr_Type parameterListAll , Func_Type exactSolFunc ):
144inputMeshP1_(),
145inputMeshP12_(),
146outputMesh_(),
147errorElementsMv_(),
148errorEstimationMv_(0),
149domainsP1_(0)
150{
151 parameterListAll_ = parameterListAll;
152 this->dim_ = parameterListAll->sublist("Parameter").get("Dimension",2);;
153 this->problemType_ = problemType;
154
155 this->FEType1_ = "P1";
156 this->FEType2_ = parameterListAll->sublist("Parameter").get("Discretization","P1");
157
158 exactSolFunc_ = exactSolFunc;
159 exactSolInput_ = true;
160
161 exactSolPInput_ = false;
162
163 tol_= parameterListAll->sublist("Mesh Refinement").get("Toleranz",0.001);
164 theta_ = parameterListAll->sublist("Mesh Refinement").get("Theta",0.35);
165 markingStrategy_ = parameterListAll->sublist("Mesh Refinement").get("RefinementType","Uniform");
166 maxIter_ = parameterListAll->sublist("Mesh Refinement").get("MaxIter",3);
167
168 writeRefinementInfo_ = parameterListAll->sublist("Mesh Refinement").get("Write Refinement Info",true);
169 writeMeshQuality_ = parameterListAll->sublist("Mesh Refinement").get("Write Mesh Quality",true);
170
171 coarseningCycle_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening Cycle",0);
172 coarseningM_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening m",1);
173 coarseningN_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening n" ,1);
174
175 refinementMode_ = parameterListAll->sublist("Mesh Refinement").get("Refinement Mode" ,"Regular");
176
177 exactSolPFunc_ = dummyFuncSol;
178}
179
188template <class SC, class LO, class GO, class NO>
189AdaptiveMeshRefinement<SC,LO,GO,NO>::AdaptiveMeshRefinement(std::string problemType, ParameterListPtr_Type parameterListAll , Func_Type exactSolFuncU ,Func_Type exactSolFuncP ):
190inputMeshP1_(),
191inputMeshP12_(),
192outputMesh_(),
193errorElementsMv_(),
194errorEstimationMv_(0),
195domainsP1_(0)
196{
197 parameterListAll_ = parameterListAll;
198 this->dim_ = parameterListAll->sublist("Parameter").get("Dimension",2);;
199 this->problemType_ = problemType;
200
201 this->FEType1_ = "P1";
202 this->FEType2_ = parameterListAll->sublist("Parameter").get("Discretization","P1");
203
204 exactSolFunc_ = exactSolFuncU;
205 exactSolPFunc_ = exactSolFuncP;
206
207 exactSolInput_ = true;
208 exactSolPInput_ = true;
209 calculatePressure_ = true;
210
211 tol_= parameterListAll->sublist("Mesh Refinement").get("Toleranz",0.001);
212 theta_ = parameterListAll->sublist("Mesh Refinement").get("Theta",0.35);
213 markingStrategy_ = parameterListAll->sublist("Mesh Refinement").get("RefinementType","Uniform");
214 maxIter_ = parameterListAll->sublist("Mesh Refinement").get("MaxIter",3);
215
216 writeRefinementInfo_ = parameterListAll->sublist("Mesh Refinement").get("Write Refinement Info",true);
217 writeMeshQuality_ = parameterListAll->sublist("Mesh Refinement").get("Write Mesh Quality",true);
218
219 coarseningCycle_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening Cycle",0);
220 coarseningM_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening m",1);
221 coarseningN_ = parameterListAll->sublist("Mesh Refinement").get("Coarsening n" ,1);
222
223 refinementMode_ = parameterListAll->sublist("Mesh Refinement").get("Refinement Mode" ,"Regular");
224
225
226}
227template <class SC, class LO, class GO, class NO>
228AdaptiveMeshRefinement<SC,LO,GO,NO>::~AdaptiveMeshRefinement(){
229
230}
239template <class SC, class LO, class GO, class NO>
240typename AdaptiveMeshRefinement<SC,LO,GO,NO>::DomainPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>:: refineArea(DomainPtr_Type domainP1, vec2D_dbl_Type area, int level ){
241
242 DomainPtr_Type domainRefined(new Domain<SC,LO,GO,NO>( domainP1->getComm() , dim_ ));
243
244 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
245 inputMeshP1_->FEType_ = domainP1->getFEType();
246
247 inputMeshP1_->assignEdgeFlags();
248
249
250 MeshUnstrPtr_Type outputMesh(new MeshUnstr_Type(domainP1->getComm(), inputMeshP1_->volumeID_));
251
252 domainRefined->initWithDomain(domainP1);
253
254 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
255 inputMeshP1_->FEType_ = domainP1->getFEType();
256
257 // Error Estimation object
258 ErrorEstimation<SC,LO,GO,NO> errorEstimator (dim_, problemType_ , writeMeshQuality_);
259
260 // Refinement Factory object
261 RefinementFactory<SC,LO,GO,NO> refinementFactory( domainP1->getComm(), inputMeshP1_->volumeID_, parameterListAll_); // refinementRestriction_, refinement3DDiagonal_, restrictionLayer_);
262
263 // Estimating the error with the Discretizations Mesh.
264 int currentLevel =0;
265 while(currentLevel < level){
266 errorEstimator.tagArea(inputMeshP1_,area);
267 refinementFactory.refineMesh(inputMeshP1_,currentLevel, outputMesh, refinementMode_);
268
269 inputMeshP1_ = outputMesh;
270 currentLevel++;
271 }
272
273 domainRefined->setMesh(outputMesh);
274
275 return domainRefined;
276
277}
278
285template <class SC, class LO, class GO, class NO>
286typename AdaptiveMeshRefinement<SC,LO,GO,NO>::DomainPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>:: refineUniform(DomainPtr_Type domainP1, int level ){
287
288 DomainPtr_Type domainRefined(new Domain<SC,LO,GO,NO>( domainP1->getComm() , dim_ ));
289
290 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
291 inputMeshP1_->FEType_ = domainP1->getFEType();
292
293 inputMeshP1_->assignEdgeFlags();
294
295 MeshUnstrPtr_Type outputMesh(new MeshUnstr_Type(domainP1->getComm(), inputMeshP1_->volumeID_));
296
297 domainRefined->initWithDomain(domainP1);
298
299 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
300 inputMeshP1_->FEType_ = domainP1->getFEType();
301
302 // Error Estimation object
303 ErrorEstimation<SC,LO,GO,NO> errorEstimator (dim_, problemType_ , writeMeshQuality_);
304
305 // Refinement Factory object
306 RefinementFactory<SC,LO,GO,NO> refinementFactory( domainP1->getComm(), inputMeshP1_->volumeID_, parameterListAll_); // refinementRestriction_, refinement3DDiagonal_, restrictionLayer_);
307
308 // Estimating the error with the Discretizations Mesh.
309 int currentLevel =0;
310 while(currentLevel < level){
311 errorEstimator.tagAll(inputMeshP1_);
312 refinementFactory.refineMesh(inputMeshP1_,currentLevel, outputMesh, refinementMode_);
313
314 inputMeshP1_ = outputMesh;
315 currentLevel++;
316 }
317
318 domainRefined->setMesh(outputMesh);
319
320 return domainRefined;
321
322}
323
324template <class SC, class LO, class GO, class NO>
325typename AdaptiveMeshRefinement<SC,LO,GO,NO>::DomainPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>:: refineFlag(DomainPtr_Type domainP1, int level ,int flag){
326
327 DomainPtr_Type domainRefined(new Domain<SC,LO,GO,NO>( domainP1->getComm() , dim_ ));
328
329 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
330 inputMeshP1_->FEType_ = domainP1->getFEType();
331
332 inputMeshP1_->assignEdgeFlags();
333
334 MeshUnstrPtr_Type outputMesh(new MeshUnstr_Type(domainP1->getComm(), inputMeshP1_->volumeID_));
335
336 domainRefined->initWithDomain(domainP1);
337
338 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
339 inputMeshP1_->FEType_ = domainP1->getFEType();
340
341 // Error Estimation object
342 ErrorEstimation<SC,LO,GO,NO> errorEstimator (dim_, problemType_ , writeMeshQuality_);
343
344 // Refinement Factory object
345 RefinementFactory<SC,LO,GO,NO> refinementFactory( domainP1->getComm(), inputMeshP1_->volumeID_, parameterListAll_); // refinementRestriction_, refinement3DDiagonal_, restrictionLayer_);
346
347 // Estimating the error with the Discretizations Mesh.
348 int currentLevel =0;
349 while(currentLevel < level){
350 errorEstimator.tagFlag(inputMeshP1_,flag);
351 refinementFactory.refineMesh(inputMeshP1_,currentLevel, outputMesh, refinementMode_);
352
353 inputMeshP1_ = outputMesh;
354 currentLevel++;
355 }
356
357 domainRefined->setMesh(outputMesh);
358
359 return domainRefined;
360
361}
362
369template <class SC, class LO, class GO, class NO>
370void AdaptiveMeshRefinement<SC,LO,GO,NO>::identifyProblem(BlockMultiVectorConstPtr_Type valuesSolution){
371
372 // dofs can be determined by checking the solutions' blocks an comparing the length of the vectors to the number of unique nodes.
373 // If the length is the same: dofs =1 , if it's twice as long: dofs =2 .. and so on
374
375 // P_12 generally represents the velocity domain:
376 if(inputMeshP12_->getMapUnique()->getNodeNumElements() == valuesSolution->getBlock(0)->getDataNonConst(0).size())
377 dofs_ = 1;
378 else if(2*(inputMeshP12_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(0)->getDataNonConst(0).size())
379 dofs_ = 2;
380 else if(3*(inputMeshP12_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(0)->getDataNonConst(0).size())
381 dofs_ = 3;
382
383 if(valuesSolution->size() > 1){
384 if(inputMeshP1_->getMapUnique()->getNodeNumElements() == valuesSolution->getBlock(1)->getDataNonConst(0).size())
385 dofsP_ = 1;
386 else if(2*(inputMeshP1_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(1)->getDataNonConst(0).size())
387 dofsP_ = 2;
388 else if(3*(inputMeshP1_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(1)->getDataNonConst(0).size())
389 dofsP_ = 3;
390 calculatePressure_=true;
391 }
392
393}
394
406
407template <class SC, class LO, class GO, class NO>
408typename AdaptiveMeshRefinement<SC,LO,GO,NO>::DomainPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>::globalAlgorithm(DomainPtr_Type domainP1, DomainPtr_Type domainP12, BlockMultiVectorConstPtr_Type solution,ProblemPtr_Type problem, RhsFunc_Type rhsFunc ){
409 TEUCHOS_TEST_FOR_EXCEPTION( !hasProblemType_ , std::runtime_error, "No consideration of Problem Type. Please specify or only use: refineArea or refineUniform.");
410
411 solution_ = solution;
412
413 currentIter_ = domainsP1_.size() ;
414
415 rhsFunc_ = rhsFunc;
416
417 problem_ = problem;
418
419 comm_ = domainP1 ->getComm();
420
421 if(this->comm_->getRank() == 0 && currentIter_ < maxIter_){
422 std::cout << " -- Adaptive Mesh Refinement --" << std::endl;
423 std::cout << " " << std::endl;
424 }
425
426 maxRank_ = std::get<1>(domainP1->getMesh()->rankRange_);
427 // We save the domains of each step
428 // The P1 Mesh is always used for refinement while the P1 or P2 Mesh is used for error Estimation depending on Discretisation
429 domainsP1_.push_back(domainP1);
430 domainsP12_.push_back(domainP12);
431
432 domainP1_ = domainP1;
433 domainP12_ = domainP12;
434
435
436 // Reading Mesh from domainP1 as we always refine the P1 Mesh, here defined as inputMesh_
437 inputMeshP1_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP1->getMesh() , true);
438 inputMeshP1_->FEType_ = domainP1->getFEType();
439
440 // With the global Algorithm we create a new P1 domain with a new mesh
441 DomainPtr_Type domainRefined(new Domain<SC,LO,GO,NO>( domainP1->getComm() , dim_ ));
442
443 // Output Mesh
444 MeshUnstrPtr_Type outputMesh(new MeshUnstr_Type(domainP1->getComm(), inputMeshP1_->volumeID_));
445
446 // Init to be refined domain with inputDomain
447 domainRefined->initWithDomain(domainP1);
448
449 inputMeshP12_ = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainP12->getMesh() , true);
450 inputMeshP12_->FEType_ = domainP12->getFEType();
451
452 this->identifyProblem(solution);
453
454 // Error Estimation object
455 ErrorEstimation<SC,LO,GO,NO> errorEstimator (dim_, problemType_, writeMeshQuality_ );
456
457 // Refinement Factory object
458 RefinementFactory<SC,LO,GO,NO> refinementFactory( domainP1->getComm(), inputMeshP1_->volumeID_, parameterListAll_);
459
460 if(dim_ ==3){
461 SurfaceElementsPtr_Type surfaceTriangleElements = inputMeshP12_->getSurfaceTriangleElements(); // Surfaces
462 if(surfaceTriangleElements.is_null()){
463 surfaceTriangleElements.reset(new SurfaceElements()); // Surface
464 refinementFactory.buildSurfaceTriangleElements( inputMeshP12_->getElementsC(),inputMeshP12_->getEdgeElements(),surfaceTriangleElements, inputMeshP12_->getEdgeMap(),inputMeshP12_->getElementMap() );
465 inputMeshP12_->surfaceTriangleElements_ = surfaceTriangleElements;
466 inputMeshP1_->surfaceTriangleElements_ = surfaceTriangleElements;
467 }
468 else if(surfaceTriangleElements->numberElements() ==0){
469 refinementFactory.buildSurfaceTriangleElements( inputMeshP12_->getElementsC(),inputMeshP12_->getEdgeElements(),inputMeshP12_->getSurfaceTriangleElements() , inputMeshP12_->getEdgeMap(),inputMeshP12_->getElementMap() );
470 inputMeshP12_->surfaceTriangleElements_ = surfaceTriangleElements;
471 inputMeshP1_->surfaceTriangleElements_ = surfaceTriangleElements;
472 }
473 }
474
475 if(currentIter_ == 0 && dim_ == 2){
476 inputMeshP1_->assignEdgeFlags();
477 inputMeshP12_->assignEdgeFlags();
478 }
479 // If coarsen Mesh is false, so consequently we refine the Mesh. we go about as folows:
480 bool coarsening= false;
481 if(coarseningCycle_ > 0 && currentIter_>0){
482 if(currentIter_ % coarseningCycle_ == 0)
483 coarsening = true;
484 }
485 errorElementsMv_ =Teuchos::rcp( new MultiVector_Type( domainP12_ ->getElementMap(), 1 ) );
486 errorElementsMv_->putScalar(0.);
487 if( coarsening== true && currentIter_ < maxIter_ ){
488
489 // We start by calculating the error of the current mesh. As this is out starting point for mesh coarsening. In the previous iteration we calculated the error estimation beforehand.
490 errorElementsMv_ = errorEstimator.estimateError(inputMeshP12_, inputMeshP1_, solution, rhsFunc_, domainP12->getFEType());
491
492 errorEstimationMv_.push_back(errorElementsMv_);
493
494 int m= coarseningM_;
495 int n= coarseningN_;
496 int k = currentIter_;
497 int iterC;
498 MeshUnstrPtrArray_Type meshUnstructuredP1(currentIter_+n);
499
500 for(int i=0; i<currentIter_+1-m; i++)
501 meshUnstructuredP1[i] = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainsP1_[i]->getMesh() , true);
502
503 // We extract the error estimation of the mesh iter
504
505 MultiVectorPtr_Type errorElements; // = meshUnstructuredRefined_k->mesh_->getErrorEstimate();
506 MeshUnstrPtr_Type meshUnstructuredRefined_k ;
507 MeshUnstrPtr_Type meshUnstructuredRefined_k_1;
508 MeshUnstrPtr_Type meshUnstructuredRefined_k_m_1;
509 for(int i=0; i<m-1 ; i++){
510 meshUnstructuredRefined_k = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainsP1_[currentIter_-i]->getMesh() , true);
511 meshUnstructuredRefined_k_1 = Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainsP1_[currentIter_-1-i]->getMesh() , true);
512
513 errorElements = errorEstimator.determineCoarseningError(meshUnstructuredRefined_k,meshUnstructuredRefined_k_1,errorEstimationMv_[currentIter_-i],"backwards",markingStrategy_,theta_);
514
515 errorEstimationMv_[currentIter_-1-i]= errorElements;
516
517 }
518
519 // Setting error of: Mesh_(k-m+1) with the previous error ->downscaling errors
520 if(m>1)
521 meshUnstructuredRefined_k_m_1 = meshUnstructuredRefined_k_1;
522 else
523 meshUnstructuredRefined_k_m_1 =Teuchos::rcp_dynamic_cast<MeshUnstr_Type>( domainsP1_[currentIter_]->getMesh() , true);
524
525 // Error of Level l is at l-1
526 for(int i=0; i< n; i++){
527 iterC = k-m+i;
528 if(i==0){
529 errorElements = errorEstimator.determineCoarseningError(meshUnstructuredRefined_k_m_1,meshUnstructuredP1[iterC],errorEstimationMv_[iterC+1],"backwards",markingStrategy_,theta_);
530 }
531 else{
532 errorElements = errorEstimator.determineCoarseningError(meshUnstructuredP1[iterC-1],meshUnstructuredP1[iterC],errorEstimationMv_[iterC-1],"forwards",markingStrategy_,theta_);
533 }
534 if(iterC > errorEstimationMv_.size())
535 errorEstimationMv_.push_back(errorElements);
536 else
537 errorEstimationMv_[iterC]=errorElements;
538
539
540 errorEstimator.markElements(errorElements,theta_,markingStrategy_, meshUnstructuredP1[iterC]);
541
542 refinementFactory.refineMesh(meshUnstructuredP1[iterC],iterC, outputMesh, refinementMode_);
543
544 meshUnstructuredP1[iterC+1] = outputMesh;
545 outputMesh.reset(new MeshUnstr_Type(domainP1->getComm(), inputMeshP1_->volumeID_));
546 }
547
548 outputMesh = meshUnstructuredP1[iterC+1];
549 }
550
551 else if( currentIter_ < maxIter_ ){
552 // Estimating the error with the Discretizations Mesh.
553 errorElementsMv_ = errorEstimator.estimateError(inputMeshP12_, inputMeshP1_, solution, rhsFunc_, domainP12->getFEType());
554
555 errorEstimationMv_.push_back(errorElementsMv_);
556
557 errorEstimator.markElements(errorElementsMv_,theta_,markingStrategy_, inputMeshP1_);
558
559 refinementFactory.refineMesh(inputMeshP1_,currentIter_, outputMesh, refinementMode_);
560 }
561
562
563 // Export distribution of elements among processors
564 MultiVectorPtr_Type procNumTmp = Teuchos::rcp( new MultiVector_Type(domainP12->getElementMap() , 1 ) );
565
566 procNumTmp->putScalar(comm_->getRank());
567 MultiVectorConstPtr_Type vecDecompositionConst = procNumTmp;
568
569 // Error in Nodes
570 MultiVectorConstPtr_Type exactSolution = this->calcExactSolution();
571
572 MultiVectorConstPtr_Type exactSolutionP;
573 MultiVectorConstPtr_Type exportSolutionPMv;
574
575 if( calculatePressure_ ){
576 exportSolutionPMv = problem->getSolution()->getBlock(1);
577 if(exactSolPInput_){
578 exactSolutionP = this->calcExactSolutionP();
579 }
580 }
581
582 calcErrorNorms(exactSolution,solution->getBlock(0), exactSolutionP);
583
584 if(this->exportWithParaview_ && initExporter_==false){
585 this->initExporter( parameterListAll_);
586 }
587 this->exportSolution( inputMeshP12_, problem->getSolution()->getBlock(0), errorNodesMv_, exactSolution, exportSolutionPMv, exactSolutionP);
588 if( currentIter_< maxIter_)
589 this->exportError( inputMeshP12_, errorElementsMv_, errorH1ElementsMv_ , difH1EtaElementsMv_ , vecDecompositionConst );
590
591
592 // Determine all essential values
593 maxErrorEl.push_back(errorElementsMv_->getMax());
594 maxErrorKn.push_back(errorNodesMv_->getMax());
595 numElements.push_back(domainP12_->getElementMap()->getMaxAllGlobalIndex()+1);
596 numElementsProc.push_back(domainP12_->getElementsC()->numberElements());
597
598 numNodes.push_back(domainP12_->getMapUnique()->getMaxAllGlobalIndex()+1);
599
600 if(currentIter_ == maxIter_){
601 writeRefinementInfo();
602 exporterSol_->closeExporter();
603 exporterError_->closeExporter();
604 }
605 else
606 std::cout << " -- done -- " << std::endl;
607
608 domainRefined->setMesh(outputMesh);
609
610 return domainRefined;
611
612}
613
618
619template <class SC, class LO, class GO, class NO>
620typename AdaptiveMeshRefinement<SC,LO,GO,NO>:: MultiVectorConstPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>::calcExactSolution(){
621
622
623 MultiVectorPtr_Type exactSolution = Teuchos::rcp(new MultiVector_Type( solution_->getBlock(0)->getMap() ) );
624 Teuchos::ArrayRCP<SC> exactSolA = exactSolution->getDataNonConst(0);
625
626 vec2D_dbl_ptr_Type points = inputMeshP12_->getPointsUnique();
627
628 Teuchos::ArrayRCP<SC> exactSol(dofs_);
629 for(int i=0; i< points->size(); i++){
630 exactSolFunc_(&points->at(i).at(0),&exactSol[0]);
631 for(int j=0; j< dofs_ ; j++)
632 exactSolA[i*dofs_+j] = exactSol[j];
633
634 }
635
636 MultiVectorConstPtr_Type exactSolConst = exactSolution;
637
638 return exactSolConst;
639}
640
645
646template <class SC, class LO, class GO, class NO>
647typename AdaptiveMeshRefinement<SC,LO,GO,NO>:: MultiVectorConstPtr_Type AdaptiveMeshRefinement<SC,LO,GO,NO>::calcExactSolutionP(){
648
649 MultiVectorPtr_Type exactSolution = Teuchos::rcp(new MultiVector_Type( domainP1_->getMapUnique()));
650 Teuchos::ArrayRCP<SC> exactSolA = exactSolution->getDataNonConst(0);
651
652 vec2D_dbl_ptr_Type points = domainP1_->getPointsUnique();
653
654 Teuchos::ArrayRCP<SC> exactSol(dofsP_);
655 for(int i=0; i< points->size(); i++){
656 exactSolPFunc_(&points->at(i).at(0),&exactSol[0]);
657 exactSolA[i] = exactSol[0];
658 }
659
660 MultiVectorConstPtr_Type exactSolConst = exactSolution;
661
662 return exactSolConst;
663}
664
672
673template <class SC, class LO, class GO, class NO>
674void AdaptiveMeshRefinement<SC,LO,GO,NO>::calcErrorNorms(MultiVectorConstPtr_Type exactSolution, MultiVectorConstPtr_Type solutionP12,MultiVectorConstPtr_Type exactSolutionP){
675
676 // Calculating the error per node
677 MultiVectorPtr_Type errorValues = Teuchos::rcp(new MultiVector_Type( solution_->getBlock(0)->getMap() ) );
678 //this = alpha*A + beta*B + gamma*this
679 errorValues->update( 1., exactSolution, -1. ,solutionP12, 0.);
680
681 // Taking abs norm
682 MultiVectorConstPtr_Type errorValuesAbs = Teuchos::rcp(new MultiVector_Type( solution_->getBlock(0)->getMap()) );
683 errorValuesAbs = errorValues;
684 errorValues->abs(errorValuesAbs);
685
686 // Absolute Error in nodes
687 errorNodesMv_ = errorValues;
688
689 // ---------------------------
690 // Calculating H1 Norm
691 errorH1.push_back(std::sqrt(problem_->calculateH1Norm(errorValues)));
692
693 // ---------------------------
694 // L2 Norm
695 MultiVectorConstPtr_Type exactSolutionTmp = Teuchos::rcp(new MultiVector_Type( domainP12_ ->getMapUnique() ) );
696 Teuchos::ArrayRCP<double > exactSolutionTmpA = exactSolutionTmp->getDataNonConst(0);
697
698 MultiVectorConstPtr_Type solutionTmp = Teuchos::rcp(new MultiVector_Type( domainP12_ ->getMapUnique() ) );
699 Teuchos::ArrayRCP<double > solutionTmpA = solutionTmp->getDataNonConst(0);
700
701 Teuchos::ArrayRCP<double > exactSolutionA = exactSolution->getDataNonConst(0);
702
703 Teuchos::ArrayRCP<double > solutionP12A = solutionP12->getDataNonConst(0);
704
705 double errorL2Tmp=0;
706 for(int i=0; i< dofs_ ; i++){
707
708 MultiVectorPtr_Type errorValues = Teuchos::rcp(new MultiVector_Type( domainP12_ ->getMapUnique() ) );
709 for(int j=0; j< solutionTmpA.size(); j++){
710 solutionTmpA[j] = solutionP12A[j*dofs_+i];
711 exactSolutionTmpA[j] = exactSolutionA[j*dofs_+i];
712 }
713
714 //this = alpha*A + beta*B + gamma*this
715 errorValues->update( 1., exactSolutionTmp, -1. ,solutionTmp, 0.);
716
717 MultiVectorConstPtr_Type errorValuesAbs = Teuchos::rcp(new MultiVector_Type( domainP12_ ->getMapUnique()) );
718 errorValuesAbs = errorValues;
719
720 errorValues->abs(errorValuesAbs);
721
722 errorL2Tmp += problem_->calculateL2Norm(errorValues);
723
724 }
725
726 errorL2.push_back(std::sqrt(errorL2Tmp));
727
728 // -------------------------------
729 // Calculating Error bound epsilon
730 if(exactSolInput_ == true){
731 relError.push_back(std::sqrt(problem_->calculateH1Norm(errorValues)) / std::sqrt(problem_->calculateH1Norm(exactSolution)));
732 }
733
734 if(exactSolInput_ == true){
735 double eta = 0.;
736 Teuchos::ArrayRCP<const double > errorElement = errorElementsMv_->getData(0);
737 for(int i=0; i < errorElement.size() ; i++){
738 //eta += errorElement[i];
739 if(eta < errorElement[i])
740 eta=errorElement[i];
741 }
742 reduceAll<int, double> (*comm_, REDUCE_MAX, eta, outArg (eta));
743 //eta = std::pow(eta,2);
744
745
746 eRelError.push_back(std::sqrt(eta)/ std::sqrt(problem_->calculateH1Norm(solutionP12)));
747 }
748
749 // -------------------------------
750 // Calculating H1 Norm elementwise
751 // First we need a repeated map that is compatible with a vector solution and error
752
753 errorH1ElementsMv_ = Teuchos::rcp( new MultiVector_Type( domainP12_ ->getElementMap(), 1 ) );
754 errorH1ElementsMv_->putScalar(0.);
755 Teuchos::ArrayRCP<SC> errorH1ElementsA = errorH1ElementsMv_->getDataNonConst(0);
756
757 difH1EtaElementsMv_ = Teuchos::rcp( new MultiVector_Type( domainP12_ ->getElementMap(), 1 ) );
758 difH1EtaElementsMv_->putScalar(0.);
759
760 if(domainP12_->getElementMap()->getMaxAllGlobalIndex()< 1500){
761 ElementsPtr_Type elements = domainP12_->getElementsC();
762 MapConstPtr_Type elementMap = domainP12_->getElementMap();
763 MapConstPtr_Type mapUnique = domainP12_->getMapUnique();
764 MapConstPtr_Type mapRep = domainP12_->getMapRepeated();
765
766 vec_GO_Type repIDsVec(0);
767 for(int i=0; i<mapRep->getMaxLocalIndex()+1; i++){
768 GO gID = mapRep->getGlobalElement(i);
769 for(int d=0; d < dofs_ ; d++)
770 repIDsVec.push_back(gID*dofs_+d);
771 }
772 Teuchos::ArrayView<GO> repIDsVecArray = Teuchos::arrayViewFromVector(repIDsVec);
773 // global Ids of Elements' Nodes
774 MapPtr_Type mapRepSystem =
775 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), repIDsVecArray , 0, domainP12_->getComm()) );
776
777 MultiVectorPtr_Type mvValuesError = Teuchos::rcp( new MultiVector_Type( mapRepSystem, 1 ) );
778 Teuchos::ArrayRCP< SC > mvValuesErrorA = mvValuesError->getDataNonConst(0);
779
780 MultiVectorPtr_Type mvValuesErrorUnique = Teuchos::rcp( new MultiVector_Type( solution_->getBlock(0)->getMap(), 1 ) );
781 //Teuchos::ArrayRCP< SC > mvValuesErrorA = mvValuesError->getDataNonConst(0);
782
783
784 MultiVectorPtr_Type errorNodesRep = Teuchos::rcp( new MultiVector_Type( mapRepSystem, 1 ) );
785 Teuchos::ArrayRCP< SC > errorNodesRepA = errorNodesRep->getDataNonConst(0);
786 errorNodesRep->importFromVector(errorNodesMv_,false,"Insert");
787
788
789 for(int k=0; k< elementMap->getMaxAllGlobalIndex()+1; k++){
790 mvValuesError->putScalar(0.);
791 vec_GO_Type notOnMyProc(0);
792 vec_dbl_Type notOnMyProcValue(0);
793 if(elementMap->getLocalElement(k) != -1){
794 vec_int_Type nodeList = elements->getElement(elementMap->getLocalElement(k)).getVectorNodeList();
795 for(int j=0; j< nodeList.size(); j++){
796 for(int d=0; d < dofs_ ; d++){
797 if(mapUnique->getLocalElement(mapRep->getGlobalElement(nodeList[j])) == -1){
798 GO gID = mapRep->getGlobalElement(nodeList[j]);
799 notOnMyProc.push_back(gID*dofs_+d);
800 notOnMyProcValue.push_back(errorNodesRepA[dofs_*nodeList[j]+d]);
801 }
802
803 mvValuesErrorA[dofs_*nodeList[j]+d] = errorNodesRepA[dofs_*nodeList[j]+d];
804 }
805 }
806 }
807 Teuchos::ArrayView<GO> globalNodeArray = Teuchos::arrayViewFromVector( notOnMyProc);
808
809 // global Ids of Elements' Nodes
810 MapPtr_Type mapNodeExport =
811 Teuchos::rcp( new Map_Type(Teuchos::OrdinalTraits<GO>::invalid(), globalNodeArray, 0, domainP12_->getComm()) );
812
813 MultiVectorPtr_Type notMV = Teuchos::rcp( new MultiVector_Type( mapNodeExport, 1 ) );
814 Teuchos::ArrayRCP<SC> notMVA = notMV->getDataNonConst(0);
815 for(int i=0; i< notMVA.size(); i++)
816 notMVA[i] = notOnMyProcValue[i];
817
818 mvValuesErrorUnique->importFromVector(mvValuesError,false,"Insert");
819 mvValuesErrorUnique->importFromVector(notMV,false,"Insert");
820
821 double valueH1 = problem_->calculateH1Norm(mvValuesErrorUnique);
822 double valueL2 = 0; // problem_->calculateL2Norm(mvValuesErrorUnique);
823 if(elementMap->getLocalElement(k) != -1){
824 errorH1ElementsA[elementMap->getLocalElement(k)]= std::sqrt(valueH1 + valueL2);
825 }
826
827 }
828
829 MultiVectorConstPtr_Type errorH1 = errorH1ElementsMv_;
830
831 MultiVectorConstPtr_Type errorElements = errorElementsMv_;
832
833 difH1EtaElementsMv_->update( 1., errorElements , -1. , errorH1, 0.);
834 }
835 if( calculatePressure_== true && exactSolPInput_ == true ){
836 // Calculating the error per node
837 MultiVectorPtr_Type errorValuesP = Teuchos::rcp(new MultiVector_Type( domainP1_->getMapUnique() ) );
838
839 //this = alpha*A + beta*B + gamma*this
840 errorValuesP->update( 1., exactSolutionP, -1. ,solution_->getBlock(1), 0.);
841
842 // Taking abs norm
843 MultiVectorConstPtr_Type errorValuesPAbs = Teuchos::rcp(new MultiVector_Type( domainP1_->getMapUnique() ) );
844 errorValuesPAbs = errorValuesP;
845 errorValuesP->abs(errorValuesPAbs);
846
847 errorNodesPMv_ = errorValuesP;
848
849 double errorL2Tmp = problem_->calculateL2Norm(errorValuesP,1);
850
851 errorL2P.push_back(errorL2Tmp);
852 }
853
854
855
856}
857
858
864template <class SC, class LO, class GO, class NO>
865void AdaptiveMeshRefinement<SC,LO,GO,NO>::initExporter( ParameterListPtr_Type parameterListAll){
866
867 exporterSol_.reset(new ExporterParaViewAMR<SC,LO,GO,NO>());
868 exporterError_.reset(new ExporterParaViewAMR<SC,LO,GO,NO>());
869
870 exporterSol_->setup( "RefinementU" , domainP12_->getMesh(), domainP12_->getFEType(), parameterListAll );
871
872 if(calculatePressure_ ){
873 exporterSolP_.reset(new ExporterParaViewAMR<SC,LO,GO,NO>());
874 exporterSolP_->setup( "RefinementP" , domainP1_->getMesh(), domainP1_->getFEType(), parameterListAll );
875 }
876
877 exporterError_->setup("ErrorEstimation", domainP1_->getMesh(), "P0",parameterListAll );
878
879 initExporter_=true;
880
881}
882
894template <class SC, class LO, class GO, class NO>
895void AdaptiveMeshRefinement<SC,LO,GO,NO>::exportSolution(MeshUnstrPtr_Type mesh, MultiVectorConstPtr_Type exportSolutionMv, MultiVectorConstPtr_Type errorValues, MultiVectorConstPtr_Type exactSolutionMv,MultiVectorConstPtr_Type exportSolutionPMv,MultiVectorConstPtr_Type exactSolutionPMv){
896
897 std::string exporterType = "Scalar";
898 if(dofs_ >1 )
899 exporterType = "Vector";
900
901 if(currentIter_==0){
902
903 exporterSol_->addVariable( exportSolutionMv, "u_h", exporterType, dofs_, domainP12_->getMapUnique() );
904 exporterSol_->addVariable( exactSolutionMv, "u", exporterType, dofs_, domainP12_->getMapUnique() );
905 exporterSol_->addVariable( errorValues, "|u-u_h|", exporterType, dofs_, domainP12_->getMapUnique() );
906
907
908 if( calculatePressure_ ){
909 exporterSolP_->addVariable( exportSolutionPMv, "p_h", "Scalar", dofsP_, domainP1_->getMapUnique() );
910 if(exactSolPInput_){
911 exporterSolP_->addVariable( exactSolutionPMv, "p", "Scalar", dofsP_, domainP1_->getMapUnique() );
912 exporterSolP_->addVariable( errorNodesPMv_, "|p-p_h|", "Scalar", dofsP_, domainP1_->getMapUnique() );
913 }
914 exporterSolP_->save( (double) currentIter_);
915 }
916
917 }
918 else{
919 exporterSol_->reSetup(mesh);
920 exporterSol_->updateVariables(exportSolutionMv, "u_h");
921 exporterSol_->updateVariables( exactSolutionMv, "u" );
922 exporterSol_->updateVariables(errorValues, "|u-u_h|");
923
924 if( calculatePressure_ ){
925 exporterSolP_->reSetup(domainP1_->getMesh());
926 exporterSolP_->updateVariables( exportSolutionPMv, "p_h");
927 if(exactSolPInput_ ){
928 exporterSolP_->updateVariables( exactSolutionPMv, "p");
929 exporterSolP_->updateVariables( errorNodesPMv_, "|p-p_h|");
930 }
931
932 exporterSolP_->save( (double) currentIter_);
933 }
934 }
935
936 exporterSol_->save( (double) currentIter_);
937
938}
939
950template <class SC, class LO, class GO, class NO>
951void AdaptiveMeshRefinement<SC,LO,GO,NO>::exportError(MeshUnstrPtr_Type mesh, MultiVectorConstPtr_Type errorElConst, MultiVectorConstPtr_Type errorElConstH1 , MultiVectorConstPtr_Type difH1Eta ,MultiVectorConstPtr_Type vecDecompositionConst ){
952
953 if(currentIter_==0){
954 exporterError_->addVariable( errorElConst, "eta_T", "Scalar", 1, domainP1_->getElementMap());
955 exporterError_->addVariable( errorElConstH1, "||u-u_h||_H1(T)", "Scalar", 1, domainP1_->getElementMap());
956 exporterError_->addVariable( difH1Eta, "eta_T-||u-u_h||_H1(T)", "Scalar", 1, domainP1_->getElementMap());
957 exporterError_->addVariable( vecDecompositionConst, "Proc", "Scalar", 1, domainP1_->getElementMap());
958 }
959 else{
960 exporterError_->reSetup(mesh);
961
962 exporterError_->updateVariables(errorElConst,"eta_T");
963 exporterError_->updateVariables(errorElConstH1,"||u-u_h||_H1(T)");
964 exporterError_->updateVariables(difH1Eta, "eta_T-||u-u_h||_H1(T)");
965 exporterError_->updateVariables(vecDecompositionConst,"Proc");
966 }
967
968 exporterError_->save((double) currentIter_);
969
970
971}
972
976template <class SC, class LO, class GO, class NO>
977void AdaptiveMeshRefinement<SC,LO,GO,NO>::writeRefinementInfo(){
978
979 using std::cout;
980 using std::endl;
981
982 vec_GO_Type globalProcs(0);
983 for (int i=0; i<= maxRank_; i++)
984 globalProcs.push_back(i);
985
986 Teuchos::ArrayView<GO> globalProcArray = Teuchos::arrayViewFromVector( globalProcs);
987
988 vec_GO_Type localProc(0);
989 localProc.push_back(comm_->getRank());
990 Teuchos::ArrayView<GO> localProcArray = Teuchos::arrayViewFromVector( localProc);
991
992 MapPtr_Type mapGlobalProc =
993 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), globalProcArray, 0, comm_) );
994
995 // Global IDs of Procs
996 MapPtr_Type mapProc =
997 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), localProcArray, 0, comm_) );
998
999 MultiVectorPtr_Type exportLocalEntry = Teuchos::rcp( new MultiVector_Type( mapProc, 1 ) );
1000
1001 exportLocalEntry->putScalar( (LO) numElementsProc[currentIter_] );
1002
1003 MultiVectorPtr_Type elementList= Teuchos::rcp( new MultiVector_Type( mapGlobalProc, 1 ) );
1004 elementList->putScalar( 0 );
1005 elementList->importFromVector( exportLocalEntry, true, "Insert");
1006
1007 Teuchos::ArrayRCP<const double > elementProcList = elementList->getData(0);
1008
1009 double maxNumElementsOnProcs = numElementsProc[currentIter_];
1010 reduceAll<int, double> (*comm_, REDUCE_MAX, maxNumElementsOnProcs, outArg (maxNumElementsOnProcs));
1011
1012 double minNumElementsOnProcs = numElementsProc[currentIter_];
1013 reduceAll<int, double> (*comm_, REDUCE_MIN, minNumElementsOnProcs, outArg (minNumElementsOnProcs));
1014
1015
1016 if(comm_->getRank() == 0 && writeRefinementInfo_ == true){
1017 cout << "__________________________________________________________________________________________________________ " << endl;
1018 cout << " " << endl;
1019 cout << " Summary of Mesh Refinement" << endl;
1020 cout << "__________________________________________________________________________________________________________ " << endl;
1021 cout << " " << endl;
1022 cout << " Marking Strategy:\t" << markingStrategy_ << endl;
1023 cout << " Theta:\t\t\t" << theta_ << endl;
1024 cout << "__________________________________________________________________________________________________________ " << endl;
1025 cout << " " << endl;
1026 cout << " Tolerance:\t\t\t" << tol_ << endl;
1027 cout << " Max number of Iterations:\t" << maxIter_ << endl;
1028 cout << " Number of Processors:\t\t\t" << maxRank_ +1 << endl;
1029 cout << " Number of Refinements:\t\t" << currentIter_ << endl;
1030 cout << "__________________________________________________________________________________________________________ " << endl;
1031 cout << " " << endl;
1032 cout << " Refinementlevel|| Elements\t|| Nodes\t|| Max. estimated error " << endl;
1033 cout << "__________________________________________________________________________________________________________ " << endl;
1034 for(int i=0; i<= currentIter_; i++)
1035 cout <<" "<< i << "\t\t|| " << numElements[i] << "\t\t|| " << numNodes[i]<< "\t\t|| " << maxErrorEl[i]<< endl;
1036 cout << "__________________________________________________________________________________________________________ " << endl;
1037 cout << " " << endl;
1038 if(exactSolInput_ == true){
1039 cout << " Maximal error in nodes after Refinement. " << endl;
1040 for (int i=1; i<=currentIter_ ; i++)
1041 cout <<" "<< i << ":\t" << maxErrorKn[i] << endl;
1042 cout << "__________________________________________________________________________________________________________ " << endl;
1043 cout << " || u-u_h ||_H1\t||\t|| u-u_h ||_L2 ||" ;
1044 if( calculatePressure_== true && exactSolPInput_ == true ){
1045 cout << " \t|| p-p_h||_L2 " << endl;
1046 }
1047 else
1048 cout << endl;
1049 cout << "__________________________________________________________________________________________________________ " << endl;
1050 for (int i=1; i<=currentIter_ ; i++){
1051 std::cout <<" "<< i << ":\t"<< std::setprecision(5) << std::fixed << errorH1[i]<< "\t\t||\t" << errorL2[i] ;
1052 if( calculatePressure_== true && exactSolPInput_ == true ){
1053 std::cout << " \t \t||\t" << std::setprecision(5) << std::fixed << errorL2P[i] << std::endl;
1054 }
1055 else
1056 cout << endl;
1057 }
1058
1059 cout << "__________________________________________________________________________________________________________ " << endl;
1060 }
1061 cout << " ||u-u_h||_H1 / ||u ||_H1 \t|| eta / ||u_h ||_H1\t" << endl;
1062 cout << "__________________________________________________________________________________________________________ " << endl;
1063 for (int i=1; i<=currentIter_ ; i++){
1064 cout <<" "<< i << ":\t" << relError[i] << " \t\t||\t" << eRelError[i] << endl;
1065 }
1066 cout << "__________________________________________________________________________________________________________ " << endl;
1067 cout << " " << endl;
1068 cout << "Distribution of elements on .. " << endl;
1069 //for(int l=0; l< maxRank_ +1 ; l++)
1070 cout <<" Max Number of Elements on Processors " << std::setprecision(0) << std::fixed << maxNumElementsOnProcs << endl;
1071 cout <<" Min Number of Elements on Processors " << minNumElementsOnProcs << endl;
1072 cout << "__________________________________________________________________________________________________________ " << endl;
1073 cout << "__________________________________________________________________________________________________________ " << endl;
1074 cout << " " << endl;
1075 cout << "__________________________________________________________________________________________________________ " << endl;
1076 cout << "__________________________________________________________________________________________________________ " << endl;
1077 cout << " " << endl;
1078 }
1079
1080}
1081
1082
1083
1084}
1085
1086#endif
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
FEDD::AdaptiveMeshRefinement
Definition AdaptiveMeshRefinement_decl.hpp:39
FEDD::AdaptiveMeshRefinement::calcExactSolution
MultiVectorConstPtr_Type calcExactSolution()
Calculating exact solution for velocity if possible with exactSolFunc_.
Definition AdaptiveMeshRefinement_def.hpp:620
FEDD::AdaptiveMeshRefinement::globalAlgorithm
DomainPtr_Type globalAlgorithm(DomainPtr_Type domainP1, DomainPtr_Type domainP12, BlockMultiVectorConstPtr_Type solution, ProblemPtr_Type problem, RhsFunc_Type rhsFunc)
Global Algorithm of Mesh Refinement.
Definition AdaptiveMeshRefinement_def.hpp:408
FEDD::AdaptiveMeshRefinement::refineUniform
DomainPtr_Type refineUniform(DomainPtr_Type domainP1, int level)
Initializing problem if uniform refinement is requested.
Definition AdaptiveMeshRefinement_def.hpp:286
FEDD::AdaptiveMeshRefinement::writeRefinementInfo
void writeRefinementInfo()
Writing refinement information at the end of mesh refinement.
Definition AdaptiveMeshRefinement_def.hpp:977
FEDD::AdaptiveMeshRefinement::exportError
void exportError(MeshUnstrPtr_Type mesh, MultiVectorConstPtr_Type errorElConst, MultiVectorConstPtr_Type errorElConstH1, MultiVectorConstPtr_Type difH1Eta, MultiVectorConstPtr_Type vecDecompositionConst)
ParaViewExporter export of solutions and other error values on current mesh.
Definition AdaptiveMeshRefinement_def.hpp:951
FEDD::AdaptiveMeshRefinement::calcExactSolutionP
MultiVectorConstPtr_Type calcExactSolutionP()
Calculating exact solution for pressure if possible with exactSolPFunc_.
Definition AdaptiveMeshRefinement_def.hpp:647
FEDD::AdaptiveMeshRefinement::exportSolution
void exportSolution(MeshUnstrPtr_Type mesh, MultiVectorConstPtr_Type exportSolutionMv, MultiVectorConstPtr_Type errorValues, MultiVectorConstPtr_Type exactSolutionMv, MultiVectorConstPtr_Type exportSolutionPMv, MultiVectorConstPtr_Type exactSolutionPMv)
ParaViewExporter export of solutions and other error values on current mesh.
Definition AdaptiveMeshRefinement_def.hpp:895
FEDD::AdaptiveMeshRefinement::initExporter
void initExporter(ParameterListPtr_Type parameterListAll)
ParaViewExporter initiation. ParameterListAll contains most settings for ExporterParaView....
Definition AdaptiveMeshRefinement_def.hpp:865
FEDD::AdaptiveMeshRefinement::refineArea
DomainPtr_Type refineArea(DomainPtr_Type domainP1, vec2D_dbl_Type area, int level)
Initializing problem if only a certain area should be refined.
Definition AdaptiveMeshRefinement_def.hpp:240
FEDD::AdaptiveMeshRefinement::identifyProblem
void identifyProblem(BlockMultiVectorConstPtr_Type valuesSolution)
Identifying the problem with respect to the degrees of freedom and whether we calculate pressure....
Definition AdaptiveMeshRefinement_def.hpp:370
FEDD::AdaptiveMeshRefinement::calcErrorNorms
void calcErrorNorms(MultiVectorConstPtr_Type exactSolution, MultiVectorConstPtr_Type solutionP12, MultiVectorConstPtr_Type exactSolutionP)
Calculating error norms. If the exact solution is unknown we use approxmated error norm and error ind...
Definition AdaptiveMeshRefinement_def.hpp:674
FEDD::Domain
Definition Domain_decl.hpp:20
FEDD::ErrorEstimation
Definition ErrorEstimation_decl.hpp:35
FEDD::ErrorEstimation::tagArea
void tagArea(MeshUnstrPtr_Type meshUnstr, vec2D_dbl_Type area)
Tags only a certain Area for refinement and is independent of any error estimation.
Definition ErrorEstimation_def.hpp:473
FEDD::ErrorEstimation::tagAll
void tagAll(MeshUnstrPtr_Type meshUnstr)
Tags only a certain Area for refinement and is independent of any error estimation.
Definition ErrorEstimation_def.hpp:579
FEDD::ErrorEstimation::markElements
void markElements(MultiVectorPtr_Type errorElementMv, double theta, std::string strategy, MeshUnstrPtr_Type meshUnstr)
Function that marks the elements for refinement.
Definition ErrorEstimation_def.hpp:1137
FEDD::ErrorEstimation::determineCoarseningError
MultiVectorPtr_Type determineCoarseningError(MeshUnstrPtr_Type mesh_k, MeshUnstrPtr_Type mesh_k_m, MultiVectorPtr_Type errorElementMv_k, std::string distribution, std::string markingStrategy, double theta)
DetermineCoarseningError is the essential part of the mesh coarsening process.
Definition ErrorEstimation_def.hpp:1907
FEDD::ErrorEstimation::estimateError
MultiVectorPtr_Type estimateError(MeshUnstrPtr_Type inputMeshP12, MeshUnstrPtr_Type inputMeshP1, BlockMultiVectorConstPtr_Type valuesSolution, RhsFunc_Type rhsFunc, std::string FEType)
Main Function for a posteriori error estimation.
Definition ErrorEstimation_def.hpp:156
FEDD::ExporterParaViewAMR
Definition ExporterParaViewAMR_decl.hpp:29
FEDD::RefinementFactory
Definition RefinementFactory_decl.hpp:32
FEDD::RefinementFactory::buildSurfaceTriangleElements
void buildSurfaceTriangleElements(ElementsPtr_Type elements, EdgeElementsPtr_Type edgeElements, SurfaceElementsPtr_Type surfaceTriangleElements, MapConstPtr_Type edgeMap, MapConstPtr_Type elementMap)
Building surface triangle elements, as they are not originally part of the mesh information provided ...
Definition RefinementFactory_def.hpp:824
FEDD::RefinementFactory::refineMesh
void refineMesh(MeshUnstrPtr_Type meshP1, int iteration, MeshUnstrPtr_Type outputMesh, std::string refinementMode)
Main function of RefinementFactory, performs one complete mesh refinement, according to red-green ref...
Definition RefinementFactory_def.hpp:99
FEDD::SurfaceElements
Definition TriangleElements.hpp:18
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36