Commit 44e00260 authored by Praetorius, Simon's avatar Praetorius, Simon

Merge branch 'issue/evalAtPoint' into 'dev'

correct parallel evalAtPoint functions if point is found on multiple partitions

See merge request !27
parents 0269a0c9 ea26c5eb
......@@ -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;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment