diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index 36899fe296b3a533f9fb166e0b94412e361d0aba..d0d091cdfdfb20fadec83b962258d01ab7e0b7a7 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -273,7 +273,14 @@ namespace AMDiS { for (DofContainerSet::iterator it = vertices.begin(); it != vertices.end(); ++it) { if (meshLevel == 0) { - primals.insert(**it); + MSG("WARNING: Modified primal detection algorithm!\n"); + + double e = 1e-3; + WorldVector<double> c; + feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c); + + if (fabs(c[0]) > e && fabs(c[1]) > e && fabs(c[0] - 1.0) > e && fabs(c[1] - 1.0) > e) + primals.insert(**it); } else { double e = 1e-8; WorldVector<double> c; @@ -909,10 +916,11 @@ namespace AMDiS { lumpedData->localToDualMap[matIndexLocal] = matIndexDual; } } - + if (stokesMode) { DOFVector<double> tmpvec(pressureFeSpace, "tmpvec"); createH2vec(tmpvec); + VtkWriter::writeFile(&tmpvec, "tmpvec.vtu"); interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec); DofMap& m = interfaceDofMap[pressureFeSpace].getMap(); @@ -926,14 +934,23 @@ namespace AMDiS { VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec); VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec); - + localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1); + primalDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_primal); fetiInterfaceLumpedPreconData.subSolver = subdomain; } KSPGetPC(ksp_feti, &precon_feti); PCSetType(precon_feti, PCSHELL); if (stokesMode) { + KSPCreate(PETSC_COMM_WORLD, &fetiInterfaceLumpedPreconData.ksp_primal); + KSPSetOperators(fetiInterfaceLumpedPreconData.ksp_primal, + subdomain->getMatCoarse(0, 0), + subdomain->getMatCoarse(0, 0), + SAME_NONZERO_PATTERN); + KSPSetOptionsPrefix(fetiInterfaceLumpedPreconData.ksp_primal, "primal_"); + KSPSetFromOptions(fetiInterfaceLumpedPreconData.ksp_primal); + PCShellSetContext(precon_feti, static_cast<void*>(&fetiInterfaceLumpedPreconData)); PCShellSetApply(precon_feti, petscApplyFetiInterfaceLumpedPrecon); } else { @@ -1968,26 +1985,40 @@ namespace AMDiS { Mesh *mesh = vec.getFeSpace()->getMesh(); TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n"); + int n0 = vec.getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX); TraverseStack stack; ElInfo *elInfo = - stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); + stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET); while (elInfo) { Element *el = elInfo->getElement(); for (int i = 0; i < 3; i++) { - int d0 = el->getVertexOfEdge(i, 0); - int d1 = el->getVertexOfEdge(i, 1); + DegreeOfFreedom dof = el->getDof(i, n0); + vec[dof] += elInfo->getDet() * 0.5; + cvec[dof] += 1; + } + - WorldVector<double> c0 = elInfo->getCoord(d0); - WorldVector<double> c1 = elInfo->getCoord(d1); + /* + + for (int i = 0; i < 3; i++) { + int idx0 = el->getVertexOfEdge(i, 0); + int idx1 = el->getVertexOfEdge(i, 1); + DegreeOfFreedom dof0 = el->getDof(idx0, n0); + DegreeOfFreedom dof1 = el->getDof(idx1, n0); + + WorldVector<double> c0 = elInfo->getCoord(idx0); + WorldVector<double> c1 = elInfo->getCoord(idx1); c0 -= c1; double length = norm(c0); - vec[el->getDof(d0, 0)] += length; - vec[el->getDof(d1, 0)] += length; - cvec[el->getDof(d0, 0)] += 1; - cvec[el->getDof(d1, 0)] += 1; + vec[dof0] += length; + vec[dof1] += length; + cvec[dof0] += 1; + cvec[dof1] += 1; } - + */ + + elInfo = stack.traverseNext(elInfo); } diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc index d9c43ce36865365a45d0884eb0d6daaf10e87b46..582053c82bcbf37bd514bfa845e8480d39d8ebda 100644 --- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc +++ b/AMDiS/src/parallel/PetscSolverFetiOperators.cc @@ -299,17 +299,22 @@ namespace AMDiS { VecNestGetSubVec(xvec, 1, &x_lagrange); VecNestGetSubVec(yvec, 0, &y_interface); VecNestGetSubVec(yvec, 1, &y_lagrange); - - - // VecCopy(x_interface, y_interface); - // VecScale(y_interface, 100.0); - VecPointwiseMult(x_interface, data->h2vec, y_interface); +#if 0 + MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_primal); + KSPSolve(data->ksp_primal, data->tmp_primal, data->tmp_primal); + MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_primal, y_interface); - // Multiply with scaled Lagrange constraint matrix. MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b0); MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b1); VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1); +#else + VecCopy(x_interface, y_interface); + VecScale(y_interface, 1024.0); + // VecPointwiseMult(y_interface, x_interface, data->h2vec); + + MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0); +#endif // === Restriction of the B nodes to the boundary nodes. === @@ -352,7 +357,12 @@ namespace AMDiS { // Multiply with scaled Lagrange constraint matrix. MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, y_lagrange); - // MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface); +#if 0 + Vec tmp_interface; + VecDuplicate(y_interface, &tmp_interface); + MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, tmp_interface); + VecAXPY(y_interface, 1.0, tmp_interface); +#endif return 0; } diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h index 820b2cf9b51a58341eece6ebd44569f384198c8b..3b1cce80ab8a46680735529d367fc860156ce9a0 100644 --- a/AMDiS/src/parallel/PetscSolverFetiStructs.h +++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h @@ -121,6 +121,10 @@ namespace AMDiS { PetscSolver* subSolver; Vec h2vec; + + KSP ksp_primal; + + Vec tmp_primal; };