46#ifndef MUELU_IPCFACTORY_DEF_HPP
47#define MUELU_IPCFACTORY_DEF_HPP
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_IO.hpp>
60#include "MueLu_PerfUtils.hpp"
62#include "MueLu_Utilities.hpp"
64#include "Teuchos_ScalarTraits.hpp"
72#include "Intrepid2_HGRAD_LINE_C1_FEM.hpp"
73#include "Intrepid2_HGRAD_LINE_Cn_FEM.hpp"
78#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
80#include "Intrepid2_HGRAD_QUAD_Cn_FEM.hpp"
95#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level,ename,entry) \
96 {if (level.IsRequested(ename,this) || level.GetKeepFlag(ename,this) != 0) this->Set(level,ename,entry);}
104namespace MueLuIntrepid {
105inline std::string
tolower(
const std::string & str) {
106 std::string data(str);
107 std::transform(data.begin(), data.end(), data.begin(),
::tolower);
113 template<
class Basis,
class LOFieldContainer,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
115 std::vector<std::vector<LocalOrdinal> > &seeds,
116 const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap,
117 const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &columnMap)
129 shards::CellTopology cellTopo = basis->getBaseCellTopology();
130 int spaceDim = cellTopo.getDimension();
132 seeds.resize(spaceDim + 1);
136 LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
137 GlobalOrdinal go_invalid = Teuchos::OrdinalTraits<GO>::invalid();
140 std::vector<std::set<LocalOrdinal>> seedSets(spaceDim+1);
142 int numCells = elementToNodeMap.extent(0);
143 auto elementToNodeMap_host = Kokkos::create_mirror_view(elementToNodeMap);
144 Kokkos::deep_copy(elementToNodeMap_host,elementToNodeMap);
145 for (
int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
147 for (
int d=0; d<=spaceDim; d++)
149 int subcellCount = cellTopo.getSubcellCount(d);
150 for (
int subcord=0; subcord<subcellCount; subcord++)
152 int dofCount = basis->getDofCount(d,subcord);
153 if (dofCount == 0)
continue;
155 GO leastGlobalDofOrdinal = go_invalid;
156 LO LID_leastGlobalDofOrdinal = lo_invalid;
157 for (
int basisOrdinalOrdinal=0; basisOrdinalOrdinal<dofCount; basisOrdinalOrdinal++)
159 int basisOrdinal = basis->getDofOrdinal(d,subcord,basisOrdinalOrdinal);
160 int colLID = elementToNodeMap_host(cellOrdinal,basisOrdinal);
161 if (colLID != Teuchos::OrdinalTraits<LO>::invalid())
165 if (rowLID != lo_invalid)
167 if ((leastGlobalDofOrdinal == go_invalid) || (colGID < leastGlobalDofOrdinal))
170 leastGlobalDofOrdinal = colGID;
171 LID_leastGlobalDofOrdinal = rowLID;
176 if (leastGlobalDofOrdinal != go_invalid)
178 seedSets[d].insert(LID_leastGlobalDofOrdinal);
183 for (
int d=0; d<=spaceDim;d++)
185 seeds[d] = std::vector<LocalOrdinal>(seedSets[d].begin(),seedSets[d].end());
196template<
class Scalar,
class KokkosExecutionSpace>
197Teuchos::RCP< Intrepid2::Basis<KokkosExecutionSpace,Scalar,Scalar> >
BasisFactory(
const std::string & name,
int & degree)
201 string myerror(
"IntrepidBasisFactory: cannot parse string name '"+name+
"'");
206 size_t pos1 = name.find_first_of(
" _");
207 if(pos1==0)
throw std::runtime_error(myerror);
208 string deriv =
tolower(name.substr(0,pos1));
209 if(deriv!=
"hgrad" && deriv!=
"hcurl" && deriv!=
"hdiv")
throw std::runtime_error(myerror);
213 size_t pos2 = name.find_first_of(
" _",pos1);
214 if(pos2==0)
throw std::runtime_error(myerror);
215 string el =
tolower(name.substr(pos1,pos2-pos1));
216 if(el!=
"hex" && el!=
"line" && el!=
"poly" && el!=
"pyr" && el!=
"quad" && el!=
"tet" && el!=
"tri" && el!=
"wedge")
throw std::runtime_error(myerror);
220 string poly =
tolower(name.substr(pos2,1));
221 if(poly!=
"c" && poly!=
"i")
throw std::runtime_error(myerror);
225 degree=std::stoi(name.substr(pos2,1));
226 if(degree<=0)
throw std::runtime_error(myerror);
229 if(deriv==
"hgrad" && el==
"quad" && poly==
"c"){
230 if(degree==1)
return rcp(
new Intrepid2::Basis_HGRAD_QUAD_C1_FEM<KokkosExecutionSpace,Scalar,Scalar>());
231 else return rcp(
new Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar>(degree,Intrepid2::POINTTYPE_EQUISPACED));
233 else if(deriv==
"hgrad" && el==
"line" && poly==
"c"){
234 if(degree==1)
return rcp(
new Intrepid2::Basis_HGRAD_LINE_C1_FEM<KokkosExecutionSpace,Scalar,Scalar>());
235 else return rcp(
new Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar>(degree,Intrepid2::POINTTYPE_EQUISPACED));
239 throw std::runtime_error(myerror);
240 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
250template<
class Scalar,
class KokkosDeviceType>
251void IntrepidGetP1NodeInHi(
const Teuchos::RCP<Intrepid2::Basis<typename KokkosDeviceType::execution_space,Scalar,Scalar> >&hi_basis,
252 std::vector<size_t> & lo_node_in_hi,
253 Kokkos::DynRankView<Scalar,KokkosDeviceType> & hi_DofCoords) {
255 typedef typename KokkosDeviceType::execution_space KokkosExecutionSpace;
257 size_t degree = hi_basis->getDegree();
258 lo_node_in_hi.resize(0);
260 if(!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar> >(hi_basis).is_null()) {
262 lo_node_in_hi.insert(lo_node_in_hi.end(),{0,degree, (degree+1)*(degree+1)-1, degree*(degree+1)});
264 else if(!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar> >(hi_basis).is_null()) {
266 lo_node_in_hi.insert(lo_node_in_hi.end(),{0,degree});
269 throw std::runtime_error(
"IntrepidPCoarsenFactory: Unknown element type");
272 Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
273 hi_basis->getDofCoords(hi_DofCoords);
285template<
class LocalOrdinal,
class GlobalOrdinal,
class Node,
class LOFieldContainer>
287 RCP<
const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > & hi_columnMap,
288 LOFieldContainer & lo_elemToHiRepresentativeNode) {
294 size_t numElem = hi_elemToNode.extent(0);
295 size_t lo_nperel = candidates.size();
296 Kokkos::resize(lo_elemToHiRepresentativeNode,numElem, lo_nperel);
298 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
299 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
300 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
301 for(
size_t i=0; i<numElem; i++)
302 for(
size_t j=0; j<lo_nperel; j++) {
303 if(candidates[j].size() == 1)
304 lo_elemToHiRepresentativeNode_host(i,j) = hi_elemToNode_host(i,candidates[j][0]);
307 std::vector<GO> GID(candidates[j].size());
308 for(
size_t k=0; k<(size_t)candidates[j].size(); k++)
309 GID[k] = hi_columnMap->getGlobalElement(hi_elemToNode_host(i,candidates[j][k]));
312 size_t which = std::distance(GID.begin(),std::min_element(GID.begin(),GID.end()));
315 lo_elemToHiRepresentativeNode_host(i,j) = hi_elemToNode_host(i,candidates[j][which]);
318 Kokkos::deep_copy(lo_elemToHiRepresentativeNode, lo_elemToHiRepresentativeNode_host);
331template <
class LocalOrdinal,
class LOFieldContainer>
333 const std::vector<bool> & hi_nodeIsOwned,
334 const LOFieldContainer & lo_elemToHiRepresentativeNode,
335 LOFieldContainer & lo_elemToNode,
336 std::vector<bool> & lo_nodeIsOwned,
337 std::vector<LocalOrdinal> & hi_to_lo_map,
338 int & lo_numOwnedNodes) {
342 size_t numElem = hi_elemToNode.extent(0);
343 size_t hi_numNodes = hi_nodeIsOwned.size();
344 size_t lo_nperel = lo_elemToHiRepresentativeNode.extent(1);
345 Kokkos::resize(lo_elemToNode,numElem, lo_nperel);
348 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
349 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
350 std::vector<bool> is_low_order(hi_numNodes,
false);
351 for(
size_t i=0; i<numElem; i++)
352 for(
size_t j=0; j<lo_nperel; j++) {
353 LO
id = lo_elemToHiRepresentativeNode_host(i,j);
354 is_low_order[id] =
true;
359 size_t lo_numNodes=0;
360 hi_to_lo_map.resize(hi_numNodes,Teuchos::OrdinalTraits<LO>::invalid());
362 for(
size_t i=0; i<hi_numNodes; i++)
363 if(is_low_order[i]) {
364 hi_to_lo_map[i] = lo_numNodes;
366 if(hi_nodeIsOwned[i]) lo_numOwnedNodes++;
370 lo_nodeIsOwned.resize(lo_numNodes,
false);
371 for(
size_t i=0; i<hi_numNodes; i++) {
372 if(is_low_order[i] && hi_nodeIsOwned[i])
373 lo_nodeIsOwned[hi_to_lo_map[i]]=
true;
377 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
378 for(
size_t i=0; i<numElem; i++)
379 for(
size_t j=0; j<lo_nperel; j++)
380 lo_elemToNode_host(i,j) = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,j)];
385 bool map_ordering_test_passed=
true;
386 for(
size_t i=0; i<lo_numNodes-1; i++)
387 if(!lo_nodeIsOwned[i] && lo_nodeIsOwned[i+1])
388 map_ordering_test_passed=
false;
390 if(!map_ordering_test_passed)
391 throw std::runtime_error(
"MueLu::MueLuIntrepid::BuildLoElemToNodeViaRepresentatives failed map ordering test");
392 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
408template <
class LocalOrdinal,
class LOFieldContainer>
410 const std::vector<bool> & hi_nodeIsOwned,
411 const std::vector<size_t> & lo_node_in_hi,
412 const Teuchos::ArrayRCP<const int> & hi_isDirichlet,
413 LOFieldContainer & lo_elemToNode,
414 std::vector<bool> & lo_nodeIsOwned,
415 std::vector<LocalOrdinal> & hi_to_lo_map,
416 int & lo_numOwnedNodes) {
419 LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
422 size_t numElem = hi_elemToNode.extent(0);
423 size_t hi_numNodes = hi_nodeIsOwned.size();
425 size_t lo_nperel = lo_node_in_hi.size();
426 Kokkos::resize(lo_elemToNode,numElem, lo_nperel);
429 std::vector<bool> is_low_order(hi_numNodes,
false);
430 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
431 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
432 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
433 for(
size_t i=0; i<numElem; i++)
434 for(
size_t j=0; j<lo_nperel; j++) {
435 LO lid = hi_elemToNode_host(i,lo_node_in_hi[j]);
438 if(hi_isDirichlet[lid])
439 lo_elemToNode_host(i,j) = LOINVALID;
441 lo_elemToNode_host(i,j) = lid;
442 is_low_order[hi_elemToNode_host(i,lo_node_in_hi[j])] =
true;
448 size_t lo_numNodes=0;
449 hi_to_lo_map.resize(hi_numNodes,Teuchos::OrdinalTraits<LO>::invalid());
451 for(
size_t i=0; i<hi_numNodes; i++)
452 if(is_low_order[i]) {
453 hi_to_lo_map[i] = lo_numNodes;
455 if(hi_nodeIsOwned[i]) lo_numOwnedNodes++;
459 lo_nodeIsOwned.resize(lo_numNodes,
false);
460 for(
size_t i=0; i<hi_numNodes; i++) {
461 if(is_low_order[i] && hi_nodeIsOwned[i])
462 lo_nodeIsOwned[hi_to_lo_map[i]]=
true;
466 for(
size_t i=0; i<numElem; i++)
467 for(
size_t j=0; j<lo_nperel; j++) {
468 if(lo_elemToNode_host(i,j) != LOINVALID)
469 lo_elemToNode_host(i,j) = hi_to_lo_map[lo_elemToNode_host(i,j)];
471 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
475 bool map_ordering_test_passed=
true;
476 for(
size_t i=0; i<lo_numNodes-1; i++)
477 if(!lo_nodeIsOwned[i] && lo_nodeIsOwned[i+1])
478 map_ordering_test_passed=
false;
480 if(!map_ordering_test_passed)
481 throw std::runtime_error(
"MueLu::MueLuIntrepid::BuildLoElemToNode failed map ordering test");
494 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
495 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) {
499 typedef Xpetra::Map<LO,GO,NO> Map;
500 typedef Xpetra::Vector<GO,LO,GO,NO> GOVector;
502 GO go_invalid = Teuchos::OrdinalTraits<GO>::invalid();
503 LO lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
505 RCP<const Map> hi_domainMap = hi_importer.getSourceMap();
506 RCP<const Map> hi_columnMap = hi_importer.getTargetMap();
512 RCP<GOVector> dvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_domainMap);
514 ArrayRCP<GO> dvec_data = dvec->getDataNonConst(0);
515 for(
size_t i=0; i<hi_domainMap->getLocalNumElements(); i++) {
516 if(hi_to_lo_map[i]!=lo_invalid) dvec_data[i] = lo_domainMap.getGlobalElement(hi_to_lo_map[i]);
517 else dvec_data[i] = go_invalid;
522 RCP<GOVector> cvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_columnMap,
true);
523 cvec->doImport(*dvec,hi_importer,Xpetra::ADD);
527 Array<GO> lo_col_data(lo_columnMapLength);
529 ArrayRCP<GO> cvec_data = cvec->getDataNonConst(0);
530 for(
size_t i=0,idx=0; i<hi_columnMap->getLocalNumElements(); i++) {
531 if(hi_to_lo_map[i]!=lo_invalid) {
532 lo_col_data[idx] = cvec_data[i];
538 lo_columnMap = Xpetra::MapFactory<LO,GO,NO>::Build(lo_domainMap.lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),lo_col_data(),lo_domainMap.getIndexBase(),lo_domainMap.getComm());
549template<
class Basis,
class SCFieldContainer>
550void GenerateRepresentativeBasisNodes(
const Basis & basis,
const SCFieldContainer & ReferenceNodeLocations,
const double threshold,std::vector<std::vector<size_t> > & representative_node_candidates) {
551 typedef SCFieldContainer FC;
552 typedef typename FC::data_type SC;
555 size_t numFieldsHi = ReferenceNodeLocations.extent(0);
557 size_t numFieldsLo = basis.getCardinality();
559 FC LoValues(
"LoValues",numFieldsLo,numFieldsHi);
561 basis.getValues(LoValues, ReferenceNodeLocations , Intrepid2::OPERATOR_VALUE);
566 printf(
"** LoValues[%d,%d] **\n",(
int)numFieldsLo,(
int)numFieldsHi);
567 for(
size_t i=0; i<numFieldsLo; i++) {
568 for(
size_t j=0; j<numFieldsHi; j++)
569 printf(
"%6.4e ",LoValues(i,j));
572 printf(
"**************\n");fflush(stdout);
575 representative_node_candidates.resize(numFieldsLo);
576 auto LoValues_host = Kokkos::create_mirror_view(LoValues);
577 Kokkos::deep_copy(LoValues_host, LoValues);
578 for(
size_t i=0; i<numFieldsLo; i++) {
580 typename Teuchos::ScalarTraits<SC>::magnitudeType vmax = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero();
581 for(
size_t j=0; j<numFieldsHi; j++)
582 vmax = std::max(vmax,Teuchos::ScalarTraits<SC>::magnitude(LoValues_host(i,j)));
585 for(
size_t j=0; j<numFieldsHi; j++) {
586 if(Teuchos::ScalarTraits<SC>::magnitude(vmax - LoValues_host(i,j)) < threshold*vmax)
587 representative_node_candidates[i].push_back(j);
592 for(
size_t i=0; i<numFieldsLo; i++)
593 if(!representative_node_candidates[i].size())
594 throw std::runtime_error(
"ERROR: GenerateRepresentativeBasisNodes: No candidates found!");
607template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
609 const std::vector<bool> & hi_nodeIsOwned,
611 const std::vector<size_t> &lo_node_in_hi,
612 const Basis &lo_basis,
613 const std::vector<LocalOrdinal> & hi_to_lo_map,
614 const Teuchos::RCP<const Map> & lo_colMap,
615 const Teuchos::RCP<const Map> & lo_domainMap,
616 const Teuchos::RCP<const Map> & hi_map,
617 Teuchos::RCP<Matrix>& P)
const{
620 size_t numFieldsHi = hi_elemToNode.extent(1);
621 size_t numFieldsLo = lo_basis.getCardinality();
622 LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
623 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
624 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
625 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
626 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
629 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
630 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
631 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
634 P = rcp(
new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi));
635 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
638 size_t Nelem=hi_elemToNode.extent(0);
639 std::vector<bool> touched(hi_map->getLocalNumElements(),
false);
640 Teuchos::Array<GO> col_gid(1);
641 Teuchos::Array<SC> val(1);
642 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
643 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
644 for(
size_t i=0; i<Nelem; i++) {
645 for(
size_t j=0; j<numFieldsHi; j++) {
646 LO row_lid = hi_elemToNode_host(i,j);
647 GO row_gid = hi_map->getGlobalElement(row_lid);
648 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
649 for(
size_t k=0; k<numFieldsLo; k++) {
651 LO col_lid = hi_to_lo_map[hi_elemToNode_host(i,lo_node_in_hi[k])];
652 if(col_lid==LOINVALID)
continue;
654 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
655 val[0] = LoValues_at_HiDofs_host(k,j);
658 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
659 P->insertGlobalValues(row_gid,col_gid(),val());
661 touched[row_lid]=
true;
665 P->fillComplete(lo_domainMap,hi_map);
669template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
671 const std::vector<bool> & hi_nodeIsOwned,
674 const Basis &lo_basis,
675 const std::vector<LocalOrdinal> & hi_to_lo_map,
676 const Teuchos::RCP<const Map> & lo_colMap,
677 const Teuchos::RCP<const Map> & lo_domainMap,
678 const Teuchos::RCP<const Map> & hi_map,
679 Teuchos::RCP<Matrix>& P)
const{
682 size_t numFieldsHi = hi_elemToNode.extent(1);
683 size_t numFieldsLo = lo_basis.getCardinality();
684 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
685 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
686 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
687 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
688 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
689 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
690 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
691 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
694 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
695 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
696 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
699 P = rcp(
new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi));
700 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
703 size_t Nelem=hi_elemToNode.extent(0);
704 std::vector<bool> touched(hi_map->getLocalNumElements(),
false);
705 Teuchos::Array<GO> col_gid(1);
706 Teuchos::Array<SC> val(1);
707 for(
size_t i=0; i<Nelem; i++) {
708 for(
size_t j=0; j<numFieldsHi; j++) {
709 LO row_lid = hi_elemToNode_host(i,j);
710 GO row_gid = hi_map->getGlobalElement(row_lid);
711 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
712 for(
size_t k=0; k<numFieldsLo; k++) {
714 LO col_lid = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,k)];
715 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
716 val[0] = LoValues_at_HiDofs_host(k,j);
719 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
720 P->insertGlobalValues(row_gid,col_gid(),val());
722 touched[row_lid]=
true;
726 P->fillComplete(lo_domainMap,hi_map);
730 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
732 RCP<ParameterList> validParamList = rcp(
new ParameterList());
734#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
737#undef SET_VALID_ENTRY
739 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
741 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
742 validParamList->set< RCP<const FactoryBase> >(
"pcoarsen: element to node map", Teuchos::null,
"Generating factory of the element to node map");
743 return validParamList;
747 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
749 Input(fineLevel,
"A");
750 Input(fineLevel,
"pcoarsen: element to node map");
751 Input(fineLevel,
"Nullspace");
755 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
757 return BuildP(fineLevel, coarseLevel);
761 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
765 const std::string prefix =
"MueLu::IntrepidPCoarsenFactory(" + levelIDs +
"): ";
768 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
769 typedef Kokkos::DynRankView<double,typename Node::device_type> FC;
772 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
773 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> >(fineLevel,
"Nullspace");
774 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Acrs =
dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*A);
777 if (restrictionMode_) {
783 std::vector<LocalOrdinal> A_dirichletRows;
790 RCP<ParameterList> APparams = rcp(
new ParameterList);
791 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
792 GetOStream(
static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
794 APparams = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data",
this);
796 if (APparams->isParameter(
"graph"))
797 finalP = APparams->get< RCP<Matrix> >(
"graph");
799 const ParameterList& pL = GetParameterList();
806 int lo_degree, hi_degree;
807 RCP<Basis> hi_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>(
"pcoarsen: hi basis"),hi_degree);
808 RCP<Basis> lo_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>(
"pcoarsen: lo basis"),lo_degree);
811 GetOStream(
Statistics1) <<
"P-Coarsening from basis "<<pL.get<std::string>(
"pcoarsen: hi basis")<<
" to "<<pL.get<std::string>(
"pcoarsen: lo basis") <<std::endl;
815 const Teuchos::RCP<FCi> Pn_elemToNode = Get<Teuchos::RCP<FCi> >(fineLevel,
"pcoarsen: element to node map");
822 RCP<const Map> rowMap = A->getRowMap();
823 RCP<const Map> colMap = Acrs.getColMap();
824 RCP<const Map> domainMap = A->getDomainMap();
825 int NumProc = rowMap->getComm()->getSize();
826 assert(rowMap->isSameAs(*domainMap));
827 std::vector<bool> Pn_nodeIsOwned(colMap->getLocalNumElements(),
false);
828 LO num_owned_rows = 0;
829 for(
size_t i=0; i<rowMap->getLocalNumElements(); i++) {
830 if(rowMap->getGlobalElement(i) == colMap->getGlobalElement(i)) {
831 Pn_nodeIsOwned[i] =
true;
838 Teuchos::RCP<FCi> P1_elemToNode = rcp(
new FCi());
840 std::vector<bool> P1_nodeIsOwned;
841 int P1_numOwnedNodes;
842 std::vector<LO> hi_to_lo_map;
845 std::vector<size_t> lo_node_in_hi;
848 FCi lo_elemToHiRepresentativeNode;
851 RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> > hi_isDirichletRow, hi_isDirichletCol;
855 printf(
"[%d] isDirichletRow = ",A->getRowMap()->getComm()->getRank());
856 for(
size_t i=0;i<hi_isDirichletRow->getMap()->getLocalNumElements(); i++)
857 printf(
"%d ",hi_isDirichletRow->getData(0)[i]);
859 printf(
"[%d] isDirichletCol = ",A->getRowMap()->getComm()->getRank());
860 for(
size_t i=0;i<hi_isDirichletCol->getMap()->getLocalNumElements(); i++)
861 printf(
"%d ",hi_isDirichletCol->getData(0)[i]);
872 MueLuIntrepid::BuildLoElemToNode(*Pn_elemToNode,Pn_nodeIsOwned,lo_node_in_hi,hi_isDirichletCol->getData(0),*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
873 assert(hi_to_lo_map.size() == colMap->getLocalNumElements());
877 double threshold = 1e-10;
878 std::vector<std::vector<size_t> > candidates;
879 Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
880 hi_basis->getDofCoords(hi_DofCoords);
882 MueLu::MueLuIntrepid::GenerateRepresentativeBasisNodes<Basis,FC>(*lo_basis,hi_DofCoords,threshold,candidates);
893 RCP<const Map> P1_domainMap = MapFactory::Build(rowMap->lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),P1_numOwnedNodes,rowMap->getIndexBase(),rowMap->getComm());
897 RCP<const Map> P1_colMap;
898 if(NumProc==1) P1_colMap = P1_domainMap;
899 else MueLuIntrepid::GenerateColMapFromImport<LO,GO,NO>(*Acrs.getCrsGraph()->getImporter(),hi_to_lo_map,*P1_domainMap,P1_nodeIsOwned.size(),P1_colMap);
904 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);
906 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);
914 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P1_domainMap, fineNullspace->getNumVectors());
915 finalP->apply(*fineNullspace,*coarseNullspace,Teuchos::TRANS);
916 Set(coarseLevel,
"Nullspace", coarseNullspace);
919 if (!restrictionMode_) {
921 Set(coarseLevel,
"P", finalP);
923 APparams->set(
"graph", finalP);
927 RCP<ParameterList> params = rcp(
new ParameterList());
928 params->set(
"printLoadBalancingInfo",
true);
929 params->set(
"printCommInfo",
true);
940 Set(coarseLevel,
"R", R);
943 RCP<ParameterList> params = rcp(
new ParameterList());
944 params->set(
"printLoadBalancingInfo",
true);
945 params->set(
"printCommInfo",
true);
#define SET_VALID_ENTRY(name)
#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level, ename, entry)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Kokkos::DynRankView< LocalOrdinal, typename Node::device_type > LOFieldContainer
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
Kokkos::DynRankView< double, typename Node::device_type > SCFieldContainer
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Intrepid2::Basis< typename Node::device_type::execution_space, double, double > Basis
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
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
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
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)
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)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
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)
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)
void GenerateRepresentativeBasisNodes(const Basis &basis, const SCFieldContainer &ReferenceNodeLocations, const double threshold, std::vector< std::vector< size_t > > &representative_node_candidates)
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)
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)
std::string tolower(const std::string &str)
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)
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)
Teuchos::RCP< Intrepid2::Basis< KokkosExecutionSpace, Scalar, Scalar > > BasisFactory(const std::string &name, int °ree)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.