From ea26c5eb1f8c3692f62e94dd34d80bc94979694d Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
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<double> lambda(dim, NO_INIT);
     ElInfo *elInfo = mesh->createNewElInfo();
     WorldVector<double> 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<WorldVector<double> > 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