diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index 607eff4cf82684a5e1c4cff1ce002bac158df516..eb32b83657638428d7bc1d0570bec07c72dcec0c 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -457,10 +457,8 @@ namespace AMDiS { /// assignment operator Element& operator=(const Element& el); - /** \brief - * Checks whether the face with vertices dof[0],..,dof[DIM-1] is - * part of mel's boundary. returns the opposite vertex if true, -1 else - */ + /// Checks whether the face with vertices dof[0],..,dof[DIM-1] is + /// part of mel's boundary. returns the opposite vertex if true, -1 else int oppVertex(FixVec pdof) const; /// Refines Element's leaf data @@ -565,38 +563,30 @@ namespace AMDiS { /// to NULL for leaf elements. Element *child[2]; - /** \brief - * Vector of pointers to DOFs. These pointers must be available for elements - * vertices (for the geometric description of the mesh). There my be pointers - * for the edges, for faces and for the center of an element. They are - * ordered the following way: The first N_VERTICES entries correspond to the - * DOFs at the vertices of the element. The next ones are those at the edges, - * if present, then those at the faces, if present, and then those at the - * barycenter, if present. - */ + /// Vector of pointers to DOFs. These pointers must be available for elements + /// vertices (for the geometric description of the mesh). There my be pointers + /// for the edges, for faces and for the center of an element. They are + /// ordered the following way: The first N_VERTICES entries correspond to the + /// DOFs at the vertices of the element. The next ones are those at the edges, + /// if present, then those at the faces, if present, and then those at the + /// barycenter, if present. DegreeOfFreedom **dof; - /** \brief - * Unique global index of the element. these indices are not strictly ordered - * and may be larger than the number of elements in the binary tree (the list - * of indices may have holes after coarsening). - */ + /// Unique global index of the element. these indices are not strictly ordered + /// and may be larger than the number of elements in the binary tree (the list + /// of indices may have holes after coarsening). int index; - /** \brief - * Marker for refinement and coarsening. if mark is positive for a leaf - * element, this element is refined mark times. if mark is negative for - * a leaf element, this element is coarsened -mark times. - */ + /// Marker for refinement and coarsening. if mark is positive for a leaf + /// element, this element is refined mark times. if mark is negative for + /// a leaf element, this element is coarsened -mark times. int mark; - /** \brief - * If the element has a boundary edge on a curved boundary, this is a pointer - * to the coordinates of the new vertex that is created due to the refinement - * of the element, otherwise it is a NULL pointer. Thus coordinate - * information can be also produced by the traversal routines in the case of - * curved boundary. - */ + /// If the element has a boundary edge on a curved boundary, this is a pointer + /// to the coordinates of the new vertex that is created due to the refinement + /// of the element, otherwise it is a NULL pointer. Thus coordinate + /// information can be also produced by the traversal routines in the case of + /// curved boundary. WorldVector *newCoord; /// Pointer to the Mesh this element belongs to diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index 37c8fed40aa7a8ac666cd593bdb75743bd184859..93b6631e1764903e39d3a85eb6d04f5139abd5a0 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -699,6 +699,22 @@ namespace AMDiS { } } + bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace) + { + Element *el = edge.el; + int i0 = el->getVertexOfEdge(edge.ithObj, 0); + int i1 = el->getVertexOfEdge(edge.ithObj, 1); + DegreeOfFreedom d0 = el->getDof(i0, 0); + DegreeOfFreedom d1 = el->getDof(i1, 0); + WorldVector c0, c1; + el->getMesh()->getDofIndexCoords(d0, feSpace, c0); + el->getMesh()->getDofIndexCoords(d1, feSpace, c1); + bool xe = fabs(c0[0] - c1[0]) < 1e-8; + bool ye = fabs(c0[1] - c1[1]) < 1e-8; + bool ze = fabs(c0[2] - c1[2]) < 1e-8; + int counter = static_cast(xe) + static_cast(ye) + static_cast(ze); + return (counter == 2); + } void PetscSolverFeti::createMatAugmentedLagrange(vector &feSpaces) { @@ -709,13 +725,16 @@ namespace AMDiS { double wtime = MPI::Wtime(); - nRankEdges = 0; nOverallEdges = 0; InteriorBoundary &intBound = meshDistributor->getIntBoundary(); + std::set allEdges; for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) - if (it->rankObj.subObj == EDGE) - nRankEdges++; - + if (it->rankObj.subObj == EDGE && + testWirebasketEdge(it->rankObj, feSpaces[0]) && + allEdges.count(it->rankObj) == 0) + allEdges.insert(it->rankObj); + + nRankEdges = allEdges.size(); int rStartEdges = 0; mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges); @@ -733,30 +752,30 @@ namespace AMDiS { MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); int rowCounter = rStartEdges; - for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) { - if (it->rankObj.subObj == EDGE) { - for (int i = 0; i < feSpaces.size(); i++) { - DofContainer edgeDofs; - it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs); - - MSG("SIZE = %d\n", edgeDofs.size()); - - for (DofContainer::iterator it = edgeDofs.begin(); - it != edgeDofs.end(); it++) { - TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false) - ("Should not be primal!\n"); - - if (dirichletRows[feSpaces[i]].count(**it)) - continue; + for (std::set::iterator edgeIt = allEdges.begin(); + edgeIt != allEdges.end(); ++edgeIt) { - int col = lagrangeMap.getMatIndex(i, **it); - double value = 1.0; - MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES); - } - - rowCounter++; + for (int i = 0; i < feSpaces.size(); i++) { + DofContainer edgeDofs; + edgeIt->el->getAllDofs(feSpaces[i], *edgeIt, edgeDofs); + + MSG("SIZE = %d\n", edgeDofs.size()); + + for (DofContainer::iterator it = edgeDofs.begin(); + it != edgeDofs.end(); ++it) { + TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false) + ("Should not be primal!\n"); + + if (dirichletRows[feSpaces[i]].count(**it)) + continue; + + int col = lagrangeMap.getMatIndex(i, **it); + double value = 1.0; + MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES); } - } + + rowCounter++; + } } MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY); @@ -1013,7 +1032,38 @@ namespace AMDiS { MatShellSetOperation(tmp, MATOP_MULT, (void(*)(void))petscMultMatSchurPrimalAugmented); + schurPrimalAugmentedData.testMode = false; MatComputeExplicitOperator(tmp, &mat_schur_primal); + schurPrimalAugmentedData.testMode = true; + + { + Vec tvec0, tvec1; + VecCreateMPI(PETSC_COMM_WORLD, + primalDofMap.getRankDofs() + nRankEdges, + primalDofMap.getOverallDofs() + nOverallEdges, + &tvec0); + VecDuplicate(tvec0, &tvec1); + + if (meshDistributor->getMpiRank() == 0) + VecSetValue(tvec0, 12, 0.0, INSERT_VALUES); + + VecAssemblyBegin(tvec0); + VecAssemblyEnd(tvec1); + + MatMult(tmp, tvec0, tvec1); + + // VecView(tvec1, PETSC_VIEWER_STDOUT_WORLD); + + PetscReal n, a, b; + VecNorm(tvec1, NORM_2, &n); + VecMax(tvec1, PETSC_NULL, &a); + VecMin(tvec1, PETSC_NULL, &b); + + MSG("DIE WERTE SIND %e %e %e\n", n, a, b); + } + + + // MatView(mat_schur_primal, PETSC_VIEWER_STDOUT_WORLD); MatDestroy(&tmp); schurPrimalAugmentedData.subSolver = NULL; diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index 0d73abb4d0ee10c10f3830d53140229db2765b72..f52cf2532ac51f1a451d8cb896742a9126cea0aa 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -133,6 +133,8 @@ namespace AMDiS { void createMatAugmentedLagrange(vector &feSpaces); + bool testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace); + /// void createPreconditionerMatrix(Matrix *mat); diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc index cd9c94c38a204d6490149ad095f30b4a3ae6af33..f035848bc5cb5ec478d133d687da819a380f7a22 100644 --- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc +++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc @@ -66,6 +66,11 @@ namespace AMDiS { PetscInt muLocalSize = allLocalSize - primalLocalSize; PetscInt muSize = allSize - primalSize; + if (data->testMode) { + VecView(x, PETSC_VIEWER_STDOUT_WORLD); + MSG("LOCAL SIZE: %d %d %d\n", allLocalSize, primalLocalSize, muLocalSize); + } + VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &x_mu); VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &y_mu); @@ -76,10 +81,17 @@ namespace AMDiS { VecGetArray(x_primal, &primalValue); VecGetArray(x_mu, &muValue); - for (int i = 0; i < primalLocalSize; i++) + for (int i = 0; i < primalLocalSize; i++) { + if (data->testMode && allValue[i] != 1.0) + MSG("ONE IS in %d ith local primal value!\n", i); primalValue[i] = allValue[i]; - for (int i = 0; i < muLocalSize; i++) + } + for (int i = 0; i < muLocalSize; i++) { + if (data->testMode && allValue[primalLocalSize + i] != 0.0) + MSG("ONE IS in %d ith local mu value!\n", primalLocalSize + i); + muValue[i] = allValue[primalLocalSize + i]; + } VecRestoreArray(x, &allValue); VecRestoreArray(x_primal, &primalValue); diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h index 8fa97340ebbf310d156b691c365cbba3e1659183..ff116a68729fb6cf8e7618d30a5d141e74e39cf4 100644 --- a/AMDiS/src/parallel/PetscSolverFetiStructs.h +++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h @@ -65,6 +65,8 @@ namespace AMDiS { PetscSolver* subSolver; bool nestedVec; + + bool testMode; };