diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 8e0eb75e4a11b7927a90c2aea0d5d749a0ede12f..9850d8f2cd8a8824cb9b0c03d4e7b19b27861b8b 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -34,8 +34,6 @@ namespace AMDiS {
using namespace mtl;
- DOFMatrix *DOFMatrix::traversePtr = NULL;
-
DOFMatrix::DOFMatrix()
: rowFeSpace(NULL),
colFeSpace(NULL),
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 7a687267ffb386d39ab79b4094a9ac13b0fe4e14..2f3b082095271be758a45476525800205b9e9eb6 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -424,9 +424,6 @@ namespace AMDiS {
/// default compressed2D<double>
base_matrix_type matrix;
- /// Used while mesh traversal
- static DOFMatrix *traversePtr;
-
/// Pointers to all operators of the equation systems. Are used in the
/// assembling process.
vector<Operator*> operators;
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 601e4ede5df8ad2d35456a9068b8b2f2ed18461e..57212a3dcda8ba01296b7360aa52d952f98fa97a 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -68,7 +68,8 @@ namespace AMDiS {
template<>
- const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
+ double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
+ ElInfo *oldElInfo) const
{
FUNCNAME("DOFVector<double>::evalAtCoords()");
@@ -82,7 +83,7 @@ namespace AMDiS {
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
- static double value = 0.0;
+ double value = 0.0;
bool inside = false;
if (oldElInfo && oldElInfo->getMacroElement()) {
@@ -107,14 +108,13 @@ namespace AMDiS {
if (oldElInfo == NULL)
delete elInfo;
- if (values != NULL)
- *values = value;
return value;
- };
+ }
template<>
- const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
+ WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
+ ElInfo *oldElInfo) const
{
FUNCNAME("DOFVector<double>::evalAtCoords()");
@@ -129,9 +129,7 @@ namespace AMDiS {
ElInfo *elInfo = mesh->createNewElInfo();
- static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
- WorldVector<double> *val = (NULL != values) ? values : &Values;
-
+ WorldVector<double> value(DEFAULT_VALUE, 0.0);
bool inside = false;
if (oldElInfo && oldElInfo->getMacroElement()) {
@@ -148,7 +146,7 @@ namespace AMDiS {
mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
for (int i = 0; i < nBasFcts; i++)
uh[i] = operator[](localIndices[i]);
- *val = basFcts->evalUh(lambda, uh);
+ value = basFcts->evalUh(lambda, uh);
} else
throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
@@ -156,8 +154,8 @@ namespace AMDiS {
if (oldElInfo == NULL)
delete elInfo;
- return ((*val));
- };
+ return value;
+ }
template<>
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 16a337270809f69608cc9fa1187aeb93b08aee27..a8335726f0fcf59105a90ad9e83be9803016c9e9 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -613,13 +613,11 @@ namespace AMDiS {
/// Eval DOFVector at given point p. If oldElInfo != NULL the search for
/// the element, where p is inside, starts from oldElInfo. implemented for:
/// double, WorldVector< double >
- inline const T& evalAtPoint(WorldVector<double> &p,
- ElInfo *oldElInfo = NULL,
- T* value = NULL) const
+ T evalAtPoint(WorldVector<double> &p,
+ ElInfo *oldElInfo = NULL) const
{
FUNCNAME("DOFVector::evalAtPoint())");
TEST_EXIT(false)("Please implement your evaluation\n");
- return *value;
}
/// Determine the DegreeOfFreedom that has coords with minimal euclidean
@@ -674,12 +672,12 @@ namespace AMDiS {
BoundaryType boundaryType, Quadrature* q) const;
template<>
- const double& DOFVector<double>::evalAtPoint(
- WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
+ double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
+ ElInfo *oldElInfo) const;
template<>
- const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
- WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
+ WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
+ ElInfo *oldElInfo) const;
template<>
void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 54473490792bef3a5b58443e42ed241e20cec518..c54c1ff938495a7c7558158257b129cb389516dc 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -1599,9 +1599,9 @@ namespace AMDiS {
elObjDb.create(partitionMap, levelData);
elObjDb.updateRankData();
- unsigned long memsize = elObjDb.calculateMemoryUsage();
- MSG("Memory usage of element object database = %5.f KByte\n",
- static_cast<double>(memsize / 1024));
+ // unsigned long memsize = elObjDb.calculateMemoryUsage();
+ // MSG("Memory usage of element object database = %5.f KByte\n",
+ // static_cast<double>(memsize / 1024));
intBoundary.create(levelData, 0, elObjDb);
ParallelDebug::printBoundaryInfo(intBoundary);
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 7a8200c1c76c6858d94a0ff748adfc2885a9368d..34cb46b900302732af9063fb10b2844da2f99938 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -680,7 +680,7 @@ namespace AMDiS {
VecAssemblyBegin(tmpVec);
VecAssemblyEnd(tmpVec);
-
+
subdomain->solve(tmpVec, tmpVec);
PetscScalar *tmpValues;
@@ -1018,9 +1018,8 @@ namespace AMDiS {
MatMult(subdomain->getMatInterior(), ktest0, ktest1);
PetscScalar *valarray;
- VecGetArray(ktest1, &valarray);
-
Vec ktest2, ktest3;
+ VecGetArray(ktest1, &valarray);
VecCreateMPIWithArray(PETSC_COMM_WORLD, 1,
localDofMap.getRankDofs(), nGlobalOverallInterior,
valarray, &ktest2);
@@ -1046,14 +1045,26 @@ namespace AMDiS {
KSPSetNullSpace(ksp_feti, matNullSpace);
MatNullSpaceDestroy(&matNullSpace);
-#if (DEBUG != 0)
+ //#if (DEBUG != 0)
Vec vecSol;
- VecDuplicate(nullSpaceBasis, &vecSol);
+ VecDuplicate(nullSpaceBasis, &vecSol);
MatMult(mat_feti, nullSpaceBasis, vecSol);
- PetscReal norm;
+
+ Vec vecSol0, vecSol1;
+ VecNestGetSubVec(vecSol, 0, &vecSol0);
+ VecNestGetSubVec(vecSol, 1, &vecSol1);
+
+ PetscReal norm, norm0, norm1;
VecNorm(vecSol, NORM_2, &norm);
- MSG("NORM: %e\n", norm);
-#endif
+ VecNorm(vecSol0, NORM_2, &norm0);
+ VecNorm(vecSol1, NORM_2, &norm1);
+
+ MSG("NORM: %e (%e %e)\n", norm, norm0, norm1);
+ VecDestroy(&vecSol);
+ //#endif
+
+ // AMDiS::finalize();
+ // exit(0);
VecDestroy(&ktest0);
VecDestroy(&ktest1);
@@ -1692,6 +1703,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
+ debugNullSpace(vec);
+
// === Some temporary vectors. ===
Vec tmp_b0, tmp_b1;
@@ -2017,4 +2030,96 @@ namespace AMDiS {
&(nestMat[0]), &mat);
}
}
+
+
+ void PetscSolverFeti::debugNullSpace(SystemVector &vec)
+ {
+ FUNCNAME("PetscSolverFeti::debugNullSpace()");
+
+ TEST_EXIT(stokesMode)("This function works only for the stokes mode!\n");
+
+ Mat fetiMat;
+ createNestedFetiMat(fetiMat);
+
+ Vec vecArray[4];
+
+
+ Vec ktest0, ktest1, ktest2;
+ localDofMap.createLocalVec(ktest0);
+ localDofMap.createLocalVec(ktest1);
+ localDofMap.createVec(ktest2, nGlobalOverallInterior);
+
+ DofMap& m = localDofMap[pressureFeSpace].getMap();
+ for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
+ if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
+ int index = localDofMap.getLocalMatIndex(pressureComponent, it->first);
+ VecSetValue(ktest0, index, 1.0, INSERT_VALUES);
+ }
+ }
+ VecAssemblyBegin(ktest0);
+ VecAssemblyEnd(ktest0);
+ MatMult(subdomain->getMatInterior(), ktest0, ktest1);
+
+ PetscScalar *valarray;
+ VecGetArray(ktest0, &valarray);
+ VecCreateMPIWithArray(PETSC_COMM_WORLD, 1,
+ localDofMap.getRankDofs(), nGlobalOverallInterior,
+ valarray, &vecArray[0]);
+
+ Vec ktest3;
+ VecGetArray(ktest1, &valarray);
+ VecCreateMPIWithArray(PETSC_COMM_WORLD, 1,
+ localDofMap.getRankDofs(), nGlobalOverallInterior,
+ valarray, &ktest3);
+
+ primalDofMap.createVec(vecArray[1]);
+ VecSet(vecArray[1], 0.0);
+
+ interfaceDofMap.createVec(vecArray[2]);
+ VecSet(vecArray[2], 1.0);
+
+ lagrangeMap.createVec(vecArray[3]);
+ MatMult(subdomain->getMatInteriorCoarse(1), vecArray[2], ktest2);
+ VecAXPY(ktest2, 1.0, ktest3);
+ MatMult(mat_lagrange_scaled, ktest2, vecArray[3]);
+ VecScale(vecArray[3], -1.0);
+
+
+ Vec nullSpaceBasis;
+ VecCreateNest(mpiCommGlobal, 4, PETSC_NULL, vecArray, &nullSpaceBasis);
+
+ Vec vecSol;
+ VecDuplicate(nullSpaceBasis, &vecSol);
+
+ MatMult(fetiMat, nullSpaceBasis, vecSol);
+ PetscReal norm;
+ VecNorm(vecSol, NORM_2, &norm);
+ MSG("Null space norm: %e\n", norm);
+
+ Vec vec0;
+ Vec vec1;
+ Vec vec2;
+ Vec vec3;
+ VecNestGetSubVec(vecSol, 0, &vec0);
+ VecNestGetSubVec(vecSol, 1, &vec1);
+ VecNestGetSubVec(vecSol, 2, &vec2);
+ VecNestGetSubVec(vecSol, 3, &vec3);
+
+ PetscReal norm0, norm1, norm2, norm3;
+ VecNorm(vec0, NORM_2, &norm0);
+ VecNorm(vec1, NORM_2, &norm1);
+ VecNorm(vec2, NORM_2, &norm2);
+ VecNorm(vec3, NORM_2, &norm3);
+ MSG("Splitted null space norm: %e %e %e %e\n", norm0, norm1, norm2, norm3);
+
+
+ recoverSolution(vec0, vec1, vec);
+ recoverInterfaceSolution(vec2, vec);
+
+ VtkWriter::writeFile(&vec, "nullspace.vtu");
+
+ MatDestroy(&fetiMat);
+ VecDestroy(&nullSpaceBasis);
+ VecDestroy(&vecSol);
+ }
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index f32f281d034b615831542d557449898fd71e6062..465f44dcfb86206a423241554a9f307646028e4a 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -221,6 +221,9 @@ namespace AMDiS {
/// For debugging only!
void createNestedFetiMat(Mat &mat);
+ /// For debugging only!
+ void debugNullSpace(SystemVector &vec);
+
protected:
/// Mapping from primal DOF indices to a global index of primals.
ParallelDofMapping primalDofMap;