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