Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
ErrorEstimation_def.hpp
1#ifndef ErrorEstimation_def_hpp
2#define ErrorEstimation_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 "feddlib/core/LinearAlgebra/MultiVector.hpp"
13#include <chrono>
21
22using Teuchos::reduceAll;
23using Teuchos::REDUCE_SUM;
24using Teuchos::REDUCE_MAX;
25using Teuchos::REDUCE_MIN;
26using Teuchos::outArg;
27using std::cout;
28using std::endl;
29using std::setprecision;
30
31namespace FEDD {
32
33
41template <class SC, class LO, class GO, class NO>
42ErrorEstimation<SC,LO,GO,NO>:: ErrorEstimation(int dim, std::string problemType, bool writeMeshQuality)
43{
44 this->dim_ = dim;
45 this->problemType_ = problemType;
46 this->writeMeshQuality_ = writeMeshQuality;
47}
48
49template <class SC, class LO, class GO, class NO>
50ErrorEstimation<SC,LO,GO,NO>::~ErrorEstimation(){
51}
52
59template <class SC, class LO, class GO, class NO>
60void ErrorEstimation<SC,LO,GO,NO>::identifyProblem(BlockMultiVectorConstPtr_Type valuesSolution){
61
62 // dofs can be determined by checking the solutions' blocks an comparing the length of the vectors to the number of unique nodes.
63 // If the length is the same: dofs =1 , if it's twice as long: dofs =2 .. and so on
64
65 // P_12 generally represents the velocity domain:
66 if(inputMesh_->getMapUnique()->getNodeNumElements() == valuesSolution->getBlock(0)->getDataNonConst(0).size())
67 dofs_ = 1;
68 else if(2*(inputMesh_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(0)->getDataNonConst(0).size())
69 dofs_ = 2;
70 else if(3*(inputMesh_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(0)->getDataNonConst(0).size())
71 dofs_ = 3;
72
73 // If ValuesSolution has more than one block its typically because the have also a pressure solution
74 if(valuesSolution->size() > 1){
75 if(inputMeshP1_->getMapUnique()->getNodeNumElements() == valuesSolution->getBlock(1)->getDataNonConst(0).size())
76 dofsP_ = 1;
77 else if(2*(inputMeshP1_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(1)->getDataNonConst(0).size())
78 dofsP_ = 2;
79 else if(3*(inputMeshP1_->getMapUnique()->getNodeNumElements()) == valuesSolution->getBlock(1)->getDataNonConst(0).size())
80 dofsP_ = 3;
81
82 calculatePressure_ = true;
83 }
84 else
85 dofsP_=0;
86
87}
88
89
96
97template <class SC, class LO, class GO, class NO>
98void ErrorEstimation<SC,LO,GO,NO>::makeRepeatedSolution(BlockMultiVectorConstPtr_Type valuesSolution){
99 // We extraxt so solution into dofs-many vectors an use them to make the repeated solution an put it back together in a block multivector later
100 BlockMultiVectorPtr_Type valuesSolutionVel = Teuchos::rcp( new BlockMultiVector_Type(dofs_ ) );
101
102 Teuchos::ArrayRCP< SC > valuesVel = valuesSolution->getBlock(0)->getDataNonConst(0);
103
104 MultiVectorPtr_Type mvValues = Teuchos::rcp( new MultiVector_Type(inputMesh_->getMapUnique(), 1 ) );
105 Teuchos::ArrayRCP< SC > mvValuesA = mvValues->getDataNonConst(0);
106
107 for(int d=0; d< dofs_ ; d++){
108 MultiVectorPtr_Type mvValuesRep = Teuchos::rcp( new MultiVector_Type(inputMesh_->getMapRepeated(), 1 ) );
109 for(int i=0; i< mvValuesA.size(); i++){
110 mvValuesA[i] = valuesVel[i*dofs_+d];
111
112 }
113 mvValuesRep->importFromVector(mvValues,false,"Insert");
114 valuesSolutionVel->addBlock(mvValuesRep,d);
115
116 }
117 if(calculatePressure_ == true){
118
119 BlockMultiVectorPtr_Type valuesSolutionPre = Teuchos::rcp( new BlockMultiVector_Type(dofsP_ ) );
120
121 Teuchos::ArrayRCP< SC > valuesPre = valuesSolution->getBlock(1)->getDataNonConst(0);
122
123 MultiVectorPtr_Type mvValues = Teuchos::rcp( new MultiVector_Type(inputMeshP1_->getMapUnique(), 1 ) );
124 Teuchos::ArrayRCP< SC > mvValuesA = mvValues->getDataNonConst(0);
125
126 for(int d=0; d< dofsP_ ; d++){
127 MultiVectorPtr_Type mvValuesRep = Teuchos::rcp( new MultiVector_Type(inputMeshP1_->getMapRepeated(), 1 ) );
128 for(int i=0; i< mvValuesA.size(); i++){
129 mvValuesA[i] = valuesPre[i*dofsP_+d];
130
131 }
132 mvValuesRep->importFromVector(mvValues,false,"Insert");
133 valuesSolutionPre->addBlock(mvValuesRep,d);
134
135 }
136
137 valuesSolutionRepPre_ = valuesSolutionPre;
138
139 }
140
141 valuesSolutionRepVel_ = valuesSolutionVel;
142
143}
144
155template <class SC, class LO, class GO, class NO>
156typename ErrorEstimation<SC,LO,GO,NO>::MultiVectorPtr_Type ErrorEstimation<SC,LO,GO,NO>::estimateError(MeshUnstrPtr_Type inputMeshP12, MeshUnstrPtr_Type inputMeshP1, BlockMultiVectorConstPtr_Type valuesSolution, RhsFunc_Type rhsFunc, std::string FETypeV){
157
158 // Setting the InputMeshes
159 inputMesh_ = inputMeshP12;
160 inputMeshP1_ = inputMeshP1;
161
162 // Identifying the Problem with respect to the dofs_ and splitting the Solution
163 this->identifyProblem(valuesSolution);
164 this->makeRepeatedSolution(valuesSolution);
165
166
167 // Setting FETypes
168 FEType1_ = "P1"; // Pressure FEType
169 FEType2_ = FETypeV; // Velocity or u FEType
170
171 if(FEType2_ != "P1" && FEType2_ != "P2"){
172 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, "Error Estimation only works for Triangular P1 or P2 Elements");
173 }
174
175
176 ElementsPtr_Type elements = inputMesh_->getElementsC();
177 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
178 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
179 vec2D_dbl_ptr_Type points = inputMesh_->getPointsRepeated();
180 int myRank = inputMesh_->comm_->getRank();
181 surfaceElements_ = inputMesh_->getSurfaceTriangleElements();
182
183 this->dim_ = inputMesh_->dim_;
184
185
186 // 3D Residual Based Error Estimation work when using Surfaces
187 // - New Surface Set
188 // -
189 MESH_TIMER_START(errorEstimation," Error Estimation ");
190
191 edgeElements->matchEdgesToElements(elementMap);
192
193 vec2D_GO_Type elementsOfEdgesGlobal = edgeElements->getElementsOfEdgeGlobal();
194 vec2D_LO_Type elementsOfEdgesLocal = edgeElements->getElementsOfEdgeLocal();
195
196 // Vector size of elements for the resdual based error estimate
197 MultiVectorPtr_Type errorElementMv = Teuchos::rcp(new MultiVector_Type( elementMap) );
198 Teuchos::ArrayRCP<SC> errorElement = errorElementMv->getDataNonConst(0);
199
200 double maxErrorEl;
201 double maxErrorElLoc=0.0;
202
203
204 if(this->dim_ == 2){
205
206 // Edge Numbers of Element
207 vec_int_Type edgeNumbers(3);
208
209
210 // The Jump is over edges or surfaces is calculated beforehand, as it involves communication
211 // In the function itselfe is decided which kind of jump is calculated.
212 vec_dbl_Type u_Jump = calculateJump();
213
214
215 // Calculating diameter of elements
216 vec_dbl_Type areaTriangles(elements->numberElements());
217 vec_dbl_Type rho_T(elements->numberElements());
218 vec_dbl_Type C_T(elements->numberElements());
219 vec_dbl_Type h_T = calcDiamTriangles(elements,points, areaTriangles, rho_T, C_T);
220
221 // The divU Part and residual of the Element are calculated elementwise, as the are independet of other processors
222 double divUElement=0;
223 double resElement=0;
224
225 for(int k=0; k< elements->numberElements();k++){
226 edgeNumbers = edgeElements->getEdgesOfElement(k); // edges of Element k
227 resElement = determineResElement(elements->getElement(k), rhsFunc); // calculating the residual of element k
228 if(problemType_ == "Stokes" || problemType_ == "NavierStokes"){ // If the Problem is a (Navier)Stokes Problem we calculate divU (maybe start a hierarchy
229 divUElement = determineDivU(elements->getElement(k));
230 }
231
232 errorElement[k] = std::sqrt(1./2*(u_Jump[edgeNumbers[0]]+u_Jump[edgeNumbers[1]]+u_Jump[edgeNumbers[2]])+divUElement+ h_T[k]*h_T[k]*resElement ); //divUElement ;
233
234 if(maxErrorElLoc < errorElement[k] )
235 maxErrorElLoc = errorElement[k];
236 }
237
238 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxErrorElLoc, outArg (maxErrorElLoc));
239
240 if( writeMeshQuality_ ){
241 // We asses the Mesh Quality
242 double maxh_T, minh_T;
243 double maxC_T, minC_T;
244 double maxrho_T, minrho_T;
245 double maxArea_T, minArea_T;
246
247 auto it = max_element(h_T.begin(), h_T.end()); //
248 maxh_T = h_T[distance(h_T.begin(), it)];
249 it = min_element(h_T.begin(), h_T.end()); //
250 minh_T = h_T[distance(h_T.begin(), it)];
251 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxh_T, outArg (maxh_T));
252 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minh_T, outArg (minh_T));
253
254 it = max_element(rho_T.begin(), rho_T.end()); //
255 maxrho_T = rho_T[distance(rho_T.begin(), it)];
256 it = min_element(rho_T.begin(), rho_T.end()); //
257 minrho_T = rho_T[distance(rho_T.begin(), it)];
258 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxrho_T, outArg (maxrho_T));
259 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minrho_T, outArg (minrho_T));
260
261 it = max_element(areaTriangles.begin(), areaTriangles.end()); //
262 maxArea_T = areaTriangles[distance(areaTriangles.begin(), it)];
263 it = min_element(areaTriangles.begin(), areaTriangles.end()); //
264 minArea_T = areaTriangles[distance(areaTriangles.begin(), it)];
265 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxArea_T, outArg (maxArea_T));
266 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minArea_T, outArg (minArea_T));
267
268
269 it = max_element(C_T.begin(), C_T.end()); //
270 maxC_T = C_T[distance(C_T.begin(), it)];
271 it = min_element(C_T.begin(), C_T.end()); //
272 minC_T = C_T[distance(C_T.begin(), it)];
273 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxC_T, outArg (maxC_T));
274 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minC_T, outArg (minC_T));
275
276 if(inputMesh_->getComm()->getRank() == 0){
277 cout << "\tMesh Quality Assessment 2D of current mesh" << endl;
278 cout << "\t__________________________________________________________________________________________________________ " << endl;
279 cout << " " << endl;
280 cout << " \tCircumdiameter h_T: \t\t" <<"max. = " << setprecision(5) << maxh_T << "\tmin. = " << setprecision(5) << minh_T << endl;
281 cout << " \tIncircumdiameter rho_T:\t\t" <<"max. = " << setprecision(5) <<maxrho_T << "\tmin. = " <<setprecision(5) << minrho_T << endl;
282 cout << " \tArea of Triangles:\t \t" <<"max. = " << setprecision(5) << maxArea_T << "\tmin. = " << setprecision(5) << minArea_T << endl;
283 cout << " \tShape parameter:\t \t" <<"max. = " << setprecision(5) << maxC_T << "\tmin. = " << setprecision(5) <<minC_T << endl;
284 cout << " \tThe maximal Error of Elements is " << maxErrorElLoc << endl;
285 cout << "\t__________________________________________________________________________________________________________ " << endl;
286 }
287 }
288
289
290
291 }
292
293 else if(this->dim_ == 3){
294
295 this->updateElementsOfSurfaceLocalAndGlobal(edgeElements, this->surfaceElements_);
296 this->buildTriangleMap();
297 this->surfaceElements_->matchSurfacesToElements(elementMap);
298 SurfaceElementsPtr_Type surfaceElements = this->surfaceElements_;
299
300 // Jump per Element over edges
301 vec_dbl_Type errorSurfacesInterior(4);
302 // Edge Numbers of Element
303 vec_int_Type surfaceNumbers(4);
304 // tmp values u_1,u_2 of gradient u of two elements corresponing with surface
305 vec_dbl_Type u1(3),u2(3);
306
307 // gradient of u elementwise
308 vec_dbl_Type u_Jump = calculateJump();
309
310 // h_E,min defined as in "A posteriori error estimation for anisotropic tetrahedral and triangular finite element meshes"
311 // This is the minimal per element, later we also need the adjacent elements' minimal height in order to determine h_E
312 vec_dbl_Type volTetraeder = determineVolTet(elements, points);
313
314 // Calculating diameter of elements
315 vec_dbl_Type areaTriangles(surfaceElements->numberElements());
316 vec_dbl_Type rho_Tri(surfaceElements->numberElements());
317 vec_dbl_Type C_Tri(surfaceElements->numberElements());
318 vec_dbl_Type h_Tri = calcDiamTriangles3D(surfaceElements,points, areaTriangles, rho_Tri, C_Tri);
319
320 vec_dbl_Type h_T = calcDiamTetraeder(elements,points, volTetraeder);
321 vec_dbl_Type rho_T = calcRhoTetraeder(elements,surfaceElements, volTetraeder, areaTriangles);
322
323 // necessary entities
324 vec_dbl_Type p1(3),p2(3); // normal Vector of Surface
325 vec_dbl_Type v_E(3); // normal Vector of Surface
326 double norm_v_E;
327
328 int elTmp1,elTmp2;
329
330
331 // Adjustment for 3D Implememtation
332 // In case we have more than one proc we need to exchange information via the interface.
333 // We design a multivector consisting of u's x and y values, to import and export it easier to other procs only for interface elements
334
335 // The divU Part and residual of the Element are calculated elementwise, as the are independet of other processors
336 double divUElement=0;
337 double resElement=0;
338
339 // Then we determine the jump over the edges, if the element we need for this is not on our proc, we import the solution u
340 for(int k=0; k< elements->numberElements();k++){
341 surfaceNumbers = surfaceElements->getSurfacesOfElement(k); // surfaces of Element k
342
343 resElement = determineResElement(elements->getElement(k), rhsFunc); // calculating the residual of element k
344 if(problemType_ == "Stokes" || problemType_ == "NavierStokes") // If the Problem is a (Navier)Stokes Problem we calculate divU (maybe start a hierarchy
345 divUElement = determineDivU(elements->getElement(k));
346 errorElement[k] = std::sqrt(1./2*(u_Jump[surfaceNumbers[0]]+u_Jump[surfaceNumbers[1]]+u_Jump[surfaceNumbers[2]]+u_Jump[surfaceNumbers[3]])+divUElement +h_T[k]*h_T[k]*resElement );
347 if(maxErrorElLoc < errorElement[k] )
348 maxErrorElLoc = errorElement[k];
349 }
350
351 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxErrorElLoc, outArg (maxErrorElLoc));
352
353 if( writeMeshQuality_ ){
354 double maxh_T, minh_T, maxh_Tri, minh_Tri;
355 double maxC_T, minC_T, maxC_Tri, minC_Tri;
356 double maxrho_T, minrho_T, maxrho_Tri, minrho_Tri;
357 double maxArea_T, minArea_T;
358 double maxVol_T, minVol_T;
359
360 // h_Triangles
361 auto it = max_element(h_Tri.begin(), h_Tri.end()); //
362 maxh_Tri = h_Tri[distance(h_Tri.begin(), it)];
363 it = min_element(h_Tri.begin(), h_Tri.end()); //
364
365 minh_Tri = h_Tri[distance(h_Tri.begin(), it)];
366 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxh_Tri, outArg (maxh_Tri));
367 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minh_Tri, outArg (minh_Tri));
368
369 // h_Tetraeder
370 it = max_element(h_T.begin(), h_T.end()); //
371 maxh_T = h_T[distance(h_T.begin(), it)];
372 it = min_element(h_T.begin(), h_T.end()); //
373
374 minh_T = h_T[distance(h_T.begin(), it)];
375 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxh_T, outArg (maxh_T));
376 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minh_T, outArg (minh_T));
377
378 // rho_Tri
379 it = max_element(rho_Tri.begin(), rho_Tri.end()); //
380 maxrho_Tri = rho_Tri[distance(rho_Tri.begin(), it)];
381 it = min_element(rho_Tri.begin(), rho_Tri.end()); //
382 minrho_Tri = rho_Tri[distance(rho_Tri.begin(), it)];
383 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxrho_Tri, outArg (maxrho_Tri));
384 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minrho_Tri, outArg (minrho_Tri));
385
386 // rho_Tetraeder
387 it = max_element(rho_T.begin(), rho_T.end()); //
388 maxrho_T = rho_T[distance(rho_T.begin(), it)];
389 it = min_element(rho_T.begin(), rho_T.end()); //
390 minrho_T = rho_T[distance(rho_T.begin(), it)];
391 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxrho_T, outArg (maxrho_T));
392 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minrho_T, outArg (minrho_T));
393
394 // Area Triangles
395 it = max_element(areaTriangles.begin(), areaTriangles.end()); //
396 maxArea_T = areaTriangles[distance(areaTriangles.begin(), it)];
397 it = min_element(areaTriangles.begin(), areaTriangles.end()); //
398 minArea_T = areaTriangles[distance(areaTriangles.begin(), it)];
399 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxArea_T, outArg (maxArea_T));
400 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minArea_T, outArg (minArea_T));
401
402 // Volume Tetraeder
403 it = max_element(volTetraeder.begin(), volTetraeder.end()); //
404 maxVol_T = volTetraeder[distance(volTetraeder.begin(), it)];
405
406 it = min_element(volTetraeder.begin(), volTetraeder.end()); //
407 minVol_T = volTetraeder[distance(volTetraeder.begin(), it)];
408
409 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxVol_T, outArg (maxVol_T));
410 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minVol_T, outArg (minVol_T));
411
412 // C_Tri
413 it = max_element(C_Tri.begin(), C_Tri.end()); //
414 maxC_Tri = C_Tri[distance(C_Tri.begin(), it)];
415
416 it = min_element(C_Tri.begin(), C_Tri.end()); //
417 minC_Tri = C_Tri[distance(C_Tri.begin(), it)];
418
419 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxC_Tri, outArg (maxC_Tri));
420 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minC_Tri, outArg (minC_Tri));
421
422 // C_T
423 vec_dbl_Type C_T(elements->numberElements());
424 for(int i=0; i< h_T.size(); i++){
425 C_T[i] = h_T[i] / rho_T[i];
426 }
427 it = max_element(C_T.begin(), C_T.end()); //
428 maxC_T = C_T[distance(C_T.begin(), it)];
429
430 it = min_element(C_T.begin(), C_T.end()); //
431 minC_T = C_T[distance(C_T.begin(), it)];
432
433 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, maxC_T, outArg (maxC_T));
434 reduceAll<int, double> (*inputMesh_->getComm(), REDUCE_MAX, minC_T, outArg (minC_T));
435
436
437 if(inputMesh_->getComm()->getRank() == 0){
438 cout << " \t-- Mesh Quality Assessment 3D of current mesh level -- \t" << endl;
439 cout << "\t__________________________________________________________________________________________________________ " << endl;
440 cout << " " << endl;
441 cout << " \tCircumdiameter h_T:\t\t\t" <<"max. = " << setprecision(5) << maxh_T << " min. = " << setprecision(5) <<minh_T << endl;
442 cout << " \tIncircumdiameter rho_T:\t\t\t" <<"max. = " <<setprecision(5) << maxrho_T << " min. = " << setprecision(5) <<minrho_T << endl;
443 cout << " \tCircumdiameter h_Tri:\t\t\t" <<"max. = " <<setprecision(5) << maxh_Tri << " min. = " << setprecision(5) <<minh_Tri << endl;
444 cout << " \tIncircumdiameter rho_Tri:\t\t" <<"max. = " << setprecision(5) <<maxrho_Tri << " min. = " << setprecision(5) <<minrho_Tri << endl;
445 cout << " \tArea of Triangles: \t\t\t" <<"max. = " << setprecision(5) <<maxArea_T << " min. = " << setprecision(5) <<minArea_T << endl;
446 cout << " \tVolume of Tetraeder: \t\t\t" <<"max. = " << setprecision(5) <<maxVol_T << " min. = " << setprecision(5) <<minVol_T << endl;
447 cout << " \tShape parameter Tetraeder: \t\t" <<"max. = " << setprecision(5) << maxC_T << " min. = " << setprecision(5) <<minC_T << endl;
448 cout << " \tShape parameter Triangles: \t\t" <<"max. = " << setprecision(5) << maxC_Tri << " min. = " << setprecision(5) <<minC_Tri << endl;
449 cout << " \tThe maximal Error of Elements is \t" << maxErrorElLoc << endl;
450 cout << "\t__________________________________________________________________________________________________________ " << endl;
451 }
452 }
453
454
455 }
456
457 errorEstimation_ = errorElementMv;
458
459 MESH_TIMER_STOP(errorEstimation);
460
461 return errorElementMv;
462
463}
464
465
472template <class SC, class LO, class GO, class NO>
473void ErrorEstimation<SC,LO,GO,NO>::tagArea( MeshUnstrPtr_Type inputMeshP1,vec2D_dbl_Type area){
474
475 inputMesh_ = inputMeshP1;
476
477 ElementsPtr_Type elements = inputMesh_->getElementsC();
478 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
479 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
480 int myRank = inputMesh_->comm_->getRank();
481
482 int numberEl = elements->numberElements();
483 int numberPoints =dim_+1;
484
485 vec2D_dbl_ptr_Type vecPoints= inputMesh_->getPointsRepeated();
486
487 int taggedElements=0;
488
489 if(inputMesh_->getComm()->getRank() == 0){
490 cout << "\t__________________________________________________________________________________________________________ " << endl;
491 cout << " " << endl;
492 cout << " \tThe area you requested for Refinement is :" << endl ;
493 cout << "\tx in [" << area[0][0] << ", " << area[0][1] << "] " << endl;
494 if(this->dim_>1)
495 cout << "\ty in [" << area[1][0] << ", " << area[1][1] << "] " << endl;
496 if(this->dim_>2)
497 cout << "\tz in [" << area[2][0] << ", " << area[2][1] << "] " << endl;
498 cout << "\t__________________________________________________________________________________________________________ " << endl;
499 }
500 vec_int_Type edgeNum(6);
501 edgeElements->matchEdgesToElements(elementMap);
502
503 vec_dbl_Type point(this->dim_);
504 int edgeNumber = 3*(this->dim_-1);
505
506 LO p1ID, p2ID;
507 vec_dbl_Type P1,P2;
508
509 for(int i=0; i<numberEl; i++){
510 edgeNum = edgeElements->getEdgesOfElement(i);
511 // Checking whether Nodes of Element are part of the to be tagged area
512 for(int k=0; k<numberPoints; k++){
513 if(vecPoints->at(elements->getElement(i).getNode(k)).at(0) >= area[0][0] && vecPoints->at(elements->getElement(i).getNode(k)).at(0) <= area[0][1]){
514 if(vecPoints->at(elements->getElement(i).getNode(k)).at(1) >= area[1][0] && vecPoints->at(elements->getElement(i).getNode(k)).at(1) <= area[1][1]){
515 if(this->dim_>2){
516 if(vecPoints->at(elements->getElement(i).getNode(k)).at(2) >= area[2][0] && vecPoints->at(elements->getElement(i).getNode(k)).at(2) <= area[2][1]){
517 elements->getElement(i).tagForRefinement();
518 k=numberPoints;
519 taggedElements++;
520 }
521 }
522 else if(this->dim_ == 2){
523 elements->getElement(i).tagForRefinement();
524 k=numberPoints;
525 taggedElements++;
526 }
527
528 }
529 }
530 }
531 for(int k=0; k<edgeNumber ; k++){
532 LO p1ID =edgeElements->getElement(edgeNum[k]).getNode(0);
533 LO p2ID =edgeElements->getElement(edgeNum[k]).getNode(1);
534 P1 = vecPoints->at(p1ID);
535 P2 = vecPoints->at(p2ID);
536
537 for (int d=0; d<this->dim_; d++){
538 point[d]= ( (P1)[d] + (P2)[d] ) / 2.;
539 }
540
541 if(point.at(0) >= area[0][0] && point.at(0) <= area[0][1] && !elements->getElement(i).isTaggedForRefinement()){
542 if(point.at(1) >= area[1][0] && point.at(1) <= area[1][1]){
543 if(this->dim_>2){
544 if(point.at(2) >= area[2][0] && point.at(2) <= area[2][1]){
545 elements->getElement(i).tagForRefinement();
546 taggedElements++;
547 }
548 }
549 else if(this->dim_ == 2){
550 elements->getElement(i).tagForRefinement();
551 taggedElements++;
552 }
553
554 }
555 }
556
557 }
558
559 }
560 reduceAll<int, int> (*inputMesh_->getComm(), REDUCE_MAX, taggedElements, outArg (taggedElements));
561
562 if(inputMesh_->getComm()->getRank()==0){
563 cout << "\t__________________________________________________________________________________________________________ " << endl;
564 cout << " " << endl;
565 cout << "\t With the 'tagArea tool' " << taggedElements << " Elements were tagged for Refinement " << endl;
566 cout << "\t__________________________________________________________________________________________________________ " << endl;
567 }
568
569}
570
577
578template <class SC, class LO, class GO, class NO>
579void ErrorEstimation<SC,LO,GO,NO>::tagAll( MeshUnstrPtr_Type inputMeshP1){
580
581 inputMesh_ = inputMeshP1;
582
583 ElementsPtr_Type elements = inputMesh_->getElementsC();
584 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
585 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
586 int myRank = inputMesh_->comm_->getRank();
587
588 int numberEl = elements->numberElements();
589 int numberPoints =dim_+1;
590
591 vec2D_dbl_ptr_Type vecPoints= inputMesh_->getPointsRepeated();
592
593 int taggedElements=0;
594
595 for(int i=0; i<numberEl; i++){
596 elements->getElement(i).tagForRefinement();
597 taggedElements++;
598 }
599
600 reduceAll<int, int> (*inputMesh_->getComm(), REDUCE_MAX, taggedElements, outArg (taggedElements));
601
602 if(inputMesh_->getComm()->getRank()==0){
603 cout << "\t__________________________________________________________________________________________________________ " << endl;
604 cout << " " << endl;
605 cout << "\tWith tag all:' " << taggedElements << " Elements were tagged for Refinement " << endl;
606 cout << "\t__________________________________________________________________________________________________________ " << endl;
607 }
608
609}
610
611// Tagging all elements adjacent to the flag.
612template <class SC, class LO, class GO, class NO>
613void ErrorEstimation<SC,LO,GO,NO>::tagFlag( MeshUnstrPtr_Type inputMeshP1,int flag){
614
615 inputMesh_ = inputMeshP1;
616
617 ElementsPtr_Type elements = inputMesh_->getElementsC();
618 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
619 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
620 int myRank = inputMesh_->comm_->getRank();
621
622 int numberEl = elements->numberElements();
623 int numberPoints =dim_+1;
624
625 vec2D_dbl_ptr_Type vecPoints= inputMesh_->getPointsRepeated();
626
627 int taggedElements=0;
628
629 for(int i=0; i<numberEl; i++){
630 FiniteElement fe = elements->getElement( i );
631 if(fe.getFlag() == flag){
632 elements->getElement(i).tagForRefinement();
633 taggedElements++;
634 }
635 else{
636 ElementsPtr_Type subEl = fe.getSubElements(); // might be null
637 for (int surface=0; surface<fe.numSubElements(); surface++) {
638 FiniteElement feSub = subEl->getElement( surface );
639 if(feSub.getFlag() == flag){
640 elements->getElement(i).tagForRefinement();
641 taggedElements++;
642 }
643 }
644 }
645 }
646
647 reduceAll<int, int> (*inputMesh_->getComm(), REDUCE_MAX, taggedElements, outArg (taggedElements));
648
649 if(inputMesh_->getComm()->getRank()==0){
650 cout << "\t__________________________________________________________________________________________________________ " << endl;
651 cout << " " << endl;
652 cout << "\tWith tag all:' " << taggedElements << " Elements were tagged for Refinement " << endl;
653 cout << "\t__________________________________________________________________________________________________________ " << endl;
654 }
655
656}
657
661template <class SC, class LO, class GO, class NO>
662vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::calculateJump(){
663
664 int surfaceNumbers = dim_+1 ; // Number of (dim-1) - dimensional surfaces of Element (triangle has 3 edges)
665
666 // Necessary mesh objects
667 ElementsPtr_Type elements = inputMesh_->getElementsC();
668 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
669 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
670 int myRank = inputMesh_->comm_->getRank();
671 SurfaceElementsPtr_Type surfaceTriangleElements = this->surfaceElements_;
672 vec2D_dbl_ptr_Type points = inputMesh_->pointsRep_;
673
674 int numberSurfaceElements=0;
675 if(dim_==2){
676 numberSurfaceElements = edgeElements->numberElements();
677 }
678 else if(dim_==3){
679 numberSurfaceElements = surfaceTriangleElements->numberElements();
680
681 }
682
683 vec_dbl_Type surfaceElementsError(numberSurfaceElements);
684
685 // We calculate the Gradient for all edges or surfaces in direction of N. This is done for the gradient of the velocity u or the pressure p.
686 vec3D_dbl_Type u_El = calcNPhi("Gradient", dofs_, FEType2_);
687
688 // We generally do not need to calculate this part of the jump as the pressure is steady over the edges resulting in a zero pressure Jump.
689 vec3D_dbl_Type p_El;
690 //if(calculatePressure_ == true)
691 //p_El = calcNPhi("None", dofsP_, FEType1_);
692
693 double quadWeightConst =1.;
694
695 if(this->dim_ == 2){
696
697 // necessary entities
698 // calculating diameter of elements
699 vec_dbl_Type p1_2(2); // normal Vector of Surface
700
701 double h_E ; // something with edge
702 vec_dbl_Type v_E(2); // normal Vector of edges
703 double norm_v_E;
704
705 for(int k=0; k< edgeElements->numberElements();k++){
706 // Normalenvektor bestimmen:
707 p1_2[0] = points->at(edgeElements->getElement(k).getNode(0)).at(0) - points->at(edgeElements->getElement(k).getNode(1)).at(0);
708 p1_2[1] = points->at(edgeElements->getElement(k).getNode(0)).at(1) - points->at(edgeElements->getElement(k).getNode(1)).at(1);
709
710 double jumpQuad=0;
711 for(int j=0; j<dofs_ ; j++){
712 for(int i=0; i<u_El[k].size();i++){
713 if( calculatePressure_ == false)
714 jumpQuad += std::pow(u_El[k][i][j],2);
715 if(calculatePressure_ == true){
716 jumpQuad += std::pow( u_El[k][i][j] ,2) ; //u_El[k][i][j] - p_El[k][i][0],2)
717 }
718 }
719 }
720
721 h_E = std::sqrt(std::pow(p1_2[0],2)+std::pow(p1_2[1],2));
722
723 if(this->FEType2_ =="P2")
724 quadWeightConst = h_E /6.;
725 else
726 quadWeightConst = h_E;
727
728 surfaceElementsError[k] =h_E *quadWeightConst*jumpQuad;
729
730
731 }
732 }
733
734 if(this->dim_==3){
735 // Edge Numbers of Element
736 vec_int_Type surfaceElementsIds(surfaceNumbers);
737
738 vec_dbl_Type areaTriangles(numberSurfaceElements);
739 vec_dbl_Type tmpVec1(numberSurfaceElements);
740 vec_dbl_Type tmpVec2(numberSurfaceElements);
741 // necessary entities
742 vec_dbl_Type p1(3),p2(3); // normal Vector of Surface
743 vec_dbl_Type v_E(3); // normal Vector of Surface
744 double norm_v_E;
745
746 vec_dbl_Type h_Tri = calcDiamTriangles3D(surfaceTriangleElements,points, areaTriangles, tmpVec1, tmpVec2);
747
748 for(int k=0; k< surfaceTriangleElements->numberElements();k++){
749 FiniteElement surfaceTmp = surfaceTriangleElements->getElement(k);
750 double jumpQuad=0.;
751 for(int j=0; j<dofs_ ; j++){
752 for(int i=0; i<u_El[k].size();i++){
753 if( calculatePressure_ == false)
754 jumpQuad += std::pow(u_El[k][i][j],2);
755 if(calculatePressure_ == true)
756 jumpQuad += std::pow(u_El[k][i][j],2); // -p_El[k][i][0]
757 }
758 }
759 if(this->FEType2_ =="P2")
760 quadWeightConst = areaTriangles[k];
761 else
762 quadWeightConst = areaTriangles[k];
763
764
765 surfaceElementsError[k] = h_Tri[k]*jumpQuad*quadWeightConst;
766 }
767
768 }
769 return surfaceElementsError;
770
771}
772
781template <class SC, class LO, class GO, class NO>
782vec3D_dbl_Type ErrorEstimation<SC,LO,GO,NO>::calcNPhi(std::string phiDerivative, int dofsSol, std::string FEType){
783
784 int surfaceNumbers = dim_+1 ; // Number of (dim-1) - dimensional surfaces of Element (triangle has 3 edges)
785
786 // necessary mesh objects
787 vec2D_dbl_ptr_Type points = inputMesh_->pointsRep_;
788 ElementsPtr_Type elements = inputMesh_->getElementsC();
789 EdgeElementsPtr_Type edgeElements = inputMesh_->getEdgeElements();
790 MapConstPtr_Type elementMap = inputMesh_->getElementMap();
791 SurfaceElementsPtr_Type surfaceTriangleElements = this->surfaceElements_;
792 vec_int_Type surfaceElementsIds(surfaceNumbers);
793 MapConstPtr_Type surfaceMap;
794
795 // number of surface elements depending on dimension
796 int numberSurfaceElements=0;
797 if(dim_==2){
798 numberSurfaceElements = edgeElements->numberElements();
799 surfaceMap= inputMesh_->edgeMap_;
800 }
801 else if(dim_==3){
802 numberSurfaceElements = surfaceTriangleElements->numberElements();
803 surfaceMap= this->surfaceTriangleMap_;
804 }
805
806 // the int-version of the fe discretisation
807 int intFE = 1;
808 int numNodes= dim_+1;
809 int quadPSize = 1;
810
811 if(FEType2_ == "P2"){
812 quadPSize=3; // Number of Surface Quadrature Points
813 }
814 if(FEType == "P2"){
815 numNodes=6;
816 intFE =2;
817 if(dim_ ==3){
818 numNodes=10;
819 }
820 }
821
822 // We determine u for each Quad Point, so we need to determine how many Points we have
823 vec3D_dbl_Type vec_El(numberSurfaceElements,vec2D_dbl_Type(quadPSize,vec_dbl_Type(dofsSol))); // return value
824 vec3D_dbl_Type vec_El1(numberSurfaceElements,vec2D_dbl_Type(quadPSize,vec_dbl_Type(dofsSol))); // as we calculate a jump over a surface we generally have two solutions for each side
825 vec3D_dbl_Type vec_El2(numberSurfaceElements,vec2D_dbl_Type(quadPSize,vec_dbl_Type(dofsSol)));
826
827 vec3D_dbl_Type vec_El_Exp(numberSurfaceElements,vec2D_dbl_Type(quadPSize,vec_dbl_Type(dofsSol))); // the values of our elements as procs we communicate
828
829 // vectors for late map building
830 vec_GO_Type elementImportIDs(0);
831 vec_GO_Type elementExportIDs(0);
832 vec_LO_Type surfaceIDsLocal(0);
833
834 // gradient of u elementwise
835 vec_LO_Type kn1(numNodes);
836
837 SC detB1;
838 SC absDetB1;
839 SmallMatrix<SC> B1(dim_);
840 SmallMatrix<SC> Binv1(dim_);
841 int index0,index;
842
843 edgeElements->matchEdgesToElements(elementMap);
844
845 // elementsOfSurfaceGlobal and -Local for determining the communication
846 vec2D_GO_Type elementsOfSurfaceGlobal;
847 vec2D_LO_Type elementsOfSurfaceLocal;
848
849 if(dim_==2){
850 elementsOfSurfaceGlobal = edgeElements->getElementsOfEdgeGlobal();
851 elementsOfSurfaceLocal = edgeElements->getElementsOfEdgeLocal();
852 }
853 else if(dim_ ==3){
854 elementsOfSurfaceGlobal = surfaceTriangleElements->getElementsOfSurfaceGlobal();
855 elementsOfSurfaceLocal = surfaceTriangleElements->getElementsOfSurfaceLocal();
856 }
857
858 // the normal vector and its norm
859 vec_dbl_Type v_E(dim_);
860 double norm_v_E=1.;
861
862 // We loop through all surfaces in order to calculate nabla u or p on each surface depending on which element we look at.
863 // The jump is calculated via two surfaces. Consequently we have two values per edge/surface which we either have or need to import/export.
864 for(int k=0; k< numberSurfaceElements;k++){
865 // We only need to calculate the jump for interior egdes/surfaces, which are characetrized by the fact that they are connected to two elements
866 bool interiorElement= false;
867 if(dim_ == 2){
868 if((edgeElements->getElement(k).getFlag() == 10 || edgeElements->getElement(k).getFlag() == 0))//&& elementsOfSurfaceLocal.at(k).size()>1)
869 interiorElement=true;
870 }
871 else if(dim_ == 3){
872 if((surfaceTriangleElements->getElement(k).getFlag() == 10 || surfaceTriangleElements->getElement(k).getFlag() == 0))//&& elementsOfSurfaceLocal.at(k).size()>1)
873 interiorElement=true;
874 }
875 if(interiorElement == true && elementsOfSurfaceLocal.at(k).size() > 1){
876 // Per edge we have depending on discretization quadrature points and weights
877 vec_dbl_Type quadWeights(quadPSize);
878 vec2D_dbl_Type quadPoints;
879
880 if(dim_ == 2){
881 quadPoints = getQuadValues(dim_, FEType2_ , "Surface", quadWeights, edgeElements->getElement(k));
882 v_E.at(0) = points->at(edgeElements->getElement(k).getNode(0)).at(1) - points->at(edgeElements->getElement(k).getNode(1)).at(1);
883 v_E.at(1) = -(points->at(edgeElements->getElement(k).getNode(0)).at(0) - points->at(edgeElements->getElement(k).getNode(1)).at(0));
884 norm_v_E = std::sqrt(std::pow(v_E[0],2)+std::pow(v_E[1],2));
885 }
886 else if(dim_ == 3){
887 quadPoints = getQuadValues(dim_, FEType2_ , "Surface", quadWeights, surfaceTriangleElements->getElement(k));
888 // Normalenvektor bestimmen:
889 vec_dbl_Type p1(3),p2(3);
890 p1[0] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(0) - points->at(surfaceTriangleElements->getElement(k).getNode(1)).at(0);
891 p1[1] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(1) - points->at(surfaceTriangleElements->getElement(k).getNode(1)).at(1);
892 p1[2] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(2) - points->at(surfaceTriangleElements->getElement(k).getNode(1)).at(2);
893
894 p2[0] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(0) - points->at(surfaceTriangleElements->getElement(k).getNode(2)).at(0);
895 p2[1] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(1) - points->at(surfaceTriangleElements->getElement(k).getNode(2)).at(1);
896 p2[2] = points->at(surfaceTriangleElements->getElement(k).getNode(0)).at(2) - points->at(surfaceTriangleElements->getElement(k).getNode(2)).at(2);
897
898 v_E[0] = p1[1]*p2[2] - p1[2]*p2[1];
899 v_E[1] = p1[2]*p2[0] - p1[0]*p2[2];
900 v_E[2] = p1[0]*p2[1] - p1[1]*p2[0];
901
902 norm_v_E = std::sqrt(std::pow(v_E[0],2)+std::pow(v_E[1],2)+std::pow(v_E[2],2));
903 }
904
905 vec_LO_Type elementsIDs(0);
906 // Case that both necessary element are on the same Proc
907 if(elementsOfSurfaceLocal.at(k).at(0) != -1){
908 elementsIDs.push_back(elementsOfSurfaceLocal.at(k).at(0));
909 }
910 if(elementsOfSurfaceLocal.at(k).at(1) != -1){
911 elementsIDs.push_back(elementsOfSurfaceLocal.at(k).at(1));
912 }
913 for(int i=0; i<elementsIDs.size(); i++){
914
915 // We extract the nodes of the elements the surface is connected to
916
917 kn1= elements->getElement(elementsIDs[i]).getVectorNodeListNonConst();
918
919 // Transformation Matrices
920 // We need this inverse Matrix to also transform the quadrature points of our surface back to the reference element
921 index0 = kn1[0];
922 for (int s=0; s<dim_; s++) {
923 index = kn1[s+1];
924 for (int t=0; t<dim_; t++) {
925 B1[t][s] = points->at(index).at(t) -points->at(index0).at(t);
926 }
927 }
928
929 detB1 = B1.computeInverse(Binv1);
930 detB1 = std::fabs(detB1);
931
932 vec2D_dbl_Type quadPointsT1(quadPSize,vec_dbl_Type(dim_));
933 for(int l=0; l< quadPSize; l++){
934 for(int p=0; p< dim_ ; p++){
935 for(int q=0; q< dim_; q++){
936 quadPointsT1[l][p] += Binv1[p][q]* (quadPoints[l][q] - points->at(elements->getElement(elementsIDs[i]).getNode(0)).at(q)) ;
937 }
938 }
939
940 }
941 // We make the distinction between a gradient jump calculation or a simple jump calculation
942 for(int l=0; l< quadPSize; l++){
943 if(phiDerivative == "Gradient"){
944 vec2D_dbl_Type deriPhi1 = gradPhi(dim_, intFE, quadPointsT1[l]);
945 vec2D_dbl_Type deriPhiT1(numNodes,vec_dbl_Type(dim_));
946 for(int q=0; q<dim_; q++){
947 for(int p=0;p<numNodes; p++){
948 for(int s=0; s< dim_ ; s++)
949 deriPhiT1[p][q] += (deriPhi1[p][s]*Binv1[s][q]);
950 }
951 }
952 vec2D_dbl_Type u1_Tmp(dofsSol, vec_dbl_Type(dim_));
953 Teuchos::ArrayRCP< SC > valuesSolRep;
954 for( int d =0; d < dofsSol ; d++){
955 valuesSolRep = valuesSolutionRepVel_->getBlock(d)->getDataNonConst(0);
956 for(int s=0; s<dim_; s++){
957 for(int t=0; t< numNodes ; t++){
958 u1_Tmp[d][s] += valuesSolRep[kn1[t]]*deriPhiT1[t][s] ;
959 }
960 }
961
962 }
963 if(i==0){
964 for( int d =0; d < dofsSol ; d++){
965 for(int s=0; s<dim_; s++){
966 vec_El1[k][l][d] += (u1_Tmp[d][s])*(v_E[s]/norm_v_E);
967 }
968 vec_El1[k][l][d] = std::sqrt(quadWeights[l])* vec_El1[k][l][d];
969 }
970 }
971 else {
972 for( int d =0; d < dofsSol ; d++){
973 for(int s=0; s<dim_; s++){
974 vec_El2[k][l][d] += (u1_Tmp[d][s])*(v_E[s]/norm_v_E);
975 }
976 vec_El2[k][l][d] = std::sqrt(quadWeights[l])* vec_El2[k][l][d];
977 }
978 }
979
980
981 }
982
983 if(phiDerivative == "None"){
984 vec_dbl_Type phiV = phi(dim_, intFE, quadPointsT1[l]);
985 vec_dbl_Type vec_Tmp(dofsSol);
986 Teuchos::ArrayRCP< SC > valuesSolRep;
987 for( int d =0; d < dofsSol ; d++){
988 valuesSolRep = valuesSolutionRepPre_->getBlock(d)->getDataNonConst(0);
989 for(int t=0; t< phiV.size() ; t++){
990 vec_Tmp[d] += valuesSolRep[kn1[t]]*phiV[t] ;
991 }
992 }
993 if(i==0){
994 for( int d =0; d < dofsSol ; d++){
995 for(int s=0; s<dim_; s++){
996 vec_El1[k][l][d] += (vec_Tmp[d])*(v_E[s]/norm_v_E);
997 }
998 vec_El1[k][l][d] = std::sqrt(quadWeights[l])* vec_El1[k][l][d];
999 }
1000 }
1001 else {
1002 for( int d =0; d < dofsSol ; d++){
1003 for(int s=0; s<dim_; s++){
1004 vec_El2[k][l][d] += (vec_Tmp[d])*(v_E[s]/norm_v_E);
1005 }
1006 vec_El2[k][l][d] = std::sqrt(quadWeights[l])* vec_El2[k][l][d];
1007 }
1008 }
1009
1010 }
1011
1012 }
1013 }
1014 // If one of out local elementsOfSurface is -1 it means, that one element connected to the surface is on another processor an we later need to exchange information
1015 if(elementsOfSurfaceLocal.at(k).at(0) == -1 || elementsOfSurfaceLocal.at(k).at(1) == -1){
1016
1017 elementImportIDs.push_back(surfaceMap->getGlobalElement(k));
1018
1019 for(int l=0; l< vec_El1[k].size() ;l++){
1020 for( int d =0; d < dofsSol ; d++)
1021 vec_El_Exp[k][l][d] = vec_El1[k][l][d];
1022 }
1023 }
1024 }
1025 }
1026 inputMesh_->getComm()->barrier();
1027
1028
1029 // We construct map from the previously extracted elementImportIDs
1030 sort(elementImportIDs.begin(),elementImportIDs.end());
1031
1032 vec_GO_Type::iterator ip = unique( elementImportIDs.begin(), elementImportIDs.end());
1033
1034 elementImportIDs.resize(distance(elementImportIDs.begin(), ip));
1035 Teuchos::ArrayView<GO> globalElementArrayImp = Teuchos::arrayViewFromVector( elementImportIDs);
1036
1037 // Map which represents the surface ids that need to import element information
1038 MapPtr_Type mapElementImport =
1039 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), globalElementArrayImp, 0, inputMesh_->getComm()) );
1040
1041
1042 int maxRank = std::get<1>(inputMesh_->rankRange_);
1043 MapPtr_Type mapElementsUnique = mapElementImport;
1044
1045
1046 if(maxRank>0 && mapElementsUnique->getMaxAllGlobalIndex() > 0)
1047 mapElementsUnique = mapElementImport->buildUniqueMap( inputMesh_->rankRange_ );
1048
1049 // In case of a vector values solution we need to import/export dofs-many entries and all of those entries for each of the quadpoints
1050
1051 MultiVectorPtr_Type valuesU_x_expU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1052 MultiVectorPtr_Type valuesU_y_expU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1053 MultiVectorPtr_Type valuesU_z_expU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1054
1055 MultiVectorPtr_Type valuesU_x_impU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1056 MultiVectorPtr_Type valuesU_y_impU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1057 MultiVectorPtr_Type valuesU_z_impU = Teuchos::rcp( new MultiVector_Type( mapElementsUnique, 1 ) );
1058
1059
1060 MultiVectorPtr_Type valuesU_x_imp = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1061 MultiVectorPtr_Type valuesU_y_imp = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1062 MultiVectorPtr_Type valuesU_z_imp = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1063
1064 MultiVectorPtr_Type valuesU_x_impF = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1065 MultiVectorPtr_Type valuesU_y_impF = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1066 MultiVectorPtr_Type valuesU_z_impF = Teuchos::rcp( new MultiVector_Type( mapElementImport, 1 ) );
1067
1068 // Array Pointers which will contain the to be exchanges information
1069 Teuchos::ArrayRCP< SC > entriesU_x_imp = valuesU_x_imp->getDataNonConst(0);
1070 Teuchos::ArrayRCP< SC > entriesU_y_imp = valuesU_y_imp->getDataNonConst(0);
1071 Teuchos::ArrayRCP< SC > entriesU_z_imp = valuesU_z_imp->getDataNonConst(0);
1072
1073 Teuchos::ArrayRCP< SC > entriesU_x_impF = valuesU_x_impF->getDataNonConst(0);
1074 Teuchos::ArrayRCP< SC > entriesU_y_impF = valuesU_y_impF->getDataNonConst(0);
1075 Teuchos::ArrayRCP< SC > entriesU_z_impF = valuesU_z_impF->getDataNonConst(0);
1076
1077 // We exchange values per quad point
1078 for(int i=0; i<quadPSize; i++){
1079 for(int j=0; j<elementImportIDs.size(); j++){
1080 entriesU_x_imp[j] = vec_El_Exp[surfaceMap->getLocalElement(elementImportIDs[j])][i][0] ;
1081 if(dofsSol == 2)
1082 entriesU_y_imp[j] = vec_El_Exp[surfaceMap->getLocalElement(elementImportIDs[j])][i][1] ;
1083 else if(dofsSol == 3)
1084 entriesU_z_imp[j] = vec_El_Exp[surfaceMap->getLocalElement(elementImportIDs[j])][i][2] ;
1085 }
1086
1087 valuesU_x_impU->importFromVector(valuesU_x_imp, false, "Insert");
1088 valuesU_y_impU->importFromVector(valuesU_y_imp, false, "Insert");
1089 valuesU_z_impU->importFromVector(valuesU_z_imp, false, "Insert");
1090
1091 valuesU_x_expU->exportFromVector(valuesU_x_imp, false, "Insert");
1092 valuesU_y_expU->exportFromVector(valuesU_y_imp, false, "Insert");
1093 valuesU_z_expU->exportFromVector(valuesU_z_imp, false, "Insert");
1094
1095 valuesU_x_impF->importFromVector(valuesU_x_impU, false, "Insert");
1096 valuesU_y_impF->importFromVector(valuesU_y_impU, false, "Insert");
1097 valuesU_z_impF->importFromVector(valuesU_z_impU, false, "Insert");
1098
1099 valuesU_x_impF->exportFromVector(valuesU_x_expU, false, "Insert");
1100 valuesU_y_impF->exportFromVector(valuesU_y_expU, false, "Insert");
1101 valuesU_z_impF->exportFromVector(valuesU_z_expU, false, "Insert");
1102
1103 for(int j=0; j<elementImportIDs.size(); j++){
1104 vec_El2[surfaceMap->getLocalElement(elementImportIDs[j])][i][0] =entriesU_x_impF[j];
1105 if(dofsSol == 2)
1106 vec_El2[surfaceMap->getLocalElement(elementImportIDs[j])][i][1] =entriesU_y_impF[j];
1107 else if(dofsSol ==3)
1108 vec_El2[surfaceMap->getLocalElement(elementImportIDs[j])][i][2] = entriesU_z_impF[j];
1109 }
1110 }
1111 for(int i=0; i< numberSurfaceElements ; i++){
1112 for(int j=0; j<quadPSize; j++){
1113 for(int k=0; k< dofsSol ; k++){
1114 vec_El[i][j][k] = std::fabs(vec_El1[i][j][k]) - std::fabs(vec_El2[i][j][k]);
1115
1116 }
1117 }
1118 }
1119
1120
1121 return vec_El;
1122
1123}
1124
1135
1136template <class SC, class LO, class GO, class NO>
1137void ErrorEstimation<SC,LO,GO,NO>::markElements(MultiVectorPtr_Type errorElementMv, double theta, std::string strategy, MeshUnstrPtr_Type meshP1){
1138
1139
1140 Teuchos::ArrayRCP<SC> errorEstimation = errorElementMv->getDataNonConst(0);
1141
1142 ElementsPtr_Type elements = meshP1->getElementsC();
1143
1144 this->markingStrategy_ = strategy;
1145
1146 theta_ = theta;
1147 // As we decide which element to tag based on the maximum error in the elements globally, we need to communicated this maxErrorEl
1148 auto it = std::max_element(errorEstimation.begin(), errorEstimation.end()); // c++11
1149 double maxErrorEl= errorEstimation[std::distance(errorEstimation.begin(), it)];
1150
1151 // Maximum-strategy for marking the elements as proposed by Verfürth in "A posterio error estimation"
1152 reduceAll<int, double> (*meshP1->getComm(), REDUCE_MAX, maxErrorEl, outArg (maxErrorEl));
1153 int flagCount=0;
1154 if( markingStrategy_ == "Maximum"){
1155 for(int k=0; k< elements->numberElements() ; k++){
1156 if( errorEstimation[k] > theta_ * maxErrorEl){
1157 elements->getElement(k).tagForRefinement();
1158 flagCount ++;
1159 }
1160 }
1161 }
1162 // Equilibirum/Doerfler-strategy for marking the elements as proposed by Verfürth in "A posterio error estimation"
1163 else if(markingStrategy_ == "Doerfler"){
1164 double thetaSumTmp=0.;
1165 double thetaSum=0.;
1166 double muMax=0.;
1167 double sigmaT=0.;
1168 vec_bool_Type flagged(elements->numberElements());
1169 for(int k=0; k< elements->numberElements(); k++){
1170 thetaSumTmp = thetaSumTmp + std::pow(errorEstimation[k],2);
1171 flagged[k]=false;
1172 }
1173 reduceAll<int, double> (*meshP1->getComm(), REDUCE_SUM, thetaSumTmp, outArg (thetaSum));
1174 while(sigmaT < theta_*thetaSum){
1175 muMax=0.;
1176 for(int k=0; k< elements->numberElements(); k++){
1177 if(muMax < errorEstimation[k] && flagged[k]==false ){
1178 muMax = errorEstimation[k];
1179 }
1180 }
1181 reduceAll<int, double> (*meshP1->getComm(), REDUCE_MAX, muMax, outArg (muMax));
1182 for(int k=0; k< elements->numberElements(); k++){
1183 if(muMax == errorEstimation[k] && flagged[k]==false ){
1184 flagged[k]=true;
1185 sigmaT = sigmaT + std::pow(errorEstimation[k],2);
1186 }
1187 }
1188 reduceAll<int, double> (*meshP1->getComm(), REDUCE_MAX, sigmaT, outArg (sigmaT));
1189 }
1190
1191 for(int k=0; k< elements ->numberElements() ; k++){
1192 if( flagged[k]==true){
1193 elements->getElement(k).tagForRefinement();
1194 flagCount++;
1195 }
1196 }
1197
1198 }
1199
1200 // If no strategy is chosen or we choose uniform refinement ourselves, all elements are marked for refinement
1201 else{
1202 for(int k=0; k< elements->numberElements() ; k++){
1203 elements->getElement(k).tagForRefinement();
1204 flagCount++;
1205 }
1206 }
1207 reduceAll<int, int> (*meshP1->getComm(), REDUCE_SUM, flagCount , outArg (flagCount));
1208
1209 if(meshP1->getComm()->getRank() == 0){
1210 cout << " " << endl;
1211 cout << " \tA-posteriori Error Estimation for " << problemType_ << "-problem tagged " << flagCount << " Elements for adaptive Refinement with " << markingStrategy_ << "-Strategy and Theta= " << theta_ << endl;
1212 cout << " " << endl;
1213 }
1214
1215}
1216
1225template <class SC, class LO, class GO, class NO>
1226double ErrorEstimation<SC,LO,GO,NO>::determineDivU(FiniteElement element){
1227
1228// Quad Points and weights of third order
1229 double divElement =0.;
1230
1231 int dim = this->dim_;
1232 SC* paras ; //= &(funcParameter[0]);
1233
1234 int t=1;
1235 if(dim == 2)
1236 t =3;
1237 else if(dim == 3)
1238 t=5;
1239
1240 vec_dbl_Type QuadW(t);
1241 vec2D_dbl_Type QuadPts = getQuadValues(dim, FEType2_, "Element", QuadW, element);
1242
1243
1244 vec2D_dbl_ptr_Type points = inputMesh_->getPointsRepeated();
1245
1246 // Transformation Matrices
1247 // We determine deltaU for the Elements. If FEType=P1 deltaU equals 0. If FEType = P2 we get a constant solution:
1248 double deltaU =0.;
1249 vec_LO_Type nodeList = element.getVectorNodeListNonConst();
1250
1251 SC detB1;
1252 SmallMatrix<SC> B1(dim);
1253 SmallMatrix<SC> Binv1(dim);
1254
1255 // We need this inverse Matrix to also transform the quad Points back to the reference element
1256 int index0 = nodeList[0];
1257 for (int s=0; s<dim; s++) {
1258 int index = nodeList[s+1];
1259 for (int t=0; t<dim; t++) {
1260 B1[t][s] = points->at(index).at(t) - points->at(index0).at(t);
1261 }
1262 }
1263
1264
1265
1266 detB1 = B1.computeInverse(Binv1);
1267 detB1 = std::fabs(detB1);
1268
1269 int intFE = 1;
1270 if(this->FEType2_ == "P2")
1271 intFE =2;
1272
1273
1274 double divElementTmp=0.;
1275 for (UN w=0; w< QuadW.size(); w++){
1276 vec2D_dbl_Type gradPhiV =gradPhi(dim, intFE, QuadPts[w]);
1277 vec2D_dbl_Type gradPhiT(gradPhiV.size(),vec_dbl_Type(dim));
1278 for(int q=0; q<dim; q++){
1279 for(int p=0;p<nodeList.size(); p++){
1280 for(int s=0; s< dim ; s++)
1281 gradPhiT[p][q] += (gradPhiV[p][s]*Binv1[s][q]);
1282 }
1283 }
1284 vec_dbl_Type uTmp(gradPhiV.size());
1285 for(int j=0; j< dofs_ ; j++){
1286 Teuchos::ArrayRCP<SC> valuesSolutionRep = valuesSolutionRepVel_->getBlock(j)->getDataNonConst(0);
1287 for(int i=0; i< nodeList.size(); i++){
1288 uTmp[i] += valuesSolutionRep[nodeList[i]]*gradPhiT[i][j];
1289 }
1290 }
1291 divElementTmp=0.;
1292 for(int i=0; i< nodeList.size(); i++){
1293 divElementTmp += uTmp[i];
1294 }
1295 divElementTmp = std::pow(divElementTmp,2);
1296 divElementTmp *= QuadW[w];
1297
1298 divElement += divElementTmp;
1299 }
1300 divElement = divElement *detB1;
1301
1302
1303
1304 return divElement;
1305
1306
1307}
1308
1317
1318template <class SC, class LO, class GO, class NO>
1319double ErrorEstimation<SC,LO,GO,NO>::determineResElement(FiniteElement element, RhsFunc_Type rhsFunc){
1320
1321// Quad Points and weights of third order
1322
1323 vec_dbl_Type resElement(dofs_);
1324
1325 int dim = this->dim_;
1326 SC* paras ;
1327
1328 int t=1;
1329 if(dim == 2)
1330 t =3;
1331 else if(dim == 3)
1332 t=5;
1333
1334 vec_dbl_Type QuadW(t);
1335 vec2D_dbl_Type QuadPts = getQuadValues(dim, FEType2_, "Element", QuadW, element);
1336
1337 int intFE = 1;
1338 if(this->FEType2_ == "P2")
1339 intFE =2;
1340
1341 vec2D_dbl_ptr_Type points = inputMesh_->getPointsRepeated();
1342
1343 // Transformation Matrices
1344 // We determine deltaU for the Elements. If FEType=P1 deltaU equals 0. If FEType = P2 we get a constant solution:
1345 vec_LO_Type nodeList = element.getVectorNodeListNonConst();
1346
1347 SC detB1;
1348 SmallMatrix<SC> B1(dim);
1349 SmallMatrix<SC> Binv1(dim);
1350
1351 // We need this inverse Matrix to also transform the quad Points of on edge back to the reference element
1352 int index0 = nodeList[0];
1353 for (int s=0; s<dim; s++) {
1354 int index = nodeList[s+1];
1355 for (int t=0; t<dim; t++) {
1356 B1[t][s] = points->at(index).at(t) - points->at(index0).at(t);
1357 }
1358 }
1359
1360
1361 detB1 = B1.computeInverse(Binv1);
1362 detB1 = std::fabs(detB1);
1363
1364 vec_dbl_Type valueFunc(dofs_);
1365
1366 vec2D_dbl_Type quadPointsTrans(QuadW.size(),vec_dbl_Type(dim));
1367
1368 for(int i=0; i< QuadW.size(); i++){
1369 for(int j=0; j< dim; j++){
1370 for(int k=0; k< dim; k++){
1371 quadPointsTrans[i][j] += B1[j][k]* QuadPts[i].at(k) ;
1372 }
1373 quadPointsTrans[i][j] += points->at(element.getNode(0)).at(j);
1374 }
1375 }
1376
1377 vec_dbl_Type nablaP(dim);
1378 if(problemType_ == "Stokes" || problemType_ == "NavierStokes"){
1379 Teuchos::ArrayRCP<SC> valuesSolutionRepPre = valuesSolutionRepPre_->getBlock(0)->getDataNonConst(0);
1380 vec2D_dbl_Type deriPhi = gradPhi(dim, 1, QuadPts[0]);
1381 vec2D_dbl_Type deriPhiT(dim+1,vec_dbl_Type(dim));
1382 for(int q=0; q<dim; q++){
1383 for(int p=0;p<dim+1; p++){
1384 for(int s=0; s< dim ; s++)
1385 deriPhiT[p][q] += (deriPhi[p][s]*Binv1[s][q]);
1386
1387 }
1388 }
1389 for(int s=0; s<dim; s++){
1390 for(int t=0; t< dim+1 ; t++){
1391 nablaP[s] += valuesSolutionRepPre[nodeList[t]]*deriPhiT[t][s] ;
1392 }
1393 }
1394 }
1395 // ####################################################################
1396 vec2D_dbl_Type nablaU(dim,vec_dbl_Type(QuadW.size()));
1397 if(problemType_ == "NavierStokes"){
1398 for (UN w=0; w< QuadW.size(); w++){
1399 vec2D_dbl_Type deriPhi = gradPhi(dim, intFE, QuadPts[w]);
1400 vec2D_dbl_Type deriPhiT(nodeList.size(),vec_dbl_Type(dim));
1401 for(int q=0; q<dim; q++){
1402 for(int p=0;p<nodeList.size(); p++){
1403 for(int s=0; s< dim ; s++)
1404 deriPhiT[p][q] += (deriPhi[p][s]*Binv1[s][q]);
1405
1406 }
1407 }
1408 Teuchos::ArrayRCP< SC > valuesSolRep;
1409 vec2D_dbl_Type vec_Tmp(nodeList.size(),vec_dbl_Type(dofs_));
1410 for( int d =0; d < dofs_ ; d++){
1411 valuesSolRep = valuesSolutionRepVel_->getBlock(d)->getDataNonConst(0);
1412 for(int t=0; t< nodeList.size() ; t++){
1413 vec_Tmp[t][d] = valuesSolRep[nodeList[t]];
1414 }
1415 }
1416
1417 vec2D_dbl_Type u1_Tmp(dofs_, vec_dbl_Type(dim_));
1418
1419 for( int d =0; d < dofs_ ; d++){
1420 valuesSolRep = valuesSolutionRepVel_->getBlock(d)->getDataNonConst(0);
1421 for(int s=0; s<dim; s++){
1422 for(int t=0; t< nodeList.size() ; t++){
1423 u1_Tmp[d][s] += valuesSolRep[nodeList[t]]*deriPhiT[t][s] ;
1424 }
1425 }
1426 }
1427 vec_dbl_Type phiV = phi(dim, intFE, QuadPts[w]);
1428 for( int d =0; d < dofs_ ; d++){
1429 for(int t=0; t< nodeList.size() ; t++){
1430 for(int s=0; s<dim; s++){
1431 nablaU[d][w] += vec_Tmp[t][s]*phiV[t]*u1_Tmp[s][d];
1432 }
1433 }
1434 }
1435
1436 }
1437 }
1438
1439 vec_dbl_Type deltaU(dofs_);
1440 if(this->FEType2_ == "P2"){
1441 vec_dbl_Type deltaPhi(nodeList.size());
1442 if(this->dim_ == 2){
1443 deltaPhi={8, 4, 4, -8, 0, -8 };
1444 }
1445 else if(this->dim_ == 3)
1446 deltaPhi={12, 4, 4, 4, -8, 0, -8, -8, 0, 0};
1447
1448 for(int j=0 ; j< dofs_; j++){
1449 Teuchos::ArrayRCP<SC> valuesSolutionRep = valuesSolutionRepVel_->getBlock(j)->getDataNonConst(0);
1450 for(int i=0; i< nodeList.size() ; i++){
1451 deltaU[j] += deltaPhi[i]*valuesSolutionRep[nodeList[i]];
1452 }
1453 }
1454 }
1455 for (UN w=0; w< QuadW.size(); w++){
1456 rhsFunc(&quadPointsTrans[w][0], &valueFunc[0] ,paras);
1457 for(int j=0 ; j< dofs_; j++){
1458 resElement[j] += QuadW[w] * std::pow(valueFunc[j] + deltaU[j] - nablaU[j][w] - nablaP[j] ,2);
1459 }
1460 }
1461 double resElementValue =0.;
1462 for(int j=0; j< dofs_ ; j++)
1463 resElementValue += resElement[j] *detB1;
1464
1465 return resElementValue;
1466
1467
1468}
1469
1484template <class SC, class LO, class GO, class NO>
1485vec2D_dbl_Type ErrorEstimation<SC,LO,GO,NO>::getQuadValues(int dim, std::string FEType, std::string Type, vec_dbl_Type &QuadW, FiniteElement surface){
1486
1487 vec2D_dbl_Type QuadPts(QuadW.size(), vec_dbl_Type(dim));
1488 vec2D_dbl_ptr_Type points = inputMesh_->pointsRep_;
1489 if(Type == "Element"){
1490 if(this->dim_ == 2){
1491
1492 double a = 1/6.;
1493 QuadPts[0][0] = 0.5;
1494 QuadPts[0][1] = 0.5;
1495
1496 QuadPts[1][0] = 0.;
1497 QuadPts[1][1] = 0.5;
1498
1499 QuadPts[2][0] = 0.5;
1500 QuadPts[2][1] = 0.;
1501
1502 QuadW[0] = a;
1503 QuadW[1] = a;
1504 QuadW[2] = a;
1505 }
1506 else if(this->dim_ == 3){
1507
1508
1509 double a = .25;
1510 double b = 1./6.;
1511 double c = .5;
1512 QuadPts[0][0] = a;
1513 QuadPts[0][1] = a;
1514 QuadPts[0][2]= a;
1515
1516 QuadPts[1][0] = b;
1517 QuadPts[1][1] = b;
1518 QuadPts[1][2] = b;
1519
1520 QuadPts[2][0] = b;
1521 QuadPts[2][1] = b;
1522 QuadPts[2][2]= c;
1523
1524 QuadPts[3][0] = b;
1525 QuadPts[3][1] = c;
1526 QuadPts[3][2] = b;
1527
1528 QuadPts[4][0] = c;
1529 QuadPts[4][1] = b;
1530 QuadPts[4][2] = b;
1531
1532 QuadW[0] = -2./15.;
1533 QuadW[1] = 3./40.;
1534 QuadW[2] = 3./40.;
1535 QuadW[3] = 3./40.;
1536 QuadW[4] = 3./40.;
1537 }
1538 }
1539 else if (Type =="Surface"){
1540 if(dim==2){
1541 double x0 = points->at(surface.getNode(0)).at(0);
1542 double y0 = points->at(surface.getNode(0)).at(1);
1543 double x1 = points->at(surface.getNode(1)).at(0);
1544 double y1 = points->at(surface.getNode(1)).at(1);
1545
1546
1547 if(FEType == "P1"){
1548
1549 QuadPts[0][0] = (x0+x1)/2.;
1550 QuadPts[0][1] = (y0+y1)/2.;
1551
1552 QuadW[0] = 1.;
1553 }
1554 else if(FEType == "P2"){
1555
1556 QuadPts[0][0] = x0;
1557 QuadPts[0][1] = y0;
1558 QuadPts[1][0] = (x0+x1)/2.;
1559 QuadPts[1][1] = (y0+y1)/2.;
1560 QuadPts[2][0] = x1;
1561 QuadPts[2][1] = y1;
1562
1563 QuadW[0] = 1.;
1564 QuadW[1] = 4.;
1565 QuadW[2] = 1.;
1566 }
1567
1568 }
1569 else if(dim==3){
1570 // Here we choose as quadpoints the midpoints of the triangle sides
1571 double x0 = points->at(surface.getNode(0)).at(0);
1572 double y0 = points->at(surface.getNode(0)).at(1);
1573 double z0 = points->at(surface.getNode(0)).at(2);
1574 double x1 = points->at(surface.getNode(1)).at(0);
1575 double y1 = points->at(surface.getNode(1)).at(1);
1576 double z1 = points->at(surface.getNode(1)).at(2);
1577 double x2 = points->at(surface.getNode(2)).at(0);
1578 double y2 = points->at(surface.getNode(2)).at(1);
1579 double z2 = points->at(surface.getNode(2)).at(2);
1580
1581 if(FEType == "P1"){
1582 // As nabla phi is a constant function, quad points don't really matter in that case
1583 QuadPts[0][0] = 1/3.;
1584 QuadPts[0][1] = 1/3.;
1585 QuadPts[0][2] = 1/3.;
1586
1587 QuadW[0] = 1.;
1588 }
1589 else if(FEType == "P2"){
1590 QuadPts[0][0] = (x0+x1)/2.;
1591 QuadPts[0][1] = (y0+y1)/2.;
1592 QuadPts[0][2] = (z0+z1)/2.;
1593 QuadPts[1][0] = (x0+x2)/2.;
1594 QuadPts[1][1] = (y0+y2)/2.;
1595 QuadPts[1][2] = (z0+z2)/2.;
1596 QuadPts[2][0] = (x1+x2)/2.;
1597 QuadPts[2][1] = (y1+y2)/2.;
1598 QuadPts[2][2] = (z1+z2)/2.;
1599
1600 QuadW[0] = 1/3.;
1601 QuadW[1] = 1/3.;
1602 QuadW[2] = 1/3.;
1603 }
1604
1605 }
1606 }
1607
1608 return QuadPts;
1609
1610}
1611
1620
1621template <class SC, class LO, class GO, class NO>
1622vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::determineVolTet(ElementsPtr_Type elements,vec2D_dbl_ptr_Type points){
1623
1624 vec_dbl_Type volumeTetrahedra(elements->numberElements());
1625
1626 vec_dbl_Type p1(3),p2(3), p3(3), v_K(3);
1627
1628 vec2D_dbl_Type p(4,vec_dbl_Type(3));
1629
1630 for(int k=0; k< elements->numberElements() ; k++){
1631
1632 // Calculating edges of Tetraeder
1633 vec_LO_Type nodeList = elements->getElement(k).getVectorNodeListNonConst();
1634
1635 for(int i=0; i<4; i++){
1636 p[i] = points->at(nodeList[i]);
1637 }
1638
1639 p1[0] = p[0][0] - p[1][0];
1640 p1[1] = p[0][1] - p[1][1];
1641 p1[2] = p[0][2] - p[1][2];
1642
1643 p2[0] = p[0][0] - p[2][0];
1644 p2[1] = p[0][1] - p[2][1];
1645 p2[2] = p[0][2] - p[2][2];
1646
1647 p3[0] = p[0][0] - p[3][0];
1648 p3[1] = p[0][1] - p[3][1];
1649 p3[2] = p[0][2] - p[3][2];
1650
1651
1652 v_K[0] = p1[1]*p2[2] - p1[2]*p2[1];
1653 v_K[1] = p1[2]*p2[0] - p1[0]*p2[2];
1654 v_K[2] = p1[0]*p2[1] - p1[1]*p2[0];
1655
1656
1657 volumeTetrahedra[k] = std::fabs(p3[0] * v_K[0] + p3[1] * v_K[1] +p3[2] * v_K[2]) / 6. ;
1658 }
1659
1660
1661
1662 return volumeTetrahedra;
1663}
1664
1675
1676template <class SC, class LO, class GO, class NO>
1677vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::calcDiamTetraeder(ElementsPtr_Type elements,vec2D_dbl_ptr_Type points, vec_dbl_Type volTet){
1678
1679 vec_dbl_Type diamElements(elements->numberElements());
1680
1681 double a,b,c,A,B,C;
1682
1683 vec2D_dbl_Type p(4,vec_dbl_Type(this->dim_));
1684 for(int k=0; k< elements->numberElements() ; k++){
1685
1686 // Calculating edges of Tetraeder
1687 vec_LO_Type nodeList = elements->getElement(k).getVectorNodeListNonConst();
1688
1689 for(int i=0; i<4; i++){
1690 p[i] = points->at(nodeList[i]);
1691 }
1692
1693 a = std::sqrt(std::pow(p[0][0] - p[1][0],2)+ std::pow(p[0][1] - p[1][1],2) +std::pow(p[0][2] - p[1][2],2));
1694 b = std::sqrt(std::pow(p[0][0] - p[2][0],2)+ std::pow(p[0][1] - p[2][1],2) +std::pow(p[0][2] - p[2][2],2));
1695 c = std::sqrt(std::pow(p[0][0] - p[3][0],2)+ std::pow(p[0][1] - p[3][1],2) +std::pow(p[0][2] - p[3][2],2));
1696
1697 A = std::sqrt(std::pow(p[3][0] - p[2][0],2)+ std::pow(p[3][1] - p[2][1],2) +std::pow(p[3][2] - p[2][2],2));
1698 B = std::sqrt(std::pow(p[1][0] - p[3][0],2)+ std::pow(p[1][1] - p[3][1],2) +std::pow(p[1][2] - p[3][2],2));
1699 C = std::sqrt(std::pow(p[1][0] - p[2][0],2)+ std::pow(p[1][1] - p[2][1],2) +std::pow(p[1][2] - p[2][2],2));
1700
1701
1702 diamElements[k] = std::sqrt(((a*A+b*B+c*C)*(-a*A+b*B+c*C)*(a*A-b*B+c*C)*(a*A+b*B-c*C))) / (12*volTet[k]) ;
1703
1704 }
1705
1706
1707 return diamElements;
1708}
1709
1721
1722template <class SC, class LO, class GO, class NO>
1723vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::calcRhoTetraeder(ElementsPtr_Type elements,SurfaceElementsPtr_Type surfaceTriangleElements, vec_dbl_Type volTet, vec_dbl_Type areaTriangles){
1724 vec_dbl_Type rhoElements(elements->numberElements());
1725
1726
1727 vec_LO_Type surfaceOfEl(4);
1728
1729 for(int k=0; k< elements->numberElements() ; k++){
1730 // Calculating edges of Tetraeder
1731 surfaceOfEl = surfaceTriangleElements->getSurfacesOfElement(k);
1732
1733 rhoElements[k] = (6*volTet[k]) / (areaTriangles[surfaceOfEl[0]] + areaTriangles[surfaceOfEl[1]] +areaTriangles[surfaceOfEl[2]] +areaTriangles[surfaceOfEl[3]]);
1734
1735 }
1736
1737
1738 return rhoElements;
1739}
1740
1749
1750template <class SC, class LO, class GO, class NO>
1751vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::calcDiamTriangles3D(SurfaceElementsPtr_Type surfaceTriangleElements,vec2D_dbl_ptr_Type points,vec_dbl_Type& areaTriangles, vec_dbl_Type& rho_T, vec_dbl_Type& C_T){
1752
1753 vec_dbl_Type diamElements(surfaceTriangleElements->numberElements());
1754
1755 FiniteElement elementTmp;
1756
1757 vec_dbl_Type vecTmp(3),vecTmp1(3),vecTmp2(3);
1758 for(int k=0; k< surfaceTriangleElements->numberElements() ; k++){
1759 double lengthA, lengthB, lengthC,s1;
1760
1761 elementTmp = surfaceTriangleElements->getElement(k);
1762
1763 vecTmp[0] = points->at(elementTmp.getNode(0)).at(0) - points->at(elementTmp.getNode(1)).at(0);
1764 vecTmp[1] = points->at(elementTmp.getNode(0)).at(1) - points->at(elementTmp.getNode(1)).at(1);
1765 vecTmp[2] = points->at(elementTmp.getNode(0)).at(2) - points->at(elementTmp.getNode(1)).at(2);
1766 lengthA = std::sqrt(std::pow(vecTmp[0],2)+std::pow(vecTmp[1],2)+std::pow(vecTmp[2],2));
1767
1768 vecTmp1[0] = points->at(elementTmp.getNode(0)).at(0) - points->at(elementTmp.getNode(2)).at(0);
1769 vecTmp1[1] = points->at(elementTmp.getNode(0)).at(1) - points->at(elementTmp.getNode(2)).at(1);
1770 vecTmp1[2] = points->at(elementTmp.getNode(0)).at(2) - points->at(elementTmp.getNode(2)).at(2);
1771 lengthB = std::sqrt(std::pow(vecTmp1[0],2)+std::pow(vecTmp1[1],2)+std::pow(vecTmp1[2],2));
1772
1773 vecTmp2[0] = points->at(elementTmp.getNode(1)).at(0) - points->at(elementTmp.getNode(2)).at(0);
1774 vecTmp2[1] = points->at(elementTmp.getNode(1)).at(1) - points->at(elementTmp.getNode(2)).at(1);
1775 vecTmp2[2] = points->at(elementTmp.getNode(1)).at(2) - points->at(elementTmp.getNode(2)).at(2);
1776 lengthC = std::sqrt(std::pow(vecTmp2[0],2)+std::pow(vecTmp2[1],2)+std::pow(vecTmp2[2],2));
1777
1778 s1 = (lengthA+lengthB+lengthC)/2.;
1779 areaTriangles[k] = std::sqrt(s1*(s1-lengthA)*(s1-lengthB)*(s1-lengthC));
1780
1781 diamElements[k] = 2*(lengthA *lengthB *lengthC)/(4*areaTriangles[k]);
1782
1783 rho_T[k] = 4*(areaTriangles[k]) / (lengthA+lengthB+lengthC);
1784
1785 C_T[k] = diamElements[k] / rho_T[k];
1786 }
1787 return diamElements;
1788}
1789
1790
1801
1802template <class SC, class LO, class GO, class NO>
1803vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::determineAreaTriangles(ElementsPtr_Type elements,EdgeElementsPtr_Type edgeElements, SurfaceElementsPtr_Type surfaceElements, vec2D_dbl_ptr_Type points){
1804
1805 vec_dbl_Type areaTriangles(surfaceElements->numberElements());
1806 vec_dbl_Type vecTmp(3),vecTmp1(3),vecTmp2(3);
1807
1808 FiniteElement surfaceTmp;
1809 for(int k=0; k< surfaceElements->numberElements() ; k++){
1810
1811 double lengthA, lengthB, lengthC,s1;
1812
1813 surfaceTmp= surfaceElements->getElement(k);
1814 vecTmp[0] = points->at(surfaceTmp.getNode(0)).at(0) - points->at(surfaceTmp.getNode(1)).at(0);
1815 vecTmp[1] = points->at(surfaceTmp.getNode(0)).at(1) - points->at(surfaceTmp.getNode(1)).at(1);
1816 vecTmp[2] = points->at(surfaceTmp.getNode(0)).at(2) - points->at(surfaceTmp.getNode(1)).at(2);
1817 lengthA = std::sqrt(std::pow(vecTmp[0],2)+std::pow(vecTmp[1],2)+std::pow(vecTmp[2],2));
1818
1819
1820 vecTmp1[0] = points->at(surfaceTmp.getNode(0)).at(0) - points->at(surfaceTmp.getNode(2)).at(0);
1821 vecTmp1[1] = points->at(surfaceTmp.getNode(0)).at(1) - points->at(surfaceTmp.getNode(2)).at(1);
1822 vecTmp1[2] = points->at(surfaceTmp.getNode(0)).at(2) - points->at(surfaceTmp.getNode(2)).at(2);
1823 lengthB = std::sqrt(std::pow(vecTmp1[0],2)+std::pow(vecTmp1[1],2)+std::pow(vecTmp1[2],2));
1824
1825 vecTmp2[0] = points->at(surfaceTmp.getNode(1)).at(0) - points->at(surfaceTmp.getNode(2)).at(0);
1826 vecTmp2[1] = points->at(surfaceTmp.getNode(1)).at(1) - points->at(surfaceTmp.getNode(2)).at(1);
1827 vecTmp2[2] = points->at(surfaceTmp.getNode(1)).at(2) - points->at(surfaceTmp.getNode(2)).at(2);
1828 lengthC = std::sqrt(std::pow(vecTmp2[0],2)+std::pow(vecTmp2[1],2)+std::pow(vecTmp2[2],2));
1829
1830 s1 = (lengthA+lengthB+lengthC)/2.;
1831 areaTriangles[k] = std::sqrt(s1*(s1-lengthA)*(s1-lengthB)*(s1-lengthC));
1832 }
1833 return areaTriangles;
1834}
1835
1846
1847template <class SC, class LO, class GO, class NO>
1848vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>:: calcDiamTriangles(ElementsPtr_Type elements,vec2D_dbl_ptr_Type points, vec_dbl_Type& areaTriangles, vec_dbl_Type& rho_T, vec_dbl_Type& C_T){
1849
1850 vec_dbl_Type diamElements(elements->numberElements());
1851
1852
1853 FiniteElement elementTmp;
1854
1855 vec_dbl_Type vecTmp(2),vecTmp1(2),vecTmp2(2);
1856 for(int k=0; k< elements->numberElements() ; k++){
1857
1858
1859 double lengthA, lengthB, lengthC,s1;
1860 elementTmp = elements->getElement(k);
1861 vecTmp[0] = points->at(elementTmp.getNode(0)).at(0) - points->at(elementTmp.getNode(1)).at(0);
1862 vecTmp[1] = points->at(elementTmp.getNode(0)).at(1) - points->at(elementTmp.getNode(1)).at(1);
1863 lengthA = std::sqrt(std::pow(vecTmp[0],2)+std::pow(vecTmp[1],2));
1864
1865
1866 vecTmp1[0] = points->at(elementTmp.getNode(0)).at(0) - points->at(elementTmp.getNode(2)).at(0);
1867 vecTmp1[1] = points->at(elementTmp.getNode(0)).at(1) - points->at(elementTmp.getNode(2)).at(1);
1868 lengthB = std::sqrt(std::pow(vecTmp1[0],2)+std::pow(vecTmp1[1],2));
1869
1870 vecTmp2[0] = points->at(elementTmp.getNode(1)).at(0) - points->at(elementTmp.getNode(2)).at(0);
1871 vecTmp2[1] = points->at(elementTmp.getNode(1)).at(1) - points->at(elementTmp.getNode(2)).at(1);
1872 lengthC = std::sqrt(std::pow(vecTmp2[0],2)+std::pow(vecTmp2[1],2));
1873
1874 s1 = (lengthA+lengthB+lengthC)/2.;
1875 double area = std::sqrt(s1*(s1-lengthA)*(s1-lengthB)*(s1-lengthC));
1876
1877 areaTriangles[k] = area;
1878
1879 diamElements[k] = 2*(lengthA *lengthB *lengthC)/(4*area);
1880 rho_T[k] = 4*(area) / (lengthA+lengthB+lengthC);
1881
1882 C_T[k] = diamElements[k] / rho_T[k];
1883
1884 }
1885 return diamElements;
1886}
1887
1888
1889
1890
1891
1906template <class SC, class LO, class GO, class NO>
1907typename ErrorEstimation<SC,LO,GO,NO>::MultiVectorPtr_Type ErrorEstimation<SC,LO,GO,NO>::determineCoarseningError(MeshUnstrPtr_Type mesh_k, MeshUnstrPtr_Type mesh_k_m, MultiVectorPtr_Type errorElementMv_k, std::string distribution, std::string markingStrategy, double theta){
1908
1909 theta_ =theta;
1910 markingStrategy_ = markingStrategy;
1911
1912 // We determine which meshes we need to focus on.
1913 // Mesh of level k ist mesh_k and the mesh at the beginning of 'iteration'
1914 ElementsPtr_Type elements_k = mesh_k->getElementsC();
1915 // Mesh of level k-m is then mesh_k_m, an the mesh that helps us determine the new coarsening error
1916 ElementsPtr_Type elements_k_m = mesh_k_m->getElementsC();
1917
1918 MultiVectorPtr_Type errorElementMv_k_m = Teuchos::rcp(new MultiVector_Type( mesh_k_m->getElementMap()) );
1919 Teuchos::ArrayRCP<SC> errorElement_k_m = errorElementMv_k_m->getDataNonConst(0);
1920
1921 Teuchos::ArrayRCP<SC> errorElement_k = errorElementMv_k->getDataNonConst(0);
1922
1923 // We determine the error of mesh level k-m with the error of mesh level k
1924 if(distribution == "backwards"){
1925 for(int i=0; i< errorElement_k.size(); i++){
1926 errorElement_k_m[elements_k->getElement(i).getPredecessorElement()] += errorElement_k[i];
1927 }
1928 }
1929 if(distribution == "forwards"){
1930 for(int i=0; i< errorElement_k_m.size(); i++){
1931 errorElement_k_m[i] = errorElement_k[elements_k_m->getElement(i).getPredecessorElement()];
1932 }
1933 }
1934 return errorElementMv_k_m;
1935
1936 // As a result we coarsened the mesh at level k= iteration with the factor 'm'
1937
1938
1939}
1940
1948
1949template <class SC, class LO, class GO, class NO>
1950void ErrorEstimation<SC,LO,GO,NO>::updateElementsOfSurfaceLocalAndGlobal(EdgeElementsPtr_Type edgeElements, SurfaceElementsPtr_Type surfaceTriangleElements){
1951
1952 // The vector contains the Information which elements contain the edge
1953 vec2D_GO_Type elementsOfEdgeGlobal = edgeElements->getElementsOfEdgeGlobal();
1954 vec2D_LO_Type elementsOfEdgeLocal = edgeElements->getElementsOfEdgeLocal();
1955
1956 vec2D_GO_Type elementsOfSurfaceGlobal = surfaceTriangleElements->getElementsOfSurfaceGlobal();
1957
1958 GO nEl;
1959 vec_LO_Type edgeTmp1, edgeTmp2, edgeTmp3;
1960 LO id1, id2, id3;
1961
1962 vec2D_LO_Type edgeList(0,vec_LO_Type(2));
1963
1964 for(int i=0; i<edgeElements->numberElements(); i++){
1965 edgeList.push_back({edgeElements->getElement(i).getNode(0),edgeElements->getElement(i).getNode(1)});
1966 }
1967
1968 for(int i=0; i < surfaceTriangleElements->numberElements() ; i++){
1969 FiniteElement surfaceTmp = surfaceTriangleElements->getElement(i);
1970 if(surfaceTmp.isInterfaceElement() && (surfaceTmp.getFlag() == 0 || surfaceTmp.getFlag() == 10)){
1971
1972 vec_int_Type tmpElements(0);
1973 //this->surfaceTriangleElements_->setElementsOfSurfaceLocalEntry(i,-1);
1974
1975 edgeTmp1={surfaceTmp.getNode(0),surfaceTmp.getNode(1)};
1976 edgeTmp2={surfaceTmp.getNode(1),surfaceTmp.getNode(2)};
1977 //edgeTmp3={surfaceTmp.getNode(1),surfaceTmp.getNode(2)};
1978
1979 sort(edgeTmp1.begin(),edgeTmp1.end());
1980 sort(edgeTmp2.begin(),edgeTmp2.end());
1981 //sort(edgeTmp3.begin(),edgeTmp3.end());
1982
1983 auto it1 = find( edgeList.begin(), edgeList.end() ,edgeTmp1 );
1984 id1 = distance( edgeList.begin() , it1 );
1985
1986 auto it2 = find(edgeList.begin(), edgeList.end() ,edgeTmp2 );
1987 id2 = distance(edgeList.begin() , it2 );
1988
1989 for(int j=0 ; j< elementsOfEdgeGlobal[id1].size() ; j++)
1990 tmpElements.push_back( elementsOfEdgeGlobal[id1][j]);
1991
1992 for(int j=0 ; j< elementsOfEdgeGlobal[id2].size() ; j++)
1993 tmpElements.push_back( elementsOfEdgeGlobal[id2][j]);
1994
1995
1996 sort(tmpElements.begin(),tmpElements.end());
1997 bool found =false;
1998
1999 //cout << "Surface Element da " << elementsOfSurfaceGlobal[i][0] << " tmp elements in question " << tmpElements[0] << " " ;
2000 for(int j=0; j< tmpElements.size()-1; j++){
2001 if((tmpElements[j] == tmpElements[j+1] )&& (tmpElements[j] != elementsOfSurfaceGlobal[i][0]) && (found==false)) {
2002 nEl = tmpElements[j];
2003 surfaceTriangleElements->setElementsOfSurfaceGlobalEntry(i,nEl);
2004 found = true;
2005 }
2006 }
2007 if(found == false){
2008 cout << " No Element Found for edges " << id1 << " " << id2 << " on Proc " << inputMesh_->getComm()->getRank() << " und element schon da=" << elementsOfSurfaceGlobal[i][0] << endl;
2009 cout << " Tmp1= ";
2010 for(int j=0; j< tmpElements.size(); j++){
2011 cout << " " << tmpElements[j];
2012 }
2013 cout<< " Flag Surface " << surfaceTmp.getFlag() << endl;
2014 }
2015 }
2016 }
2017
2018 vec2D_LO_Type elementsOfSurfaceLocal = surfaceTriangleElements->getElementsOfSurfaceLocal();
2019 elementsOfSurfaceGlobal = surfaceTriangleElements->getElementsOfSurfaceGlobal();
2020}
2021
2030
2031template <class SC, class LO, class GO, class NO>
2032void ErrorEstimation<SC,LO,GO,NO>::buildTriangleMap(){
2033
2034 int maxRank = std::get<1>(inputMesh_->rankRange_);
2035 const int myRank = inputMesh_->getComm()->getRank();
2036
2037 vec_GO_Type globalProcs(0);
2038 for (int i=0; i<= maxRank; i++)
2039 globalProcs.push_back(i);
2040
2041 Teuchos::ArrayView<GO> globalProcArray = Teuchos::arrayViewFromVector( globalProcs);
2042
2043 vec_GO_Type localProc(0);
2044 localProc.push_back(inputMesh_->getComm()->getRank());
2045 Teuchos::ArrayView<GO> localProcArray = Teuchos::arrayViewFromVector( localProc);
2046
2047 MapPtr_Type mapGlobalProc =
2048 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), globalProcArray, 0, inputMesh_->getComm()) );
2049
2050 MapPtr_Type mapProc =
2051 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), localProcArray, 0, inputMesh_->getComm()) );
2052
2053
2054
2055 vec2D_int_Type interfaceSurfacesLocalId(1,vec_int_Type(1));
2056
2057
2058 MultiVectorLOPtr_Type exportLocalEntry = Teuchos::rcp( new MultiVectorLO_Type( mapProc, 1 ) );
2059
2060 // (A) First we determine a Map only for the interface Nodes
2061 // This will reduce the size of the Matrix we build later significantly if only look at the interface edges
2062 int numSurfaces= surfaceElements_->numberElements();
2063 vec2D_GO_Type inzidenzIndices(0,vec_GO_Type(2)); // Vector that stores global IDs of each edge (in Repeated Sense)
2064 vec_LO_Type localSurfaceIndex(0); // stores the local ID of surfaces in question
2065 vec_GO_Type id(2);
2066 int surfacesUnique=0;
2067
2068 vec2D_dbl_ptr_Type points = inputMesh_->getPointsRepeated();
2069
2070 vec2D_GO_Type elementsOfSurfaceGlobal = surfaceElements_->getElementsOfSurfaceGlobal();
2071
2072 vec_GO_Type elementRep(0);
2073
2074 int interfaceNum=0;
2075 for(int i=0; i<numSurfaces; i++ ){
2076 if(surfaceElements_->getElement(i).isInterfaceElement()){
2077
2078 id[0] = elementsOfSurfaceGlobal[i][0];
2079 id[1] = elementsOfSurfaceGlobal[i][1];
2080
2081
2082
2083 sort(id.begin(),id.end());
2084 inzidenzIndices.push_back(id);
2085
2086 localSurfaceIndex.push_back(i);
2087 interfaceNum++;
2088
2089 }
2090
2091 else{
2092 surfacesUnique++;
2093 }
2094
2095 for(int j=0; j < elementsOfSurfaceGlobal[i].size() ; j++)
2096 elementRep.push_back(elementsOfSurfaceGlobal[i][j]);
2097
2098 }
2099
2100 sort(elementRep.begin(),elementRep.end());
2101 vec_GO_Type::iterator ip = unique( elementRep.begin(), elementRep.end());
2102 elementRep.resize(distance(elementRep.begin(), ip));
2103
2104 Teuchos::ArrayView<GO> elementRepArray = Teuchos::arrayViewFromVector( elementRep);
2105
2106 MapPtr_Type elementMapRep =
2107 Teuchos::rcp( new Map_Type( Teuchos::OrdinalTraits<GO>::invalid(), elementRepArray, 0, inputMesh_->getComm()) );
2108
2109
2110
2111 // This Matrix is row based, where the row is based on mapInterfaceNodesUnqiue
2112 // We then add a '1' Entry when two global Node IDs form an edge
2113 MatrixPtr_Type inzidenzMatrix = Teuchos::rcp( new Matrix_Type(inputMesh_->getElementMap(), 40 ) );
2114 Teuchos::Array<GO> index(1);
2115 Teuchos::Array<GO> col(1);
2116 Teuchos::Array<SC> value(1, Teuchos::ScalarTraits<SC>::one() );
2117
2118 for(int i=0; i<inzidenzIndices.size(); i++ ){
2119 index[0] = inzidenzIndices[i][0];
2120 col[0] = inzidenzIndices[i][1];
2121 inzidenzMatrix->insertGlobalValues(index[0], col(), value());
2122 }
2123 inzidenzMatrix->fillComplete();
2124
2125
2126 // ---------------------------------------------------
2127 // 2. Set unique edges IDs ---------------------------
2128 // Setting the IDs of Edges that are uniquely on one
2129 // Processor
2130 // ---------------------------------------------------
2131 exportLocalEntry->putScalar( (LO) surfacesUnique );
2132
2133 MultiVectorLOPtr_Type newSurfacesUniqueGlobal= Teuchos::rcp( new MultiVectorLO_Type( mapGlobalProc, 1 ) );
2134 newSurfacesUniqueGlobal->putScalar( (LO) 0 );
2135 newSurfacesUniqueGlobal->importFromVector( exportLocalEntry, false, "Insert");
2136
2137 // offset EdgesUnique for proc and globally
2138 Teuchos::ArrayRCP< const LO > newSurfacesList = newSurfacesUniqueGlobal->getData(0);
2139
2140 GO procOffsetSurface=0;
2141 for(int i=0; i< myRank; i++)
2142 procOffsetSurface= procOffsetSurface + newSurfacesList[i];
2143
2144 // global IDs for map
2145 vec_GO_Type vecGlobalIDsSurfaces(numSurfaces);
2146
2147 // Step 1: adding unique global edge IDs
2148 int count=0;
2149 for(int i=0; i< numSurfaces; i++){
2150 if(!surfaceElements_->getElement(i).isInterfaceElement()){
2151 vecGlobalIDsSurfaces.at(i) = procOffsetSurface+count;
2152 count++;
2153 }
2154 }
2155
2156 // Now we add the repeated ids, by first turning interfaceEdgesTag into a map
2157 // Offset for interface IDS:
2158 GO offsetInterface=0;
2159 for(int i=0; i< maxRank+1; i++)
2160 offsetInterface= offsetInterface + newSurfacesList[i];
2161
2162 // (D) Now we count the row entries on each processor an set global IDs
2163
2164 Teuchos::ArrayView<const LO> indices;
2165 Teuchos::ArrayView<const SC> values;
2166 vec2D_GO_Type inzidenzIndicesUnique(0,vec_GO_Type(2)); // Vector that stores only both global IDs if the first is part of my unique Interface Nodes
2167 MapConstPtr_Type colMap = inzidenzMatrix->getMap("col");
2168 MapConstPtr_Type rowMap = inzidenzMatrix->getMap("row");
2169 int numRows = rowMap->getNodeNumElements();
2170 int uniqueSurfaces =0;
2171 for(int i=0; i<numRows; i++ ){
2172 inzidenzMatrix->getLocalRowView(i, indices,values);
2173 uniqueSurfaces = uniqueSurfaces+indices.size();
2174 vec_GO_Type surfaceTmp(2);
2175 for(int j=0; j<indices.size(); j++){
2176 surfaceTmp[0] = rowMap->getGlobalElement(i);
2177 surfaceTmp[1] = colMap->getGlobalElement(indices[j]);
2178 inzidenzIndicesUnique.push_back(surfaceTmp);
2179 }
2180 }
2181
2182 exportLocalEntry->putScalar( uniqueSurfaces);
2183 MultiVectorLOPtr_Type newSurfaceInterfaceGlobal= Teuchos::rcp( new MultiVectorLO_Type( mapGlobalProc, 1 ) );
2184 newSurfaceInterfaceGlobal->putScalar( (LO) 0 );
2185 newSurfaceInterfaceGlobal->importFromVector( exportLocalEntry,true, "Insert");
2186
2187 // offset EdgesUnique for proc and globally
2188 Teuchos::ArrayRCP< const LO > numUniqueInterface = newSurfaceInterfaceGlobal->getData(0);
2189
2190 procOffsetSurface=0;
2191 for(int i=0; i< myRank; i++)
2192 procOffsetSurface= procOffsetSurface + numUniqueInterface[i];
2193
2194 int numInterfaceSurface=0;
2195
2196 vec_GO_Type uniqueInterfaceIDsList_(inzidenzIndicesUnique.size());
2197 for(int i=0; i< uniqueInterfaceIDsList_.size(); i++)
2198 uniqueInterfaceIDsList_[i] = procOffsetSurface + i;
2199
2200 MatrixPtr_Type indMatrix = Teuchos::rcp( new Matrix_Type(inputMesh_->getElementMap(), 40 ) );
2201
2202 for(int i=0; i<inzidenzIndicesUnique.size(); i++ ){
2203 index[0] = inzidenzIndicesUnique[i][0];
2204 col[0] = inzidenzIndicesUnique[i][1];
2205 Teuchos::Array<SC> value2(1,uniqueInterfaceIDsList_[i]);
2206 indMatrix->insertGlobalValues(index[0], col(), value2());
2207 }
2208 indMatrix->fillComplete();
2209
2210 MatrixPtr_Type importMatrix = Teuchos::rcp( new Matrix_Type(elementMapRep, 40 ) );
2211
2212 importMatrix->importFromVector(indMatrix,false,"Insert");
2213 importMatrix->fillComplete();
2214
2215 // Determine global indices
2216 GO surfaceID=0;
2217 colMap = importMatrix->getMap("col");
2218 rowMap = importMatrix->getMap("row");
2219 LO valueID=0;
2220 bool found = false;
2221 GO entry =0;
2222 for(int i=0; i<inzidenzIndices.size(); i++ ){
2223
2224 importMatrix->getLocalRowView(rowMap->getLocalElement(inzidenzIndices[i][0]), indices,values); // Indices and values connected to node i / row i in Matrix
2225 // Entries in 'indices' represent the local entry in 'colmap
2226 // with 'getGlobalElement' we know the global Node ID that belongs to the first Node that form an edge
2227 // vector in with entries only for edges belonging to node i;
2228 vec2D_GO_Type indicesTmp(indices.size(),vec_GO_Type(2));
2229 vec_GO_Type indTmp(2);
2230 for(int j=0; j<indices.size(); j++){
2231 indTmp[0] = colMap->getGlobalElement(indices[j]);
2232 indTmp[1] = values[j];
2233 indicesTmp.push_back(indTmp); // vector with the indices and values belonging to node i
2234 }
2235
2236 found = false;
2237 for(int k=0; k<indicesTmp.size();k++){
2238 if(inzidenzIndices[i][1] == indicesTmp[k][0]){
2239 entry =k;
2240 k = indicesTmp.size();
2241 surfaceID = indicesTmp[entry][1];
2242 vecGlobalIDsSurfaces.at(localSurfaceIndex[i]) = offsetInterface + surfaceID;
2243 found =true;
2244 }
2245 }
2246 if(found == false)
2247 cout << " Asking for row " << rowMap->getLocalElement(inzidenzIndices[i][0]) << " for Edge [" << inzidenzIndices[i][0] << ", " << inzidenzIndices[i][1] << "], on Proc " << myRank << " but no Value found " <<endl;
2248 }
2249
2250
2251 Teuchos::RCP<std::vector<GO>> surfacesGlobMapping = Teuchos::rcp( new std::vector<GO>( vecGlobalIDsSurfaces ) );
2252 Teuchos::ArrayView<GO> surfacesGlobMappingArray = Teuchos::arrayViewFromVector( *surfacesGlobMapping);
2253
2254 this->surfaceTriangleMap_.reset(new Map<LO,GO,NO>(Teuchos::OrdinalTraits<GO>::invalid(), surfacesGlobMappingArray, 0, inputMesh_->getComm()) );
2255}
2256
2266
2267template <class SC, class LO, class GO, class NO>
2268vec2D_dbl_Type ErrorEstimation<SC,LO,GO,NO>::gradPhi(int dim,
2269 int intFE,
2270 vec_dbl_Type &p){
2271
2272 int numNodes = dim+1;
2273 if(intFE == 2){
2274 numNodes=6;
2275 if(dim==3)
2276 numNodes=10;
2277 }
2278
2279 vec2D_dbl_Type value(numNodes,vec_dbl_Type(dim));
2280
2281 if (dim==2) {
2282 switch (intFE) {
2283 case 1://P1
2284 value[0][0]= -1.;
2285 value[0][1]= -1.;
2286
2287 value[1][0]= 1.;
2288 value[1][1]= 0.;
2289
2290 value[2][0]= 0.;
2291 value[2][1]= 1.;
2292
2293 break;
2294 case 2://P2
2295 value[0][0]= 1. - 4.*(1 - p[0] - p[1]);
2296 value[0][1]= 1. - 4.*(1 - p[0] - p[1]);
2297
2298 value[1][0]= 4.*p[0] - 1;
2299 value[1][1]= 0.;
2300
2301 value[2][0]= 0.;
2302 value[2][1]= 4.*p[1] - 1;
2303
2304 value[3][0]= 4 * (1. - 2*p[0] - p[1]);
2305 value[3][1]= -4 * p[0];
2306
2307 value[4][0]= 4.*p[1];
2308 value[4][1]= 4.*p[0];
2309
2310 value[5][0]= - 4.*p[1];
2311 value[5][1]= 4 * (1. - p[0] - 2*p[1]);
2312
2313 break;
2314 }
2315 }
2316 else if(dim==3) {
2317 switch (intFE) {
2318 case 1://P1
2319
2320 value[0][0]= -1.;
2321 value[0][1]= -1.;
2322 value[0][2]= -1.;
2323
2324 value[1][0]= 1.;
2325 value[1][1]= 0.;
2326 value[1][2]= 0.;
2327
2328 value[2][0]= 0.;
2329 value[2][1]= 1.;
2330 value[2][2]= 0.;
2331
2332 value[3][0]= 0.;
2333 value[3][1]= 0.;
2334 value[3][2]= 1.;
2335
2336 break;
2337
2338 case 2://P2
2339
2340 value[0][0]= -3. + 4.*p[0] + 4.*p[1] + 4.*p[2];
2341 value[0][1]= -3. + 4.*p[0] + 4.*p[1] + 4.*p[2];
2342 value[0][2]= -3. + 4.*p[0] + 4.*p[1] + 4.*p[2];
2343
2344 value[1][0]= 4.*p[0] - 1;
2345 value[1][1]= 0.;
2346 value[1][2]= 0.;
2347
2348 value[2][0]= 0.;
2349 value[2][1]= 4.*p[1] - 1;
2350 value[2][2]= 0.;
2351
2352 value[3][0]= 0.;
2353 value[3][1]= 0.;
2354 value[3][2]= 4.*p[2] - 1;
2355
2356 value[4][0]= 4. - 8.*p[0] - 4.*p[1] - 4.*p[2];
2357 value[4][1]= - 4.*p[0];
2358 value[4][2]= - 4.*p[0];
2359
2360 value[5][0]= 4.*p[1];
2361 value[5][1]= 4.*p[0];
2362 value[5][2]= 0.;
2363
2364 value[6][0]= - 4.*p[1];
2365 value[6][1]= 4. - 4.*p[0] - 8.*p[1] - 4.*p[2];
2366 value[6][2]= - 4.*p[1];
2367
2368 value[7][0]= - 4.*p[2];
2369 value[7][1]= - 4.*p[2];
2370 value[7][2]= 4. - 4.*p[0] - 4.*p[1] - 8.*p[2];
2371
2372 value[8][0]= 4.*p[2];
2373 value[8][1]= 0.;
2374 value[8][2]= 4.*p[0];
2375
2376 value[9][0]= 0.;
2377 value[9][1]= 4.*p[2];
2378 value[9][2]= 4.*p[1];
2379
2380
2381 break;
2382 }
2383 }
2384 return value;
2385
2386}
2387
2397
2398template <class SC, class LO, class GO, class NO>
2399vec_dbl_Type ErrorEstimation<SC,LO,GO,NO>::phi(int dim,
2400 int intFE,
2401 vec_dbl_Type &p){
2402
2403 int numNodes = dim+1;
2404 if(intFE == 2){
2405 numNodes=6;
2406 if(dim==3)
2407 numNodes=10;
2408 }
2409
2410 vec_dbl_Type value(numNodes);
2411
2412 if (dim==2) {
2413 switch (intFE) {
2414 case 1://P1
2415 value[0]= p[0];
2416 value[1]= p[1];
2417 value[2]= 1-p[0]-p[1];
2418
2419 break;
2420
2421 case 2://P2
2422
2423 value[0]= -(1. - p[0]-p[1]) * (1 - 2.*(1-p[0] - p[1]));
2424 value[1] = -p[0] * (1 - 2*p[0]);
2425 value[2] = -p[1] * (1 - 2*p[1]);
2426 value[3] = 4*p[0] * (1 - p[0]-p[1]);
2427 value[4] = 4*p[0]*p[1];
2428 value[5] = 4*p[1] * (1 - p[0]-p[1]);
2429 break;
2430 }
2431 }
2432
2433 else if(dim==3){
2434 switch (intFE) {
2435 case 1://P1
2436
2437 value[0] = (1. - p[0]-p[1]-p[2]);
2438 value[1] = p[0];
2439 value[2] = p[1];
2440 value[3] = p[2];
2441
2442 break;
2443 case 2: //P2
2444 value[0] = (1. - p[0]-p[1]-p[2]) * (1 - 2*p[0] - 2*p[1] - 2*p[2]);
2445 value[1] = p[0] * (2*p[0] - 1);
2446 value[2] = p[1] * (2*p[1] - 1);
2447 value[3] = p[2] * (2*p[2] - 1);
2448 value[4] = 4*p[0] * (1 - p[0]-p[1]-p[2]);
2449 value[5] = 4*p[0]*p[1];
2450 value[6] = 4*p[1] * (1 - p[0]-p[1]-p[2]);
2451 value[7] = 4*p[2] * (1 - p[0]-p[1]-p[2]);
2452 value[8] = 4*p[0]*p[2];
2453 value[9] = 4*p[1]*p[2];
2454 break;
2455 }
2456 }
2457 return value;
2458
2459}
2460
2461
2462}
2463#endif
2464
FEDD::ErrorEstimation::determineAreaTriangles
vec_dbl_Type determineAreaTriangles(ElementsPtr_Type elements, EdgeElementsPtr_Type edgeElements, SurfaceElementsPtr_Type surfaceElements, vec2D_dbl_ptr_Type points)
Calculating the area of the triangle elements of tetrahedra.
Definition ErrorEstimation_def.hpp:1803
FEDD::ErrorEstimation::getQuadValues
vec2D_dbl_Type getQuadValues(int dim, std::string FEType, std::string Type, vec_dbl_Type &QuadW, FiniteElement surface)
Returns neccesary quadrature Values. Is distinguishes between needing Element or Surface information.
Definition ErrorEstimation_def.hpp:1485
FEDD::ErrorEstimation::calcDiamTetraeder
vec_dbl_Type calcDiamTetraeder(ElementsPtr_Type elements, vec2D_dbl_ptr_Type points, vec_dbl_Type volTet)
Calculating the circumdiameter of tetraeder.
Definition ErrorEstimation_def.hpp:1677
FEDD::ErrorEstimation::determineDivU
double determineDivU(FiniteElement element)
Function that that determines || div(u) ||_T for a Element T.
Definition ErrorEstimation_def.hpp:1226
FEDD::ErrorEstimation::calculateJump
vec_dbl_Type calculateJump()
Part of the error estimator that calculates the jump part of the estimation. What kind of jump is cal...
Definition ErrorEstimation_def.hpp:662
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::phi
vec_dbl_Type phi(int dim, int intFE, vec_dbl_Type &p)
Calcutlating phi depending on quad points p.
Definition ErrorEstimation_def.hpp:2399
FEDD::ErrorEstimation::identifyProblem
void identifyProblem(BlockMultiVectorConstPtr_Type valuesSolution)
Identifying the problem with respect to the degrees of freedom and whether we calculate pressure....
Definition ErrorEstimation_def.hpp:60
FEDD::ErrorEstimation::calcDiamTriangles3D
vec_dbl_Type calcDiamTriangles3D(SurfaceElementsPtr_Type surfaceTriangleElements, vec2D_dbl_ptr_Type points, vec_dbl_Type &areaTriangles, vec_dbl_Type &rho_T, vec_dbl_Type &C_T)
Calculating the diameter of triangles.
Definition ErrorEstimation_def.hpp:1751
FEDD::ErrorEstimation::gradPhi
vec2D_dbl_Type gradPhi(int dim, int intFE, vec_dbl_Type &p)
Calcutlating the gradient of phi depending on quad points p.
Definition ErrorEstimation_def.hpp:2268
FEDD::ErrorEstimation::buildTriangleMap
void buildTriangleMap()
Build Surface Map. Contrary to building the edge map, building the surface map is somewhat simpler as...
Definition ErrorEstimation_def.hpp:2032
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::calcRhoTetraeder
vec_dbl_Type calcRhoTetraeder(ElementsPtr_Type elements, SurfaceElementsPtr_Type surfaceTriangleElements, vec_dbl_Type volTet, vec_dbl_Type areaTriangles)
Calcutlating the incircumdiameter of tetrahedra.
Definition ErrorEstimation_def.hpp:1723
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::calcNPhi
vec3D_dbl_Type calcNPhi(std::string phiDerivative, int dofsSol, std::string FEType)
Function that calculates the jump part for nabla u or p.
Definition ErrorEstimation_def.hpp:782
FEDD::ErrorEstimation::updateElementsOfSurfaceLocalAndGlobal
void updateElementsOfSurfaceLocalAndGlobal(EdgeElementsPtr_Type edgeElements, SurfaceElementsPtr_Type surfaceTriangleElements)
UpdateElementsOfSurfaceLocalAndGlobal is performed here instead of in meshRefinement,...
Definition ErrorEstimation_def.hpp:1950
FEDD::ErrorEstimation::determineResElement
double determineResElement(FiniteElement element, RhsFunc_Type rhsFunc)
Function that that determines ||\Delta u_h + f ||_(L2(T)), || \Delta u_h + f - \nabla p_h ||_T or || ...
Definition ErrorEstimation_def.hpp:1319
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::calcDiamTriangles
vec_dbl_Type calcDiamTriangles(ElementsPtr_Type elements, vec2D_dbl_ptr_Type points, vec_dbl_Type &areaTriangles, vec_dbl_Type &rho_T, vec_dbl_Type &C_T)
Calculating the diameter of triangles.
Definition ErrorEstimation_def.hpp:1848
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::ErrorEstimation::determineVolTet
vec_dbl_Type determineVolTet(ElementsPtr_Type elements, vec2D_dbl_ptr_Type points)
function, that determines volume of tetrahedra.
Definition ErrorEstimation_def.hpp:1622
FEDD::ErrorEstimation::makeRepeatedSolution
void makeRepeatedSolution(BlockMultiVectorConstPtr_Type valuesSolution)
We split the solution from the BlockMultiVector valuesSolution into one or more seperate blocks,...
Definition ErrorEstimation_def.hpp:98
FEDD::FiniteElement
Definition FiniteElement.hpp:17
FEDD::Map
Definition Map_decl.hpp:30
FEDD::SmallMatrix
This class represents a templated small Matrix of type T. Primarily created for 2x2 and 3x3 matrices....
Definition SmallMatrix.hpp:20
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36