From ea26c5eb1f8c3692f62e94dd34d80bc94979694d Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 22 Jan 2019 13:55:56 +0100 Subject: [PATCH] correct parallel evalAtPoint functions if point is found on multiple partitions --- AMDiS/src/DOFVector.hh | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 61b1c767..21ad9743 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -124,10 +124,10 @@ namespace AMDiS { { this->name = n; this->feSpace = f; - + if (this->feSpace && this->feSpace->getAdmin()) (this->feSpace->getAdmin())->addDOFIndexed(this); - + this->boundaryManager = new BoundaryManager(f); #ifdef HAVE_PARALLEL_DOMAIN_AMDIS if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL) @@ -142,7 +142,7 @@ namespace AMDiS { if ( Parallel::MeshDistributor::globalMeshDistributor != NULL) Parallel::MeshDistributor::globalMeshDistributor->removeInterchangeVector(this); #endif - + #ifdef _OPENMP #pragma omp critical #endif @@ -1090,7 +1090,7 @@ namespace AMDiS { ElInfo *elInfo = mesh->createNewElInfo(); double value = 0.0; - bool inside = false; + unsigned short inside = 0; if (oldElInfo && oldElInfo->getMacroElement()) { inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL); @@ -1101,28 +1101,24 @@ namespace AMDiS { if (oldElInfo) oldElInfo = elInfo; - if (inside) { + if (inside > 0) { basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices); ElementVector uh(nBasFcts); for (int i = 0; i < nBasFcts; i++) uh[i] = operator[](localIndices[i]); value = basFcts->evalUh(lambda, uh); - } else { -#ifdef HAVE_PARALLEL_DOMAIN_AMDIS - value = 0.0; -#else - ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry."); -#endif } - if (oldElInfo == NULL) delete elInfo; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS Parallel::mpi::globalAdd(value); + Parallel::mpi::globalAdd(inside); #endif - return value; + + TEST_EXIT(inside > 0)("Can not eval DOFVector at point p, because point is outside geometry."); + return value/inside; } @@ -1140,7 +1136,7 @@ namespace AMDiS { DimVec lambda(dim, NO_INIT); ElInfo *elInfo = mesh->createNewElInfo(); WorldVector value(DEFAULT_VALUE, 0.0); - bool inside = false; + unsigned short inside = 0; if (oldElInfo && oldElInfo->getMacroElement()) { inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL); @@ -1151,19 +1147,24 @@ namespace AMDiS { if (oldElInfo) oldElInfo = elInfo; - if (inside) { + if (inside > 0) { basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices); mtl::dense_vector > uh(nBasFcts); for (int i = 0; i < nBasFcts; i++) uh[i] = operator[](localIndices[i]); value = basFcts->evalUh(lambda, uh); - } else { - ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry."); } if (oldElInfo == NULL) delete elInfo; +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + Parallel::mpi::globalAdd(value); + Parallel::mpi::globalAdd(inside); +#endif + + TEST_EXIT(inside > 0)("Can not eval DOFVector at point p, because point is outside geometry."); + value *= 1.0/inside; return value; } -- GitLab