From 098221a6cbffc2f68aabedc48202d5380689a732 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Fri, 9 Mar 2012 15:21:46 +0000 Subject: [PATCH] Have a nice weekend. --- AMDiS/src/ProblemStat.cc | 11 -- AMDiS/src/parallel/FeSpaceMapping.cc | 28 ++-- AMDiS/src/parallel/FeSpaceMapping.h | 51 +++++-- AMDiS/src/parallel/PetscSolverFeti.cc | 206 +++++++++++++------------- 4 files changed, 151 insertions(+), 145 deletions(-) diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc index 992dcf77..dfbb4a2d 100644 --- a/AMDiS/src/ProblemStat.cc +++ b/AMDiS/src/ProblemStat.cc @@ -870,17 +870,6 @@ namespace AMDiS { INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first, clock())); #endif - -#if 0 - VtkWriter::writeFile(rhs, "rhs0.vtu"); - for (int i = 0; i < nComponents; i++) { - double s = rhs->getDOFVector(i)->sum(); - s /= -static_cast(rhs->getDOFVector(i)->getUsedSize()); - MSG("MOVE RHS: %f\n", s); - add(*(rhs->getDOFVector(i)), s, *(rhs->getDOFVector(i))); - } - VtkWriter::writeFile(rhs, "rhs1.vtu"); -#endif } diff --git a/AMDiS/src/parallel/FeSpaceMapping.cc b/AMDiS/src/parallel/FeSpaceMapping.cc index 4e548659..05bfb5d4 100644 --- a/AMDiS/src/parallel/FeSpaceMapping.cc +++ b/AMDiS/src/parallel/FeSpaceMapping.cc @@ -26,17 +26,17 @@ namespace AMDiS { } - void GlobalDofMap::update(bool add) + void GlobalDofMap::update() { - for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it) - if (it->second == -1 && nonRankDofs.count(it->first) == 0) - it->second = nRankDofs++; + for (map::iterator it = dofMap.begin(); it != dofMap.end(); ++it) + if (it->second.local == -1 && nonRankDofs.count(it->first) == 0) + it->second.local = nRankDofs++; nOverallDofs = 0; rStartDofs = 0; mpi::getDofNumbering(*mpiComm, nRankDofs, rStartDofs, nOverallDofs); - if (add) + if (needGlobalMapping) addOffset(rStartDofs); if (overlap) @@ -48,8 +48,8 @@ namespace AMDiS { void GlobalDofMap::addOffset(int offset) { - for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it) - it->second += offset; + for (map::iterator it = dofMap.begin(); it != dofMap.end(); ++it) + it->second.local += offset; } void GlobalDofMap::computeNonLocalIndices() @@ -62,7 +62,7 @@ namespace AMDiS { for (DofComm::Iterator it(*sendDofs, feSpace); !it.end(); it.nextRank()) for (; !it.endDofIter(); it.nextDof()) if (dofMap.count(it.getDofIndex()) && !nonRankDofs.count(it.getDofIndex())) - stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()]); + stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].local); stdMpi.updateSendDataSize(); @@ -86,7 +86,7 @@ namespace AMDiS { int i = 0; for (; !it.endDofIter(); it.nextDof()) if (nonRankDofs.count(it.getDofIndex())) - dofMap[it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[i++]; + dofMap[it.getDofIndex()].local = stdMpi.getRecvData(it.getRank())[i++]; } } @@ -97,9 +97,9 @@ namespace AMDiS { map dofToMatIndex; - for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { + for (map::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { if (!nonRankDofs.count(it->first)) { - int globalMatIndex = it->second - rStartDofs + offset; + int globalMatIndex = it->second.local - rStartDofs + offset; dofToMatIndex[it->first] = globalMatIndex; } } @@ -143,10 +143,10 @@ namespace AMDiS { MSG("Local to global mapping on this rank: \n"); - for (DofMapping::iterator it = dofMap.begin(); it != dofMap.end(); ++it) + for (map::iterator it = dofMap.begin(); it != dofMap.end(); ++it) if (nonRankDofs.count(it->first) == 0) - MSG(" %d -> %d (rank-dof)\n", it->first, it->second); + MSG(" %d -> %d (rank-dof)\n", it->first, it->second.local); else - MSG(" %d -> %d \n", it->first, it->second); + MSG(" %d -> %d \n", it->first, it->second.local); } } diff --git a/AMDiS/src/parallel/FeSpaceMapping.h b/AMDiS/src/parallel/FeSpaceMapping.h index d75e1a2e..ae172dee 100644 --- a/AMDiS/src/parallel/FeSpaceMapping.h +++ b/AMDiS/src/parallel/FeSpaceMapping.h @@ -33,6 +33,11 @@ namespace AMDiS { using namespace std; + struct MultiIndex + { + int local, global; + }; + class GlobalDofMap { public: @@ -48,6 +53,7 @@ namespace AMDiS { feSpace(NULL), sendDofs(NULL), recvDofs(NULL), + needGlobalMapping(false), nRankDofs(0), nOverallDofs(0), rStartDofs(0), @@ -56,20 +62,23 @@ namespace AMDiS { void clear(); - DegreeOfFreedom operator[](DegreeOfFreedom d) + MultiIndex& operator[](DegreeOfFreedom d) { TEST_EXIT_DBG(dofMap.count(d))("Should not happen!\n"); return dofMap[d]; } - void insertRankDof(DegreeOfFreedom dof0) + void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1) { FUNCNAME("GlobalDofMap::insertRankDof()"); TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n"); - dofMap[dof0] = -1; + dofMap[dof0].local = dof1; + + if (dof1 != -1) + nRankDofs++; } void insert(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1) @@ -78,7 +87,7 @@ namespace AMDiS { TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n"); - dofMap[dof0] = dof1; + dofMap[dof0].local = dof1; nonRankDofs.insert(dof0); } @@ -93,12 +102,12 @@ namespace AMDiS { return dofMap.size(); } - DofMapping& getMap() + map& getMap() { return dofMap; } - void update(bool add = true); + void update(); void addOffset(int offset); @@ -124,19 +133,26 @@ namespace AMDiS { recvDofs = &pRecv; } + void setNeedGlobalMapping(bool b) + { + needGlobalMapping = b; + } + private: MPI::Intracomm* mpiComm; const FiniteElemSpace *feSpace; /// - DofMapping dofMap; + map dofMap; std::set nonRankDofs; DofComm *sendDofs; DofComm *recvDofs; + + bool needGlobalMapping; public: /// int nRankDofs, nOverallDofs, rStartDofs; @@ -157,11 +173,6 @@ namespace AMDiS { rStartDofs(-1) {} - void setMpiComm(MPI::Intracomm *m) - { - mpiComm = m; - } - T& operator[](const FiniteElemSpace* feSpace) { FUNCNAME("FeSpaceData::operator[]()"); @@ -239,15 +250,23 @@ namespace AMDiS { return rStartDofs; } - void setFeSpaces(vector &fe) + void init(MPI::Intracomm *m, + vector &fe, + bool needGlobalMapping) { + mpiComm = m; feSpaces = fe; - for (unsigned int i = 0; i < feSpaces.size(); i++) + for (unsigned int i = 0; i < feSpaces.size(); i++) { addFeSpace(feSpaces[i]); + data[feSpaces[i]].setNeedGlobalMapping(needGlobalMapping); + } } void update() { +// for (unsigned int i = 0; i < feSpaces.size(); i++) +// data[feSpaces[i]].update(); + nRankDofs = getRankDofs(feSpaces); nOverallDofs = getOverallDofs(feSpaces); rStartDofs = getStartDofs(feSpaces); @@ -258,7 +277,7 @@ namespace AMDiS { int result = 0; for (int i = 0; i < ithFeSpace; i++) result += data[feSpaces[i]].nRankDofs; - result += data[feSpaces[ithFeSpace]][index]; + result += data[feSpaces[ithFeSpace]][index].local; return result; } @@ -277,7 +296,7 @@ namespace AMDiS { int result = rStartDofs; for (int i = 0; i < ithFeSpace; i++) result += data[feSpaces[i]].nRankDofs; - result += data[feSpaces[ithFeSpace]][index]; + result += data[feSpaces[ithFeSpace]][index].local; return result; } diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc index bf407e0a..edc5167f 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.cc +++ b/AMDiS/src/parallel/PetscSolverFeti.cc @@ -224,11 +224,76 @@ namespace AMDiS { createLagrange(feSpace); createIndexB(feSpace); } + + primalDofMap.update(); + dualDofMap.update(); + lagrangeMap.update(); + localDofMap.update(); + + for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); + + MSG("nRankPrimals = %d nOverallPrimals = %d\n", + primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs); + + MSG("nRankDuals = %d nOverallDuals = %d\n", + dualDofMap[feSpace].nRankDofs, + dualDofMap[feSpace].nOverallDofs); + + MSG("nRankLagrange = %d nOverallLagrange = %d\n", + lagrangeMap[feSpace].nRankDofs, + lagrangeMap[feSpace].nOverallDofs); + } + + // === Communicate lagrangeMap to all other ranks. === + + for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) { + const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i); + + StdMpi > stdMpi(meshDistributor->getMpiComm()); + + for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); + !it.end(); it.nextRank()) { + for (; !it.endDofIter(); it.nextDof()) + if (!isPrimal(feSpace, it.getDofIndex())) { + DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()].local; + stdMpi.getSendData(it.getRank()).push_back(d); + } + } + + stdMpi.updateSendDataSize(); + + for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); + !it.end(); it.nextRank()) { + bool recvData = false; + for (; !it.endDofIter(); it.nextDof()) + if (!isPrimal(feSpace, it.getDofIndex())) { + recvData = true; + break; + } + + if (recvData) + stdMpi.recv(it.getRank()); + } + + stdMpi.startCommunication(); + + for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); + !it.end(); it.nextRank()) { + int counter = 0; + for (; !it.endDofIter(); it.nextDof()) { + if (!isPrimal(feSpace, it.getDofIndex())) { + DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++]; + lagrangeMap[feSpace].insert(it.getDofIndex(), d); + } + } + } + } } void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace) - { + { FUNCNAME("PetscSolverFeti::createPrimals()"); // === Define all vertices on the interior boundaries of the macro mesh === @@ -257,14 +322,8 @@ namespace AMDiS { primalDofMap[feSpace].setOverlap(true); primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(), meshDistributor->getRecvDofs()); - primalDofMap[feSpace].update(); - MSG("nRankPrimals = %d nOverallPrimals = %d\n", - primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs); - - TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size()) - ("Number of primals %d, but number of global primals on this rank is %d!\n", - primals.size(), primalDofMap[feSpace].size()); + primalDofMap[feSpace].update(); } @@ -335,11 +394,7 @@ namespace AMDiS { if (!isPrimal(feSpace, **it)) dualDofMap[feSpace].insertRankDof(**it); - dualDofMap[feSpace].update(false); - - MSG("nRankDuals = %d nOverallDuals = %d\n", - dualDofMap[feSpace].nRankDofs, - dualDofMap[feSpace].nOverallDofs); + dualDofMap[feSpace].update(); } @@ -351,8 +406,8 @@ namespace AMDiS { // === appropriate number of Lagrange constraints. === int nRankLagrange = 0; - DofMapping& dualMap = dualDofMap[feSpace].getMap(); - for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { + map& dualMap = dualDofMap[feSpace].getMap(); + for (map::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { if (meshDistributor->getIsRankDof(feSpace, it->first)) { lagrangeMap[feSpace].insert(it->first, nRankLagrange); int degree = boundaryDofRanks[it->first].size(); @@ -361,52 +416,6 @@ namespace AMDiS { } lagrangeMap[feSpace].nRankDofs = nRankLagrange; lagrangeMap[feSpace].update(); - - MSG("nRankLagrange = %d nOverallLagrange = %d\n", - lagrangeMap[feSpace].nRankDofs, - lagrangeMap[feSpace].nOverallDofs); - - - // === Communicate lagrangeMap to all other ranks. === - - StdMpi > stdMpi(meshDistributor->getMpiComm()); - - for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); - !it.end(); it.nextRank()) { - for (; !it.endDofIter(); it.nextDof()) - if (!isPrimal(feSpace, it.getDofIndex())) { - DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()]; - stdMpi.getSendData(it.getRank()).push_back(d); - } - } - - stdMpi.updateSendDataSize(); - - for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); - !it.end(); it.nextRank()) { - bool recvData = false; - for (; !it.endDofIter(); it.nextDof()) - if (!isPrimal(feSpace, it.getDofIndex())) { - recvData = true; - break; - } - - if (recvData) - stdMpi.recv(it.getRank()); - } - - stdMpi.startCommunication(); - - for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); - !it.end(); it.nextRank()) { - int counter = 0; - for (; !it.endDofIter(); it.nextDof()) { - if (!isPrimal(feSpace, it.getDofIndex())) { - DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++]; - lagrangeMap[feSpace].insert(it.getDofIndex(), d); - } - } - } } @@ -425,13 +434,11 @@ namespace AMDiS { if (admin->isDofFree(i) == false && isPrimal(feSpace, i) == false && dualDofMap[feSpace].isSet(i) == false) { - localDofMap[feSpace].insertRankDof(i); + localDofMap[feSpace].insertRankDof(i, nLocalInterior); nLocalInterior++; } } - localDofMap[feSpace].update(false); - TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + dualDofMap[feSpace].size() == static_cast(admin->getUsedDofs())) @@ -439,14 +446,11 @@ namespace AMDiS { // === And finally, add the global indicies of all dual nodes. === - for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin(); + for (map::iterator it = dualDofMap[feSpace].getMap().begin(); it != dualDofMap[feSpace].getMap().end(); ++it) localDofMap[feSpace].insertRankDof(it->first); - localDofMap[feSpace].update(false); - - dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + - nLocalInterior); + localDofMap[feSpace].update(); } @@ -471,13 +475,13 @@ namespace AMDiS { // === m == r, than the rank sets -1.0 for the corresponding === // === constraint. === - DofMapping &dualMap = dualDofMap[feSpace].getMap(); - for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { + map &dualMap = dualDofMap[feSpace].getMap(); + for (map::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { TEST_EXIT_DBG(boundaryDofRanks.count(it->first)) ("Should not happen!\n"); // Global index of the first Lagrange constriant for this node. - int index = lagrangeMap[feSpace][it->first]; + int index = lagrangeMap[feSpace][it->first].local; // Copy set of all ranks that contain this dual node. vector W(boundaryDofRanks[it->first].begin(), boundaryDofRanks[it->first].end()); @@ -490,7 +494,10 @@ namespace AMDiS { // Set the constraint for all components of the system. for (int k = 0; k < nComponents; k++) { int rowIndex = index * nComponents + k; - int colIndex = it->second * nComponents + k; + // int colIndex = it->second * nComponents + k; + int colIndex = + (localDofMap[feSpace].rStartDofs + + localDofMap[feSpace][it->first].local) * nComponents + k; double value = (W[i] == mpiRank ? 1.0 : -1.0); MatSetValue(mat_lagrange, rowIndex, colIndex, value, INSERT_VALUES); @@ -829,10 +836,10 @@ namespace AMDiS { { int counter = 0; - for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin(); + for (map::iterator it = primalDofMap[feSpace].getMap().begin(); it != primalDofMap[feSpace].getMap().end(); ++it) { for (int i = 0; i < nComponents; i++) { - globalIsIndex.push_back(it->second * nComponents + i); + globalIsIndex.push_back(it->second.local * nComponents + i); localIsIndex.push_back(counter++); } } @@ -873,14 +880,14 @@ namespace AMDiS { for (int i = 0; i < nComponents; i++) { DOFVector& dofVec = *(vec.getDOFVector(i)); - for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin(); + for (map::iterator it = localDofMap[feSpace].getMap().begin(); it != localDofMap[feSpace].getMap().end(); ++it) { int petscIndex = localDofMap.mapLocal(it->first, i); dofVec[it->first] = localSolB[petscIndex]; } int counter = 0; - for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin(); + for (map::iterator it = primalDofMap[feSpace].getMap().begin(); it != primalDofMap[feSpace].getMap().end(); ++it) { dofVec[it->first] = localSolPrimal[counter * nComponents + i]; counter++; @@ -904,23 +911,13 @@ namespace AMDiS { // === Create all sets and indices. === - primalDofMap.setMpiComm(mpiComm); - dualDofMap.setMpiComm(mpiComm); - lagrangeMap.setMpiComm(mpiComm); - localDofMap.setMpiComm(mpiComm); - - primalDofMap.setFeSpaces(feSpaces); - dualDofMap.setFeSpaces(feSpaces); - lagrangeMap.setFeSpaces(feSpaces); - localDofMap.setFeSpaces(feSpaces); + primalDofMap.init(mpiComm, feSpaces, true); + dualDofMap.init(mpiComm, feSpaces, false); + lagrangeMap.init(mpiComm, feSpaces, true); + localDofMap.init(mpiComm, feSpaces, false); updateDofData(); - primalDofMap.update(); - dualDofMap.update(); - lagrangeMap.update(); - localDofMap.update(); - // === Create matrices for the FETI-DP method. === int nRowsRankB = localDofMap.getRankDofs(); @@ -1051,8 +1048,8 @@ namespace AMDiS { // === For preconditioner === if (!rowPrimal && !colPrimal) { - int rowIndex = localDofMap[feSpace][*cursor]; - int colIndex = localDofMap[feSpace][col(*icursor)]; + int rowIndex = localDofMap[feSpace][*cursor].local; + int colIndex = localDofMap[feSpace][col(*icursor)].local; if (rowIndex < nLocalInterior) { if (colIndex < nLocalInterior) { @@ -1060,7 +1057,8 @@ namespace AMDiS { colsLocal.push_back(colIndex2); valuesLocal.push_back(value(*icursor)); } else { - int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) * + int colIndex2 = + (localDofMap[feSpace][col(*icursor)].local - nLocalInterior) * nComponents + j; colsLocalOther.push_back(colIndex2); @@ -1072,7 +1070,7 @@ namespace AMDiS { colsLocalOther.push_back(colIndex2); valuesLocalOther.push_back(value(*icursor)); } else { - int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) * + int colIndex2 = (localDofMap[feSpace][col(*icursor)].local - nLocalInterior) * nComponents + j; colsLocal.push_back(colIndex2); @@ -1089,10 +1087,10 @@ namespace AMDiS { if (rowPrimal) { int rowIndex = - primalDofMap[feSpace][*cursor] * nComponents + i; + primalDofMap[feSpace][*cursor].local * nComponents + i; for (unsigned int k = 0; k < cols.size(); k++) - cols[k] = primalDofMap[feSpace][cols[k]] * nComponents + j; + cols[k] = primalDofMap[feSpace][cols[k]].local * nComponents + j; MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); @@ -1118,7 +1116,7 @@ namespace AMDiS { for (unsigned int k = 0; k < colsOther.size(); k++) colsOther[k] = - primalDofMap[feSpace][colsOther[k]] * nComponents + j; + primalDofMap[feSpace][colsOther[k]].local * nComponents + j; MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); @@ -1131,10 +1129,10 @@ namespace AMDiS { switch (fetiPreconditioner) { case FETI_DIRICHLET: { - int rowIndex = localDofMap[feSpace][*cursor]; + int rowIndex = localDofMap[feSpace][*cursor].local; if (rowIndex < nLocalInterior) { - int rowIndex2 = localDofMap[feSpace][*cursor] * nComponents + i; + int rowIndex2 = localDofMap[feSpace][*cursor].local * nComponents + i; MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); @@ -1144,7 +1142,7 @@ namespace AMDiS { &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES); } else { int rowIndex2 = - (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i; + (localDofMap[feSpace][*cursor].local - nLocalInterior) * nComponents + i; MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); @@ -1159,11 +1157,11 @@ namespace AMDiS { case FETI_LUMPED: { - int rowIndex = localDofMap[feSpace][*cursor]; + int rowIndex = localDofMap[feSpace][*cursor].local; if (rowIndex >= nLocalInterior) { int rowIndex2 = - (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i; + (localDofMap[feSpace][*cursor].local - nLocalInterior) * nComponents + i; MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(), &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES); @@ -1255,7 +1253,7 @@ namespace AMDiS { for (dofIt.reset(); !dofIt.end(); ++dofIt) { int index = dofIt.getDOFIndex(); if (isPrimal(feSpace, index)) { - index = primalDofMap[feSpace][index] * nComponents + i; + index = primalDofMap[feSpace][index].local * nComponents + i; double value = *dofIt; VecSetValues(f_primal, 1, &index, &value, ADD_VALUES); } else { -- GitLab