47 #ifndef PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
48 #define PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
56 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
60 validParamList->
set<
double>(
"Advanced Dirichlet: threshold", 1e-5,
"Drop tolerance for Dirichlet detection");
61 validParamList->
set<
double>(
"Variable DOF amalgamation: threshold", 1.8e-9,
"Drop tolerance for amalgamation process");
62 validParamList->
set<
int>(
"maxDofPerNode", 1,
"Maximum number of DOFs per node");
67 return validParamList;
70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 Input(currentLevel,
"A");
76 Input(currentLevel,
"Coordinates");
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
97 RCP<dxMV> Coords = Get<RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO,
NO> > >(currentLevel,
"Coordinates");
99 int maxDofPerNode = pL.
get<
int>(
"maxDofPerNode");
100 Scalar dirDropTol = Teuchos::as<Scalar>(pL.
get<
double>(
"Advanced Dirichlet: threshold"));
101 Scalar amalgDropTol = Teuchos::as<Scalar>(pL.
get<
double>(
"Variable DOF amalgamation: threshold"));
103 bool bHasZeroDiagonal =
false;
121 std::vector<LocalOrdinal> map(A->getLocalNumRows());
122 this->buildPaddedMap(dofPresent, map, A->getLocalNumRows());
125 std::vector<LocalOrdinal> myLocalNodeIds(A->getColMap()->getLocalNumElements());
128 size_t nLocalNodes, nLocalPlusGhostNodes;
129 this->assignGhostLocalNodeIds(A->getRowMap(), A->getColMap(), myLocalNodeIds, map, maxDofPerNode, nLocalNodes, nLocalPlusGhostNodes, comm);
133 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dofPresent.
size()) != Teuchos::as<size_t>(nLocalNodes * maxDofPerNode),
MueLu::Exceptions::RuntimeError,
"VariableDofLaplacianFactory: size of provided DofPresent array is " << dofPresent.
size() <<
" but should be " << nLocalNodes * maxDofPerNode <<
" on the current processor.");
143 size_t nLocalDofs = A->getRowMap()->getLocalNumElements();
144 size_t nLocalPlusGhostDofs = A->getColMap()->getLocalNumElements();
153 if (nLocalDofs > 0) {
154 amalgRowMapGIDs[count] = myGids[0];
155 amalgColMapGIDs[count] = myGids[0];
159 for (
size_t i = 1; i < nLocalDofs; i++) {
160 if (myLocalNodeIds[i] != myLocalNodeIds[i - 1]) {
161 amalgRowMapGIDs[count] = myGids[i];
162 amalgColMapGIDs[count] = myGids[i];
167 RCP<GOVector> tempAmalgColVec = GOVectorFactory::Build(A->getDomainMap());
170 for (
size_t i = 0; i < A->getDomainMap()->getLocalNumElements(); i++)
171 tempAmalgColVecData[i] = amalgColMapGIDs[myLocalNodeIds[i]];
174 RCP<GOVector> tempAmalgColVecTarget = GOVectorFactory::Build(A->getColMap());
176 tempAmalgColVecTarget->doImport(*tempAmalgColVec, *dofImporter,
Xpetra::INSERT);
181 for (
size_t i = 0; i < myLocalNodeIds.size(); i++)
182 amalgColMapGIDs[myLocalNodeIds[i]] = tempAmalgColVecBData[i];
188 A->getRowMap()->getIndexBase(),
194 A->getRangeMap()->getIndexBase(),
205 size_t nNonZeros = 0;
206 std::vector<bool> isNonZero(nLocalPlusGhostDofs,
false);
207 std::vector<size_t> nonZeroList(nLocalPlusGhostDofs);
212 A->getLocalDiagCopy(*diagVecUnique);
219 Acrs->getAllValues(rowptr, colind, values);
230 amalgRowPtr[0] = newNzs;
232 bool doNotDrop =
false;
234 if (values.size() == 0) doNotDrop =
true;
236 for (decltype(rowptr.size()) i = 0; i < rowptr.size() - 1; i++) {
237 blockRow = std::floor<LocalOrdinal>(map[i] / maxDofPerNode);
238 if (blockRow != oldBlockRow) {
240 for (
size_t j = 0; j < nNonZeros; j++) isNonZero[nonZeroList[j]] =
false;
242 amalgRowPtr[blockRow] = newNzs;
244 for (
size_t j = rowptr[i]; j < rowptr[i + 1]; j++) {
245 if (doNotDrop ==
true ||
246 (STS::magnitude(values[j] / STS::magnitude(sqrt(STS::magnitude(diagVecData[i]) * STS::magnitude(diagVecData[colind[j]])))) >= STS::magnitude(amalgDropTol))) {
247 blockColumn = myLocalNodeIds[colind[j]];
248 if (isNonZero[blockColumn] ==
false) {
249 isNonZero[blockColumn] =
true;
250 nonZeroList[nNonZeros++] = blockColumn;
251 amalgCols[newNzs++] = blockColumn;
255 oldBlockRow = blockRow;
257 amalgRowPtr[blockRow + 1] = newNzs;
259 TEUCHOS_TEST_FOR_EXCEPTION((blockRow + 1 != Teuchos::as<LO>(nLocalNodes)) && (nLocalNodes != 0), MueLu::Exceptions::RuntimeError,
"VariableDofsPerNodeAmalgamation: error, computed # block rows (" << blockRow + 1 <<
") != nLocalNodes (" << nLocalNodes <<
")");
261 amalgCols.resize(amalgRowPtr[nLocalNodes]);
281 std::vector<bool> keep(amalgRowPtr[amalgRowPtr.
size() - 1],
true);
284 for (decltype(amalgRowPtr.
size()) i = 0; i < amalgRowPtr.
size() - 1; i++) {
287 for (decltype(maxDofPerNode) j = 0; j < maxDofPerNode; j++) {
288 if (dofPresent[ii++]) uniqueId[i] += temp;
299 for (decltype(amalgRowPtr.
size()) i = 0; i < amalgRowPtr.
size() - 1; i++) {
300 nodeIdSrcData[i] = uniqueId[i];
303 nodeIdTarget->doImport(*nodeIdSrc, *nodeImporter,
Xpetra::INSERT);
306 for (decltype(uniqueId.
size()) i = 0; i < uniqueId.
size(); i++) {
307 uniqueId[i] = nodeIdTargetData[i];
314 for (decltype(amalgRowPtr.
size()) i = 0; i < amalgRowPtr.
size() - 1; i++) {
315 for (
size_t j = amalgRowPtr[i]; j < amalgRowPtr[i + 1]; j++) {
316 if (uniqueId[i] != uniqueId[amalgCols[j]]) keep[j] =
false;
322 this->squeezeOutNnzs(amalgRowPtr, amalgCols, amalgVals, keep);
325 RCP<dxMV> ghostedCoords = dxMVf::Build(amalgColMap, Coords->getNumVectors());
327 TEUCHOS_TEST_FOR_EXCEPTION(amalgRowMap->getLocalNumElements() != Coords->getMap()->getLocalNumElements(), MueLu::Exceptions::RuntimeError,
"MueLu::VariableDofLaplacianFactory: the number of Coordinates and amalgamated nodes is inconsistent.");
335 Coords->replaceMap(amalgRowMap);
339 this->buildLaplacian(amalgRowPtr, amalgCols, lapVals, Coords->getNumVectors(), ghostedCoords);
342 for (decltype(amalgRowPtr.
size()) i = 0; i < amalgRowPtr.
size() - 1; i++) {
343 size_t j = amalgRowPtr[i];
344 this->MueLu_az_sort<LocalOrdinal>(&(amalgCols[j]), amalgRowPtr[i + 1] - j, NULL, &(lapVals[j]));
351 for (decltype(status.
size()) i = 0; i < status.
size(); i++) status[i] =
's';
352 for (decltype(status.
size()) i = 0; i < status.
size(); i++) {
353 if (dofPresent[i] ==
false) status[i] =
'p';
355 if (dirOrNot.
size() > 0) {
356 for (decltype(map.size()) i = 0; i < map.size(); i++) {
357 if (dirOrNot[i] ==
true) {
358 status[map[i]] =
'd';
362 Set(currentLevel,
"DofStatus", status);
368 for (
size_t i = 0; i < nLocalNodes; i++) {
369 lapCrsMat->insertLocalValues(i, amalgCols.view(amalgRowPtr[i], amalgRowPtr[i + 1] - amalgRowPtr[i]),
370 lapVals.
view(amalgRowPtr[i], amalgRowPtr[i + 1] - amalgRowPtr[i]));
372 lapCrsMat->fillComplete(amalgRowMap, amalgRowMap);
377 Set(currentLevel,
"A", lapMat);
380 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
381 void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildLaplacian(
const Teuchos::ArrayRCP<size_t>& rowPtr,
const Teuchos::ArrayRCP<LocalOrdinal>& cols,
Teuchos::ArrayRCP<Scalar>& vals,
const size_t& numdim,
const RCP<
Xpetra::MultiVector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType,
LocalOrdinal,
GlobalOrdinal,
Node> >& ghostedCoords)
const {
388 for (decltype(rowPtr.
size()) i = 0; i < rowPtr.
size() - 1; i++) {
391 for (
size_t j = rowPtr[i]; j < rowPtr[i + 1]; j++) {
392 if (cols[j] != Teuchos::as<LO>(i)) {
393 vals[j] = std::sqrt((x[i] - x[cols[j]]) * (x[i] - x[cols[j]]) +
394 (y[i] - y[cols[j]]) * (y[i] - y[cols[j]]));
411 for (decltype(rowPtr.
size()) i = 0; i < rowPtr.
size() - 1; i++) {
414 for (
size_t j = rowPtr[i]; j < rowPtr[i + 1]; j++) {
415 if (cols[j] != Teuchos::as<LO>(i)) {
416 vals[j] = std::sqrt((x[i] - x[cols[j]]) * (x[i] - x[cols[j]]) +
417 (y[i] - y[cols[j]]) * (y[i] - y[cols[j]]) +
418 (z[i] - z[cols[j]]) * (z[i] - z[cols[j]]));
435 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
443 size_t nRows = rowPtr.
size() - 1;
444 if (vals.
size() > 0) {
445 for (
size_t i = 0; i < nRows; i++) {
446 size_t newStart = count;
447 for (
size_t j = rowPtr[i]; j < rowPtr[i + 1]; j++) {
449 cols[count] = cols[j];
450 vals[count++] = vals[j];
453 rowPtr[i] = newStart;
456 for (
size_t i = 0; i < nRows; i++) {
457 size_t newStart = count;
458 for (
size_t j = rowPtr[i]; j < rowPtr[i + 1]; j++) {
459 if (keep[j] ==
true) {
460 cols[count++] = cols[j];
463 rowPtr[i] = newStart;
466 rowPtr[nRows] = count;
469 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
472 for (decltype(dofPresent.
size()) i = 0; i < dofPresent.
size(); i++)
473 if (dofPresent[i] == 1) map[count++] = Teuchos::as<LocalOrdinal>(i);
477 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
478 void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::assignGhostLocalNodeIds(
const Teuchos::RCP<const Map>& rowDofMap,
const Teuchos::RCP<const Map>& colDofMap, std::vector<LocalOrdinal>& myLocalNodeIds,
const std::vector<LocalOrdinal>& dofMap,
size_t maxDofPerNode,
size_t& nLocalNodes,
size_t& nLocalPlusGhostNodes,
Teuchos::RCP<
const Teuchos::Comm<int> > comm)
const {
479 size_t nLocalDofs = rowDofMap->getLocalNumElements();
480 size_t nLocalPlusGhostDofs = colDofMap->getLocalNumElements();
492 for (
size_t i = 0; i < localNodeIdsTemp->getLocalLength(); i++)
493 localNodeIdsTempData[i] = std::floor<LocalOrdinal>(dofMap[i] / maxDofPerNode);
496 localNodeIds->doImport(*localNodeIdsTemp, *importer,
Xpetra::INSERT);
508 for (
size_t i = 0; i < myProcTemp->getLocalLength(); i++)
509 myProcTempData[i] = Teuchos::as<LocalOrdinal>(comm->getRank());
526 std::vector<size_t> location(nLocalPlusGhostDofs - nLocalDofs + 1);
527 std::vector<size_t> tempId(nLocalPlusGhostDofs - nLocalDofs + 1);
528 std::vector<size_t> tempProc(nLocalPlusGhostDofs - nLocalDofs + 1);
530 size_t notProcessed = nLocalDofs;
531 size_t tempIndex = 0;
532 size_t first = tempIndex;
535 while (notProcessed < nLocalPlusGhostDofs) {
536 neighbor = myProcData[notProcessed];
538 location[tempIndex] = notProcessed;
539 tempId[tempIndex++] = localNodeIdsData[notProcessed];
540 myProcData[notProcessed] = -1 - neighbor;
542 for (
size_t i = notProcessed + 1; i < nLocalPlusGhostDofs; i++) {
543 if (myProcData[i] == neighbor) {
544 location[tempIndex] = i;
545 tempId[tempIndex++] = localNodeIdsData[i];
549 this->MueLu_az_sort<size_t>(&(tempId[first]), tempIndex - first, &(location[first]), NULL);
550 for (
size_t i = first; i < tempIndex; i++) tempProc[i] = neighbor;
554 while ((notProcessed < nLocalPlusGhostDofs) && (myProcData[notProcessed] < 0))
564 if (nLocalDofs > 0) nLocalNodes = localNodeIdsData[nLocalDofs - 1] + 1;
566 nLocalPlusGhostNodes = nLocalNodes;
567 if (nLocalDofs < nLocalPlusGhostDofs) nLocalPlusGhostNodes++;
574 for (
size_t i = nLocalDofs + 1; i < nLocalPlusGhostDofs; i++) {
575 size_t lagged = nLocalPlusGhostNodes - 1;
578 if ((tempId[i - nLocalDofs] != tempId[i - 1 - nLocalDofs]) ||
579 (tempProc[i - nLocalDofs] != tempProc[i - 1 - nLocalDofs]))
580 nLocalPlusGhostNodes++;
581 tempId[i - 1 - nLocalDofs] = lagged;
583 if (nLocalPlusGhostDofs > nLocalDofs)
584 tempId[nLocalPlusGhostDofs - 1 - nLocalDofs] = nLocalPlusGhostNodes - 1;
587 for (
size_t i = 0; i < nLocalDofs; i++)
588 myLocalNodeIds[i] = std::floor<LocalOrdinal>(dofMap[i] / maxDofPerNode);
591 for (
size_t i = nLocalDofs; i < nLocalPlusGhostDofs; i++)
592 myLocalNodeIds[location[i - nLocalDofs]] = tempId[i - nLocalDofs];
VariableDofLaplacianFactory()
Constructor.
void buildPaddedMap(const Teuchos::ArrayRCP< const LocalOrdinal > &dofPresent, std::vector< LocalOrdinal > &map, size_t nDofs) const
MueLu::DefaultLocalOrdinal LocalOrdinal
void buildLaplacian(const Teuchos::ArrayRCP< size_t > &rowPtr, const Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const size_t &numdim, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > &ghostedCoords) const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void DeclareInput(Level ¤tLevel) const
Input.
void assignGhostLocalNodeIds(const Teuchos::RCP< const Map > &rowDofMap, const Teuchos::RCP< const Map > &colDofMap, std::vector< LocalOrdinal > &myLocalNodeIds, const std::vector< LocalOrdinal > &dofMap, size_t maxDofPerNode, size_t &nLocalNodes, size_t &nLocalPlusGhostNodes, Teuchos::RCP< const Teuchos::Comm< int > > comm) const
static const NoFactory * get()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
void Build(Level ¤tLevel) const
Build an object with this factory.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void squeezeOutNnzs(Teuchos::ArrayRCP< size_t > &rowPtr, Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const std::vector< bool > &keep) const
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
ArrayView< T > view(size_type lowerOffset, size_type size) const