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