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);} 104 namespace MueLuIntrepid {
105 inline 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());
196 template<
class Scalar,
class KokkosExecutionSpace>
197 Teuchos::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);
250 template<
class Scalar,
class KokkosDeviceType>
251 void 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);
285 template<
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);
331 template <
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);
408 template <
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);
513 ArrayRCP<GO> dvec_data = dvec->getDataNonConst(0);
514 for(
size_t i=0; i<hi_domainMap->getNodeNumElements(); i++) {
515 if(hi_to_lo_map[i]!=lo_invalid) dvec_data[i] = lo_domainMap.getGlobalElement(hi_to_lo_map[i]);
516 else dvec_data[i] = go_invalid;
520 RCP<GOVector> cvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_columnMap,
true);
521 cvec->doImport(*dvec,hi_importer,Xpetra::ADD);
525 Array<GO> lo_col_data(lo_columnMapLength);
526 ArrayRCP<GO> cvec_data = cvec->getDataNonConst(0);
527 for(
size_t i=0,idx=0; i<hi_columnMap->getNodeNumElements(); i++) {
528 if(hi_to_lo_map[i]!=lo_invalid) {
529 lo_col_data[idx] = cvec_data[i];
534 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());
545 template<
class Basis,
class SCFieldContainer>
546 void GenerateRepresentativeBasisNodes(
const Basis & basis,
const SCFieldContainer & ReferenceNodeLocations,
const double threshold,std::vector<std::vector<size_t> > & representative_node_candidates) {
547 typedef SCFieldContainer FC;
548 typedef typename FC::data_type SC;
551 size_t numFieldsHi = ReferenceNodeLocations.extent(0);
553 size_t numFieldsLo = basis.getCardinality();
555 FC LoValues(
"LoValues",numFieldsLo,numFieldsHi);
557 basis.getValues(LoValues, ReferenceNodeLocations , Intrepid2::OPERATOR_VALUE);
562 printf(
"** LoValues[%d,%d] **\n",(
int)numFieldsLo,(
int)numFieldsHi);
563 for(
size_t i=0; i<numFieldsLo; i++) {
564 for(
size_t j=0; j<numFieldsHi; j++)
565 printf(
"%6.4e ",LoValues(i,j));
568 printf(
"**************\n");fflush(stdout);
571 representative_node_candidates.resize(numFieldsLo);
572 auto LoValues_host = Kokkos::create_mirror_view(LoValues);
573 Kokkos::deep_copy(LoValues_host, LoValues);
574 for(
size_t i=0; i<numFieldsLo; i++) {
576 typename Teuchos::ScalarTraits<SC>::magnitudeType vmax = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero();
577 for(
size_t j=0; j<numFieldsHi; j++)
578 vmax = std::max(vmax,Teuchos::ScalarTraits<SC>::magnitude(LoValues_host(i,j)));
581 for(
size_t j=0; j<numFieldsHi; j++) {
582 if(Teuchos::ScalarTraits<SC>::magnitude(vmax - LoValues_host(i,j)) < threshold*vmax)
583 representative_node_candidates[i].push_back(j);
588 for(
size_t i=0; i<numFieldsLo; i++)
589 if(!representative_node_candidates[i].size())
590 throw std::runtime_error(
"ERROR: GenerateRepresentativeBasisNodes: No candidates found!");
603 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
605 const std::vector<bool> & hi_nodeIsOwned,
607 const std::vector<size_t> &lo_node_in_hi,
608 const Basis &lo_basis,
609 const std::vector<LocalOrdinal> & hi_to_lo_map,
610 const Teuchos::RCP<const Map> & lo_colMap,
611 const Teuchos::RCP<const Map> & lo_domainMap,
612 const Teuchos::RCP<const Map> & hi_map,
613 Teuchos::RCP<Matrix>& P)
const{
616 size_t numFieldsHi = hi_elemToNode.extent(1);
617 size_t numFieldsLo = lo_basis.getCardinality();
618 LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
619 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
620 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
621 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
622 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
625 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
626 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
627 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
630 P = rcp(
new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi));
631 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
634 size_t Nelem=hi_elemToNode.extent(0);
635 std::vector<bool> touched(hi_map->getNodeNumElements(),
false);
636 Teuchos::Array<GO> col_gid(1);
637 Teuchos::Array<SC> val(1);
638 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
639 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
640 for(
size_t i=0; i<Nelem; i++) {
641 for(
size_t j=0; j<numFieldsHi; j++) {
642 LO row_lid = hi_elemToNode_host(i,j);
643 GO row_gid = hi_map->getGlobalElement(row_lid);
644 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
645 for(
size_t k=0; k<numFieldsLo; k++) {
647 LO col_lid = hi_to_lo_map[hi_elemToNode_host(i,lo_node_in_hi[k])];
648 if(col_lid==LOINVALID)
continue;
650 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
651 val[0] = LoValues_at_HiDofs_host(k,j);
654 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
655 P->insertGlobalValues(row_gid,col_gid(),val());
657 touched[row_lid]=
true;
661 P->fillComplete(lo_domainMap,hi_map);
665 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
667 const std::vector<bool> & hi_nodeIsOwned,
670 const Basis &lo_basis,
671 const std::vector<LocalOrdinal> & hi_to_lo_map,
672 const Teuchos::RCP<const Map> & lo_colMap,
673 const Teuchos::RCP<const Map> & lo_domainMap,
674 const Teuchos::RCP<const Map> & hi_map,
675 Teuchos::RCP<Matrix>& P)
const{
678 size_t numFieldsHi = hi_elemToNode.extent(1);
679 size_t numFieldsLo = lo_basis.getCardinality();
680 FC LoValues_at_HiDofs(
"LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
681 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
682 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
683 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
684 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
685 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
686 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
687 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
690 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
691 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
692 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
695 P = rcp(
new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi));
696 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
699 size_t Nelem=hi_elemToNode.extent(0);
700 std::vector<bool> touched(hi_map->getNodeNumElements(),
false);
701 Teuchos::Array<GO> col_gid(1);
702 Teuchos::Array<SC> val(1);
703 for(
size_t i=0; i<Nelem; i++) {
704 for(
size_t j=0; j<numFieldsHi; j++) {
705 LO row_lid = hi_elemToNode_host(i,j);
706 GO row_gid = hi_map->getGlobalElement(row_lid);
707 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
708 for(
size_t k=0; k<numFieldsLo; k++) {
710 LO col_lid = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,k)];
711 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
712 val[0] = LoValues_at_HiDofs_host(k,j);
715 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
716 P->insertGlobalValues(row_gid,col_gid(),val());
718 touched[row_lid]=
true;
722 P->fillComplete(lo_domainMap,hi_map);
726 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
728 RCP<ParameterList> validParamList = rcp(
new ParameterList());
730 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 733 #undef SET_VALID_ENTRY 735 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
737 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
738 validParamList->set< RCP<const FactoryBase> >(
"pcoarsen: element to node map", Teuchos::null,
"Generating factory of the element to node map");
739 return validParamList;
743 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
745 Input(fineLevel,
"A");
746 Input(fineLevel,
"pcoarsen: element to node map");
747 Input(fineLevel,
"Nullspace");
751 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
753 return BuildP(fineLevel, coarseLevel);
757 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
761 const std::string prefix =
"MueLu::IntrepidPCoarsenFactory(" + levelIDs +
"): ";
764 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
765 typedef Kokkos::DynRankView<double,typename Node::device_type> FC;
768 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
769 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> >(fineLevel,
"Nullspace");
770 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Acrs =
dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
>(*A);
773 if (restrictionMode_) {
779 std::vector<LocalOrdinal> A_dirichletRows;
786 RCP<ParameterList> APparams = rcp(
new ParameterList);
787 if (coarseLevel.
IsAvailable(
"AP reuse data",
this)) {
788 GetOStream(static_cast<MsgType>(
Runtime0 |
Test)) <<
"Reusing previous AP data" << std::endl;
790 APparams = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data",
this);
792 if (APparams->isParameter(
"graph"))
793 finalP = APparams->get< RCP<Matrix> >(
"graph");
795 const ParameterList& pL = GetParameterList();
802 int lo_degree, hi_degree;
803 RCP<Basis> hi_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>(
"pcoarsen: hi basis"),hi_degree);
804 RCP<Basis> lo_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>(
"pcoarsen: lo basis"),lo_degree);
807 GetOStream(
Statistics1) <<
"P-Coarsening from basis "<<pL.get<std::string>(
"pcoarsen: hi basis")<<
" to "<<pL.get<std::string>(
"pcoarsen: lo basis") <<std::endl;
811 const Teuchos::RCP<FCi> Pn_elemToNode = Get<Teuchos::RCP<FCi> >(fineLevel,
"pcoarsen: element to node map");
818 RCP<const Map> rowMap = A->getRowMap();
819 RCP<const Map> colMap = Acrs.getColMap();
820 RCP<const Map> domainMap = A->getDomainMap();
821 int NumProc = rowMap->getComm()->getSize();
822 assert(rowMap->isSameAs(*domainMap));
823 std::vector<bool> Pn_nodeIsOwned(colMap->getNodeNumElements(),
false);
824 LO num_owned_rows = 0;
825 for(
size_t i=0; i<rowMap->getNodeNumElements(); i++) {
826 if(rowMap->getGlobalElement(i) == colMap->getGlobalElement(i)) {
827 Pn_nodeIsOwned[i] =
true;
834 Teuchos::RCP<FCi> P1_elemToNode = rcp(
new FCi());
836 std::vector<bool> P1_nodeIsOwned;
837 int P1_numOwnedNodes;
838 std::vector<LO> hi_to_lo_map;
841 std::vector<size_t> lo_node_in_hi;
844 FCi lo_elemToHiRepresentativeNode;
847 RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> > hi_isDirichletRow, hi_isDirichletCol;
851 printf(
"[%d] isDirichletRow = ",A->getRowMap()->getComm()->getRank());
852 for(
size_t i=0;i<hi_isDirichletRow->getMap()->getNodeNumElements(); i++)
853 printf(
"%d ",hi_isDirichletRow->getData(0)[i]);
855 printf(
"[%d] isDirichletCol = ",A->getRowMap()->getComm()->getRank());
856 for(
size_t i=0;i<hi_isDirichletCol->getMap()->getNodeNumElements(); i++)
857 printf(
"%d ",hi_isDirichletCol->getData(0)[i]);
868 MueLuIntrepid::BuildLoElemToNode(*Pn_elemToNode,Pn_nodeIsOwned,lo_node_in_hi,hi_isDirichletCol->getData(0),*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
869 assert(hi_to_lo_map.size() == colMap->getNodeNumElements());
873 double threshold = 1e-10;
874 std::vector<std::vector<size_t> > candidates;
875 Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
876 hi_basis->getDofCoords(hi_DofCoords);
878 MueLu::MueLuIntrepid::GenerateRepresentativeBasisNodes<Basis,FC>(*lo_basis,hi_DofCoords,threshold,candidates);
889 RCP<const Map> P1_domainMap = MapFactory::Build(rowMap->lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),P1_numOwnedNodes,rowMap->getIndexBase(),rowMap->getComm());
893 RCP<const Map> P1_colMap;
894 if(NumProc==1) P1_colMap = P1_domainMap;
895 else MueLuIntrepid::GenerateColMapFromImport<LO,GO,NO>(*Acrs.getCrsGraph()->getImporter(),hi_to_lo_map,*P1_domainMap,P1_nodeIsOwned.size(),P1_colMap);
900 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);
902 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);
910 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P1_domainMap, fineNullspace->getNumVectors());
911 finalP->apply(*fineNullspace,*coarseNullspace,Teuchos::TRANS);
912 Set(coarseLevel,
"Nullspace", coarseNullspace);
915 if (!restrictionMode_) {
917 Set(coarseLevel,
"P", finalP);
919 APparams->set(
"graph", finalP);
923 RCP<ParameterList> params = rcp(
new ParameterList());
924 params->set(
"printLoadBalancingInfo",
true);
925 params->set(
"printCommInfo",
true);
936 Set(coarseLevel,
"R", R);
939 RCP<ParameterList> params = rcp(
new ParameterList());
940 params->set(
"printLoadBalancingInfo",
true);
941 params->set(
"printCommInfo",
true);
950 #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)
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)
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)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
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
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
std::string tolower(const std::string &str)
Namespace for MueLu classes and methods.
Print even more statistics.
int GetLevelID() const
Return level number.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Kokkos::DynRankView< double, typename Node::device_type > SCFieldContainer
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)
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)
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 void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
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)
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
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Kokkos::DynRankView< LocalOrdinal, typename Node::device_type > LOFieldContainer
Intrepid2::Basis< typename Node::device_type::execution_space, double, double > Basis
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 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 GenerateRepresentativeBasisNodes(const Basis &basis, const SCFieldContainer &ReferenceNodeLocations, const double threshold, std::vector< std::vector< size_t > > &representative_node_candidates)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
#define SET_VALID_ENTRY(name)