46 #ifndef MUELU_IPCFACTORY_DEF_HPP
47 #define MUELU_IPCFACTORY_DEF_HPP
59 #include "MueLu_PerfUtils.hpp"
60 #include "MueLu_Utilities.hpp"
70 #include "Intrepid2_HGRAD_LINE_C1_FEM.hpp"
71 #include "Intrepid2_HGRAD_LINE_Cn_FEM.hpp"
76 #include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
78 #include "Intrepid2_HGRAD_QUAD_Cn_FEM.hpp"
93 #define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level, ename, entry) \
95 if (level.IsRequested(ename, this) || level.GetKeepFlag(ename, this) != 0) this->Set(level, ename, entry); \
101 namespace MueLuIntrepid {
102 inline std::string
tolower(
const std::string &str) {
103 std::string data(str);
104 std::transform(data.begin(), data.end(), data.begin(),
::tolower);
109 template <
class Basis,
class LOFieldContainer,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 std::vector<std::vector<LocalOrdinal>> &seeds,
123 shards::CellTopology cellTopo = basis->getBaseCellTopology();
124 int spaceDim = cellTopo.getDimension();
126 seeds.resize(spaceDim + 1);
133 std::vector<std::set<LocalOrdinal>> seedSets(spaceDim + 1);
135 int numCells = elementToNodeMap.extent(0);
136 auto elementToNodeMap_host = Kokkos::create_mirror_view(elementToNodeMap);
137 Kokkos::deep_copy(elementToNodeMap_host, elementToNodeMap);
138 for (
int cellOrdinal = 0; cellOrdinal < numCells; cellOrdinal++) {
139 for (
int d = 0; d <= spaceDim; d++) {
140 int subcellCount = cellTopo.getSubcellCount(d);
141 for (
int subcord = 0; subcord < subcellCount; subcord++) {
142 int dofCount = basis->getDofCount(d, subcord);
143 if (dofCount == 0)
continue;
145 GO leastGlobalDofOrdinal = go_invalid;
146 LO LID_leastGlobalDofOrdinal = lo_invalid;
147 for (
int basisOrdinalOrdinal = 0; basisOrdinalOrdinal < dofCount; basisOrdinalOrdinal++) {
148 int basisOrdinal = basis->getDofOrdinal(d, subcord, basisOrdinalOrdinal);
149 int colLID = elementToNodeMap_host(cellOrdinal, basisOrdinal);
153 if (rowLID != lo_invalid) {
154 if ((leastGlobalDofOrdinal == go_invalid) || (colGID < leastGlobalDofOrdinal)) {
156 leastGlobalDofOrdinal = colGID;
157 LID_leastGlobalDofOrdinal = rowLID;
162 if (leastGlobalDofOrdinal != go_invalid) {
163 seedSets[d].insert(LID_leastGlobalDofOrdinal);
168 for (
int d = 0; d <= spaceDim; d++) {
169 seeds[d] = std::vector<LocalOrdinal>(seedSets[d].begin(), seedSets[d].end());
180 template <
class Scalar,
class KokkosExecutionSpace>
184 string myerror(
"IntrepidBasisFactory: cannot parse string name '" + name +
"'");
189 size_t pos1 = name.find_first_of(
" _");
190 if (pos1 == 0)
throw std::runtime_error(myerror);
191 string deriv =
tolower(name.substr(0, pos1));
192 if (deriv !=
"hgrad" && deriv !=
"hcurl" && deriv !=
"hdiv")
throw std::runtime_error(myerror);
196 size_t pos2 = name.find_first_of(
" _", pos1);
197 if (pos2 == 0)
throw std::runtime_error(myerror);
198 string el =
tolower(name.substr(pos1, pos2 - pos1));
199 if (el !=
"hex" && el !=
"line" && el !=
"poly" && el !=
"pyr" && el !=
"quad" && el !=
"tet" && el !=
"tri" && el !=
"wedge")
throw std::runtime_error(myerror);
203 string poly =
tolower(name.substr(pos2, 1));
204 if (poly !=
"c" && poly !=
"i")
throw std::runtime_error(myerror);
208 degree = std::stoi(name.substr(pos2, 1));
209 if (degree <= 0)
throw std::runtime_error(myerror);
212 if (deriv ==
"hgrad" && el ==
"quad" && poly ==
"c") {
214 return rcp(
new Intrepid2::Basis_HGRAD_QUAD_C1_FEM<KokkosExecutionSpace, Scalar, Scalar>());
216 return rcp(
new Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace, Scalar, Scalar>(degree, Intrepid2::POINTTYPE_EQUISPACED));
217 }
else if (deriv ==
"hgrad" && el ==
"line" && poly ==
"c") {
219 return rcp(
new Intrepid2::Basis_HGRAD_LINE_C1_FEM<KokkosExecutionSpace, Scalar, Scalar>());
221 return rcp(
new Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace, Scalar, Scalar>(degree, Intrepid2::POINTTYPE_EQUISPACED));
225 throw std::runtime_error(myerror);
236 template <
class Scalar,
class KokkosDeviceType>
238 std::vector<size_t> &lo_node_in_hi,
239 Kokkos::DynRankView<Scalar, KokkosDeviceType> &hi_DofCoords) {
240 typedef typename KokkosDeviceType::execution_space KokkosExecutionSpace;
242 size_t degree = hi_basis->getDegree();
243 lo_node_in_hi.resize(0);
245 if (!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace, Scalar, Scalar>>(hi_basis).
is_null()) {
247 lo_node_in_hi.insert(lo_node_in_hi.end(), {0, degree, (degree + 1) * (degree + 1) - 1, degree * (degree + 1)});
248 }
else if (!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace, Scalar, Scalar>>(hi_basis).is_null()) {
250 lo_node_in_hi.insert(lo_node_in_hi.end(), {0, degree});
252 throw std::runtime_error(
"IntrepidPCoarsenFactory: Unknown element type");
255 Kokkos::resize(hi_DofCoords, hi_basis->getCardinality(), hi_basis->getBaseCellTopology().getDimension());
256 hi_basis->getDofCoords(hi_DofCoords);
267 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class LOFieldContainer>
270 LOFieldContainer &lo_elemToHiRepresentativeNode) {
276 size_t numElem = hi_elemToNode.extent(0);
277 size_t lo_nperel = candidates.size();
278 Kokkos::resize(lo_elemToHiRepresentativeNode, numElem, lo_nperel);
280 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
281 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
282 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
283 for (
size_t i = 0; i < numElem; i++)
284 for (
size_t j = 0; j < lo_nperel; j++) {
285 if (candidates[j].size() == 1)
286 lo_elemToHiRepresentativeNode_host(i, j) = hi_elemToNode_host(i, candidates[j][0]);
289 std::vector<GO> GID(candidates[j].size());
290 for (
size_t k = 0; k < (size_t)candidates[j].size(); k++)
291 GID[k] = hi_columnMap->getGlobalElement(hi_elemToNode_host(i, candidates[j][k]));
294 size_t which = std::distance(GID.begin(), std::min_element(GID.begin(), GID.end()));
297 lo_elemToHiRepresentativeNode_host(i, j) = hi_elemToNode_host(i, candidates[j][which]);
300 Kokkos::deep_copy(lo_elemToHiRepresentativeNode, lo_elemToHiRepresentativeNode_host);
313 template <
class LocalOrdinal,
class LOFieldContainer>
315 const std::vector<bool> &hi_nodeIsOwned,
316 const LOFieldContainer &lo_elemToHiRepresentativeNode,
317 LOFieldContainer &lo_elemToNode,
318 std::vector<bool> &lo_nodeIsOwned,
319 std::vector<LocalOrdinal> &hi_to_lo_map,
320 int &lo_numOwnedNodes) {
324 size_t numElem = hi_elemToNode.extent(0);
325 size_t hi_numNodes = hi_nodeIsOwned.size();
326 size_t lo_nperel = lo_elemToHiRepresentativeNode.extent(1);
327 Kokkos::resize(lo_elemToNode, numElem, lo_nperel);
330 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
331 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
332 std::vector<bool> is_low_order(hi_numNodes,
false);
333 for (
size_t i = 0; i < numElem; i++)
334 for (
size_t j = 0; j < lo_nperel; j++) {
335 LO
id = lo_elemToHiRepresentativeNode_host(i, j);
336 is_low_order[id] =
true;
340 lo_numOwnedNodes = 0;
341 size_t lo_numNodes = 0;
344 for (
size_t i = 0; i < hi_numNodes; i++)
345 if (is_low_order[i]) {
346 hi_to_lo_map[i] = lo_numNodes;
348 if (hi_nodeIsOwned[i]) lo_numOwnedNodes++;
352 lo_nodeIsOwned.resize(lo_numNodes,
false);
353 for (
size_t i = 0; i < hi_numNodes; i++) {
354 if (is_low_order[i] && hi_nodeIsOwned[i])
355 lo_nodeIsOwned[hi_to_lo_map[i]] =
true;
359 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
360 for (
size_t i = 0; i < numElem; i++)
361 for (
size_t j = 0; j < lo_nperel; j++)
362 lo_elemToNode_host(i, j) = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i, j)];
366 bool map_ordering_test_passed =
true;
367 for (
size_t i = 0; i < lo_numNodes - 1; i++)
368 if (!lo_nodeIsOwned[i] && lo_nodeIsOwned[i + 1])
369 map_ordering_test_passed =
false;
371 if (!map_ordering_test_passed)
372 throw std::runtime_error(
"MueLu::MueLuIntrepid::BuildLoElemToNodeViaRepresentatives failed map ordering test");
373 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
387 template <
class LocalOrdinal,
class LOFieldContainer>
389 const std::vector<bool> &hi_nodeIsOwned,
390 const std::vector<size_t> &lo_node_in_hi,
392 LOFieldContainer &lo_elemToNode,
393 std::vector<bool> &lo_nodeIsOwned,
394 std::vector<LocalOrdinal> &hi_to_lo_map,
395 int &lo_numOwnedNodes) {
401 size_t numElem = hi_elemToNode.extent(0);
402 size_t hi_numNodes = hi_nodeIsOwned.size();
404 size_t lo_nperel = lo_node_in_hi.size();
405 Kokkos::resize(lo_elemToNode, numElem, lo_nperel);
408 std::vector<bool> is_low_order(hi_numNodes,
false);
409 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
410 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
411 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
412 for (
size_t i = 0; i < numElem; i++)
413 for (
size_t j = 0; j < lo_nperel; j++) {
414 LO lid = hi_elemToNode_host(i, lo_node_in_hi[j]);
417 if (hi_isDirichlet[lid])
418 lo_elemToNode_host(i, j) = LOINVALID;
420 lo_elemToNode_host(i, j) = lid;
421 is_low_order[hi_elemToNode_host(i, lo_node_in_hi[j])] =
true;
426 lo_numOwnedNodes = 0;
427 size_t lo_numNodes = 0;
430 for (
size_t i = 0; i < hi_numNodes; i++)
431 if (is_low_order[i]) {
432 hi_to_lo_map[i] = lo_numNodes;
434 if (hi_nodeIsOwned[i]) lo_numOwnedNodes++;
438 lo_nodeIsOwned.resize(lo_numNodes,
false);
439 for (
size_t i = 0; i < hi_numNodes; i++) {
440 if (is_low_order[i] && hi_nodeIsOwned[i])
441 lo_nodeIsOwned[hi_to_lo_map[i]] =
true;
445 for (
size_t i = 0; i < numElem; i++)
446 for (
size_t j = 0; j < lo_nperel; j++) {
447 if (lo_elemToNode_host(i, j) != LOINVALID)
448 lo_elemToNode_host(i, j) = hi_to_lo_map[lo_elemToNode_host(i, j)];
450 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
454 bool map_ordering_test_passed =
true;
455 for (
size_t i = 0; i < lo_numNodes - 1; i++)
456 if (!lo_nodeIsOwned[i] && lo_nodeIsOwned[i + 1])
457 map_ordering_test_passed =
false;
459 if (!map_ordering_test_passed)
460 throw std::runtime_error(
"MueLu::MueLuIntrepid::BuildLoElemToNode failed map ordering test");
472 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
473 void GenerateColMapFromImport(
const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> &hi_importer,
const std::vector<LocalOrdinal> &hi_to_lo_map,
const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> &lo_domainMap,
const size_t &lo_columnMapLength,
RCP<
const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> &lo_columnMap) {
493 for (
size_t i = 0; i < hi_domainMap->getLocalNumElements(); i++) {
494 if (hi_to_lo_map[i] != lo_invalid)
497 dvec_data[i] = go_invalid;
506 Array<GO> lo_col_data(lo_columnMapLength);
509 for (
size_t i = 0, idx = 0; i < hi_columnMap->getLocalNumElements(); i++) {
510 if (hi_to_lo_map[i] != lo_invalid) {
511 lo_col_data[idx] = cvec_data[i];
528 template <
class Basis,
class SCFieldContainer>
529 void GenerateRepresentativeBasisNodes(
const Basis &basis,
const SCFieldContainer &ReferenceNodeLocations,
const double threshold, std::vector<std::vector<size_t>> &representative_node_candidates) {
530 typedef SCFieldContainer FC;
531 typedef typename FC::data_type
SC;
534 size_t numFieldsHi = ReferenceNodeLocations.extent(0);
536 size_t numFieldsLo = basis.getCardinality();
538 FC LoValues(
"LoValues", numFieldsLo, numFieldsHi);
540 basis.getValues(LoValues, ReferenceNodeLocations, Intrepid2::OPERATOR_VALUE);
545 printf(
"** LoValues[%d,%d] **\n",(
int)numFieldsLo,(
int)numFieldsHi);
546 for(
size_t i=0; i<numFieldsLo; i++) {
547 for(
size_t j=0; j<numFieldsHi; j++)
548 printf(
"%6.4e ",LoValues(i,j));
551 printf(
"**************\n");fflush(stdout);
554 representative_node_candidates.resize(numFieldsLo);
555 auto LoValues_host = Kokkos::create_mirror_view(LoValues);
556 Kokkos::deep_copy(LoValues_host, LoValues);
557 for (
size_t i = 0; i < numFieldsLo; i++) {
560 for (
size_t j = 0; j < numFieldsHi; j++)
564 for (
size_t j = 0; j < numFieldsHi; j++) {
566 representative_node_candidates[i].push_back(j);
571 for (
size_t i = 0; i < numFieldsLo; i++)
572 if (!representative_node_candidates[i].size())
573 throw std::runtime_error(
"ERROR: GenerateRepresentativeBasisNodes: No candidates found!");
581 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
583 const std::vector<bool> &hi_nodeIsOwned,
585 const std::vector<size_t> &lo_node_in_hi,
586 const Basis &lo_basis,
587 const std::vector<LocalOrdinal> &hi_to_lo_map,
594 size_t numFieldsHi = hi_elemToNode.extent(1);
595 size_t numFieldsLo = lo_basis.getCardinality();
597 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs", numFieldsLo, numFieldsHi);
598 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
599 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
600 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
608 P =
rcp(
new CrsMatrixWrap(hi_map, lo_colMap, numFieldsHi));
609 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
612 size_t Nelem = hi_elemToNode.extent(0);
613 std::vector<bool> touched(hi_map->getLocalNumElements(),
false);
616 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
617 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
618 for (
size_t i = 0; i < Nelem; i++) {
619 for (
size_t j = 0; j < numFieldsHi; j++) {
620 LO row_lid = hi_elemToNode_host(i, j);
621 GO row_gid = hi_map->getGlobalElement(row_lid);
622 if (hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
623 for (
size_t k = 0; k < numFieldsLo; k++) {
625 LO col_lid = hi_to_lo_map[hi_elemToNode_host(i, lo_node_in_hi[k])];
626 if (col_lid == LOINVALID)
continue;
628 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
629 val[0] = LoValues_at_HiDofs_host(k, j);
633 P->insertGlobalValues(row_gid, col_gid(), val());
635 touched[row_lid] =
true;
639 P->fillComplete(lo_domainMap, hi_map);
643 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
645 const std::vector<bool> &hi_nodeIsOwned,
648 const Basis &lo_basis,
649 const std::vector<LocalOrdinal> &hi_to_lo_map,
656 size_t numFieldsHi = hi_elemToNode.extent(1);
657 size_t numFieldsLo = lo_basis.getCardinality();
658 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs", numFieldsLo, numFieldsHi);
659 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
660 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
661 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
662 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
663 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
664 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
665 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
673 P =
rcp(
new CrsMatrixWrap(hi_map, lo_colMap, numFieldsHi));
674 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
677 size_t Nelem = hi_elemToNode.extent(0);
678 std::vector<bool> touched(hi_map->getLocalNumElements(),
false);
681 for (
size_t i = 0; i < Nelem; i++) {
682 for (
size_t j = 0; j < numFieldsHi; j++) {
683 LO row_lid = hi_elemToNode_host(i, j);
684 GO row_gid = hi_map->getGlobalElement(row_lid);
685 if (hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
686 for (
size_t k = 0; k < numFieldsLo; k++) {
688 LO col_lid = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i, k)];
689 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
690 val[0] = LoValues_at_HiDofs_host(k, j);
694 P->insertGlobalValues(row_gid, col_gid(), val());
696 touched[row_lid] =
true;
700 P->fillComplete(lo_domainMap, hi_map);
704 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
708 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
711 #undef SET_VALID_ENTRY
713 validParamList->set<
RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
715 validParamList->set<
RCP<const FactoryBase>>(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
716 validParamList->set<
RCP<const FactoryBase>>(
"pcoarsen: element to node map", Teuchos::null,
"Generating factory of the element to node map");
717 return validParamList;
721 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
723 Input(fineLevel,
"A");
724 Input(fineLevel,
"pcoarsen: element to node map");
725 Input(fineLevel,
"Nullspace");
729 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
731 return BuildP(fineLevel, coarseLevel);
735 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
739 const std::string prefix =
"MueLu::IntrepidPCoarsenFactory(" + levelIDs +
"): ";
742 typedef Kokkos::DynRankView<LocalOrdinal, typename Node::device_type> FCi;
743 typedef Kokkos::DynRankView<double, typename Node::device_type> FC;
747 RCP<MultiVector> fineNullspace = Get<RCP<MultiVector>>(fineLevel,
"Nullspace");
750 if (restrictionMode_) {
756 std::vector<LocalOrdinal> A_dirichletRows;
764 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
765 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
769 if (APparams->isParameter(
"graph"))
779 int lo_degree, hi_degree;
780 RCP<Basis> hi_basis = MueLuIntrepid::BasisFactory<double, typename Node::device_type::execution_space>(pL.
get<std::string>(
"pcoarsen: hi basis"), hi_degree);
781 RCP<Basis> lo_basis = MueLuIntrepid::BasisFactory<double, typename Node::device_type::execution_space>(pL.
get<std::string>(
"pcoarsen: lo basis"), lo_degree);
784 GetOStream(
Statistics1) <<
"P-Coarsening from basis " << pL.
get<std::string>(
"pcoarsen: hi basis") <<
" to " << pL.
get<std::string>(
"pcoarsen: lo basis") << std::endl;
788 const Teuchos::RCP<FCi> Pn_elemToNode = Get<Teuchos::RCP<FCi>>(fineLevel,
"pcoarsen: element to node map");
798 int NumProc = rowMap->getComm()->getSize();
799 assert(rowMap->isSameAs(*domainMap));
800 std::vector<bool> Pn_nodeIsOwned(colMap->getLocalNumElements(),
false);
801 LO num_owned_rows = 0;
802 for (
size_t i = 0; i < rowMap->getLocalNumElements(); i++) {
803 if (rowMap->getGlobalElement(i) == colMap->getGlobalElement(i)) {
804 Pn_nodeIsOwned[i] =
true;
813 std::vector<bool> P1_nodeIsOwned;
814 int P1_numOwnedNodes;
815 std::vector<LO> hi_to_lo_map;
818 std::vector<size_t> lo_node_in_hi;
821 FCi lo_elemToHiRepresentativeNode;
828 printf(
"[%d] isDirichletRow = ",A->getRowMap()->getComm()->getRank());
829 for(
size_t i=0;i<hi_isDirichletRow->getMap()->getLocalNumElements(); i++)
830 printf(
"%d ",hi_isDirichletRow->getData(0)[i]);
832 printf(
"[%d] isDirichletCol = ",A->getRowMap()->getComm()->getRank());
833 for(
size_t i=0;i<hi_isDirichletCol->getMap()->getLocalNumElements(); i++)
834 printf(
"%d ",hi_isDirichletCol->getData(0)[i]);
840 if (lo_degree == 1) {
845 MueLuIntrepid::BuildLoElemToNode(*Pn_elemToNode, Pn_nodeIsOwned, lo_node_in_hi, hi_isDirichletCol->getData(0), *P1_elemToNode, P1_nodeIsOwned, hi_to_lo_map, P1_numOwnedNodes);
846 assert(hi_to_lo_map.size() == colMap->getLocalNumElements());
849 double threshold = 1e-10;
850 std::vector<std::vector<size_t>> candidates;
851 Kokkos::resize(hi_DofCoords, hi_basis->getCardinality(), hi_basis->getBaseCellTopology().getDimension());
852 hi_basis->getDofCoords(hi_DofCoords);
854 MueLu::MueLuIntrepid::GenerateRepresentativeBasisNodes<Basis, FC>(*lo_basis, hi_DofCoords, threshold, candidates);
871 P1_colMap = P1_domainMap;
873 MueLuIntrepid::GenerateColMapFromImport<LO, GO, NO>(*Acrs.getCrsGraph()->getImporter(), hi_to_lo_map, *P1_domainMap, P1_nodeIsOwned.size(), P1_colMap);
878 GenerateLinearCoarsening_pn_kirby_to_p1(*Pn_elemToNode, Pn_nodeIsOwned, hi_DofCoords, lo_node_in_hi, *lo_basis, hi_to_lo_map, P1_colMap, P1_domainMap, A->getRowMap(), finalP);
880 GenerateLinearCoarsening_pn_kirby_to_pm(*Pn_elemToNode, Pn_nodeIsOwned, hi_DofCoords, lo_elemToHiRepresentativeNode, *lo_basis, hi_to_lo_map, P1_colMap, P1_domainMap, A->getRowMap(), finalP);
888 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P1_domainMap, fineNullspace->getNumVectors());
890 Set(coarseLevel,
"Nullspace", coarseNullspace);
893 if (!restrictionMode_) {
895 Set(coarseLevel,
"P", finalP);
897 APparams->set(
"graph", finalP);
902 params->set(
"printLoadBalancingInfo",
true);
903 params->set(
"printCommInfo",
true);
914 Set(coarseLevel,
"R", R);
918 params->set(
"printLoadBalancingInfo",
true);
919 params->set(
"printCommInfo",
true);
928 #endif // MUELU_IPCFACTORY_DEF_HPP
void GenerateLoNodeInHiViaGIDs(const std::vector< std::vector< size_t > > &candidates, const LOFieldContainer &hi_elemToNode, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &hi_columnMap, LOFieldContainer &lo_elemToHiRepresentativeNode)
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const =0
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Teuchos::RCP< Intrepid2::Basis< KokkosExecutionSpace, Scalar, Scalar > > BasisFactory(const std::string &name, int °ree)
bool is_null(const boost::shared_ptr< T > &p)
void FindGeometricSeedOrdinals(Teuchos::RCP< Basis > basis, const LOFieldContainer &elementToNodeMap, std::vector< std::vector< LocalOrdinal > > &seeds, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &columnMap)
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node >> &isDirichletCol)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
static magnitudeType eps()
T & get(const std::string &name, T def_value)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Timer to be used in factories. Similar to Monitor but with additional timers.
void IntrepidGetP1NodeInHi(const Teuchos::RCP< Intrepid2::Basis< typename KokkosDeviceType::execution_space, Scalar, Scalar >> &hi_basis, std::vector< size_t > &lo_node_in_hi, Kokkos::DynRankView< Scalar, KokkosDeviceType > &hi_DofCoords)
One-liner description of what is happening.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
virtual LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const =0
std::string tolower(const std::string &str)
virtual GlobalOrdinal getIndexBase() const =0
Kokkos::DynRankView< LocalOrdinal, typename Node::device_type > LOFieldContainer
Print even more statistics.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getTargetMap() const =0
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void BuildLoElemToNode(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const std::vector< size_t > &lo_node_in_hi, const Teuchos::ArrayRCP< const int > &hi_isDirichlet, LOFieldContainer &lo_elemToNode, std::vector< bool > &lo_nodeIsOwned, std::vector< LocalOrdinal > &hi_to_lo_map, int &lo_numOwnedNodes)
Kokkos::DynRankView< double, typename Node::device_type > SCFieldContainer
void BuildLoElemToNodeViaRepresentatives(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const LOFieldContainer &lo_elemToHiRepresentativeNode, LOFieldContainer &lo_elemToNode, std::vector< bool > &lo_nodeIsOwned, std::vector< LocalOrdinal > &hi_to_lo_map, int &lo_numOwnedNodes)
#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level, ename, entry)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void GenerateLinearCoarsening_pn_kirby_to_p1(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const SCFieldContainer &hi_DofCoords, const std::vector< size_t > &lo_node_in_hi, const Basis &lo_Basis, const std::vector< LocalOrdinal > &hi_to_lo_map, const Teuchos::RCP< const Map > &lo_colMap, const Teuchos::RCP< const Map > &lo_domainMap, const Teuchos::RCP< const Map > &hi_map, Teuchos::RCP< Matrix > &P) const
virtual GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const =0
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void GenerateColMapFromImport(const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &hi_importer, const std::vector< LocalOrdinal > &hi_to_lo_map, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &lo_domainMap, const size_t &lo_columnMapLength, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &lo_columnMap)
Intrepid2::Basis< typename Node::device_type::execution_space, double, double > Basis
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void GenerateRepresentativeBasisNodes(const Basis &basis, const SCFieldContainer &ReferenceNodeLocations, const double threshold, std::vector< std::vector< size_t > > &representative_node_candidates)
int GetLevelID() const
Return level number.
void GenerateLinearCoarsening_pn_kirby_to_pm(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const SCFieldContainer &hi_DofCoords, const LOFieldContainer &lo_elemToHiRepresentativeNode, const Basis &lo_basis, const std::vector< LocalOrdinal > &hi_to_lo_map, const Teuchos::RCP< const Map > &lo_colMap, const Teuchos::RCP< const Map > &lo_domainMap, const Teuchos::RCP< const Map > &hi_map, Teuchos::RCP< Matrix > &P) const
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getSourceMap() const =0
virtual UnderlyingLib lib() const =0
#define SET_VALID_ENTRY(name)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.