43 #ifndef PANZER_DOF_MANAGER2_IMPL_HPP 44 #define PANZER_DOF_MANAGER2_IMPL_HPP 51 #include "PanzerDofMgr_config.hpp" 59 #include "Teuchos_GlobalMPISession.hpp" 62 #include "Teuchos_RCP.hpp" 63 #include "Teuchos_Array.hpp" 64 #include "Teuchos_ArrayView.hpp" 66 #include "Tpetra_Map.hpp" 67 #include "Tpetra_Export.hpp" 68 #include "Tpetra_Vector.hpp" 69 #include "Tpetra_MultiVector.hpp" 71 #include <unordered_set> 73 #define PANZER_DOFMGR_FUNC_TIME_MONITOR(a) \ 74 PANZER_FUNC_TIME_MONITOR(a) 76 #ifdef PHX_KOKKOS_DEVICE_TYPE_CUDA 77 #define PANZER_DOFMGR_REQUIRE_CUDA 91 template <
typename LocalOrdinal,
typename GlobalOrdinal>
92 class HashTieBreak :
public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
96 HashTieBreak(
const unsigned int seed = (2654435761U))
99 virtual std::size_t selectedIndex(GlobalOrdinal GID,
100 const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid)
const 103 int intkey = (int) ((GID & 0x000000007fffffffLL) +
104 ((GID & 0x7fffffff80000000LL) >> 31));
105 return std::size_t((
seed_ ^ intkey) % pid_and_lid.size());
113 template <
typename LocalOrdinal,
typename GlobalOrdinal>
114 class GreedyTieBreak :
public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
119 virtual bool mayHaveSideEffects()
const {
123 virtual std::size_t selectedIndex(GlobalOrdinal ,
124 const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid)
const 127 auto numLids = pid_and_lid.size();
128 decltype(numLids) idx = 0;
129 auto minpid = pid_and_lid[0].first;
130 decltype(minpid) minidx = 0;
131 for (idx = 0; idx < numLids; ++idx) {
132 if (pid_and_lid[idx].first < minpid) {
133 minpid = pid_and_lid[idx].first;
146 using Teuchos::ArrayRCP;
147 using Teuchos::Array;
148 using Teuchos::ArrayView;
150 template <
typename LO,
typename GO>
152 : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
155 template <
typename LO,
typename GO>
157 : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
162 template <
typename LO,
typename GO>
165 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
166 "DOFManager::setConnManager: setConnManager cannot be called after " 167 "buildGlobalUnknowns has been called");
168 connMngr_ = connMngr;
170 connMngr_->getElementBlockIds(blockOrder_);
171 for (
size_t i = 0; i < blockOrder_.size(); ++i) {
172 blockNameToID_.insert(std::map<std::string,int>::value_type(blockOrder_[i],i));
175 blockToAssociatedFP_.resize(blockOrder_.size());
182 template <
typename LO,
typename GO>
186 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
187 "DOFManager::addField: addField cannot be called after " 188 "buildGlobalUnknowns has been called");
190 fieldPatterns_.push_back(pattern);
191 fieldTypes_.push_back(type);
192 fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
195 fieldStringOrder_.push_back(str);
196 fieldAIDOrder_.push_back(numFields_);
198 for(std::size_t i=0;i<blockOrder_.size();i++) {
199 blockToAssociatedFP_[i].push_back(numFields_);
206 template <
typename LO,
typename GO>
210 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
211 "DOFManager::addField: addField cannot be called after " 212 "buildGlobalUnknowns has been called");
213 TEUCHOS_TEST_FOR_EXCEPTION((connMngr_==Teuchos::null),std::logic_error,
214 "DOFManager::addField: you must add a ConnManager before" 215 "you can associate a FP with a given block.")
218 while(!found && blocknum<blockOrder_.size()){
219 if(blockOrder_[blocknum]==blockID){
225 TEUCHOS_TEST_FOR_EXCEPTION(!found,std::logic_error,
"DOFManager::addField: Invalid block name.");
230 std::map<std::string,int>::const_iterator fpIter = fieldNameToAID_.find(str);
231 if(fpIter!=fieldNameToAID_.end())
235 fieldPatterns_.push_back(pattern);
236 fieldTypes_.push_back(type);
237 fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
239 fieldStringOrder_.push_back(str);
240 fieldAIDOrder_.push_back(numFields_);
243 blockToAssociatedFP_[blocknum].push_back(numFields_);
248 blockToAssociatedFP_[blocknum].push_back(fpIter->second);
256 template <
typename LO,
typename GO>
259 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(name);
260 if(fitr==fieldNameToAID_.end())
261 return Teuchos::null;
263 if(fitr->second<
int(fieldPatterns_.size()))
264 return fieldPatterns_[fitr->second];
266 return Teuchos::null;
270 template <
typename LO,
typename GO>
273 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(fieldName);
274 if(fitr==fieldNameToAID_.end())
275 return Teuchos::null;
279 while(!found && blocknum<blockOrder_.size()){
280 if(blockOrder_[blocknum]==blockId){
288 return Teuchos::null;
290 std::vector<int>::const_iterator itr = std::find(blockToAssociatedFP_[blocknum].begin(),
291 blockToAssociatedFP_[blocknum].end(),fitr->second);
292 if(itr!=blockToAssociatedFP_[blocknum].end()) {
293 if(fitr->second<
int(fieldPatterns_.size()))
294 return fieldPatterns_[fitr->second];
297 return Teuchos::null;
305 template<
typename LO,
typename GO>
309 std::vector<GO>& indices)
const 319 template<
typename LO,
typename GO>
323 std::vector<GO>& indices)
const 333 template<
typename LO,
typename GO>
337 std::vector<GO>& indices)
const 340 indices.resize(owned_.size() + ghosted_.size());
341 for (
size_t i(0); i < owned_.size(); ++i)
342 indices[i] = owned_[i];
343 for (
size_t i(0); i < ghosted_.size(); ++i)
344 indices[owned_.size() + i] = ghosted_[i];
352 template<
typename LO,
typename GO>
357 return owned_.size();
365 template<
typename LO,
typename GO>
370 return ghosted_.size();
378 template<
typename LO,
typename GO>
383 return owned_.size() + ghosted_.size();
387 template <
typename LO,
typename GO>
393 template <
typename LO,
typename GO>
396 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,
"DOFManager::getGIDFieldOffsets: cannot be called before " 397 "buildGlobalUnknowns has been called");
398 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
399 if(bitr==blockNameToID_.end())
400 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid block name");
401 int bid=bitr->second;
402 if(fa_fps_[bid]!=Teuchos::null)
403 return fa_fps_[bid]->localOffsets(fieldNum);
405 static const std::vector<int> empty;
409 template <
typename LO,
typename GO>
413 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,
"DOFManager::getGIDFieldOffsets: cannot be called before " 414 "buildGlobalUnknowns has been called");
415 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
416 if(bitr==blockNameToID_.end())
417 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid block name");
418 int bid=bitr->second;
419 if(fa_fps_[bid]!=Teuchos::null)
420 return fa_fps_[bid]->localOffsetsKokkos(fieldNum);
422 static const Kokkos::View<int*,PHX::Device> empty(
"panzer::DOFManager::getGIDFieldOffsetsKokkos() empty",0);
426 template <
typename LO,
typename GO>
429 gids = elementGIDs_[localElementID];
432 template <
typename LocalOrdinalT,
typename GlobalOrdinalT>
438 if(requireOrientations_){
439 fieldPatterns_.push_back(Teuchos::rcp(
new NodalFieldPattern(fieldPatterns_[0]->getCellTopology())));
443 TEUCHOS_ASSERT(fieldPatterns_.size() == fieldTypes_.size());
444 std::vector<std::pair<FieldType,Teuchos::RCP<const FieldPattern>>> tmp;
445 for (std::size_t i=0; i < fieldPatterns_.size(); ++i)
446 tmp.push_back(std::make_pair(fieldTypes_[i],fieldPatterns_[i]));
450 connMngr_->buildConnectivity(*aggFieldPattern);
453 buildGlobalUnknowns(aggFieldPattern);
457 template <
typename LO,
typename GO>
462 typedef Tpetra::Map<LO, GO, Node> Map;
464 typedef Tpetra::Import<LO,GO,Node> Import;
467 typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
471 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
472 out.setOutputToRootOnly(-1);
473 out.setShowProcRank(
true);
475 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
476 "DOFManager::buildGlobalUnknowns: buildGlobalUnknowns cannot be called again " 477 "after buildGlobalUnknowns has been called");
481 if(getOrientationsRequired()) {
482 std::size_t sz = geomPattern->getSubcellIndices(0,0).size();
484 TEUCHOS_TEST_FOR_EXCEPTION(sz==0,std::logic_error,
485 "DOFManager::buildGlobalUnknowns requires a geometric pattern including " 486 "the nodes when orientations are needed!");
492 ga_fp_ = geomPattern;
500 RCP<MultiVector> tagged_overlap_mv = buildTaggedMultiVector(ownedAccess);
501 RCP<const Map> overlap_map = tagged_overlap_mv->getMap();
503 RCP<MultiVector> overlap_mv = Tpetra::createMultiVector<GO>(overlap_map,(size_t)numFields_);
506 auto non_overlap_pair = buildGlobalUnknowns_GUN(*tagged_overlap_mv,*overlap_mv);
507 RCP<MultiVector> non_overlap_mv = non_overlap_pair.first;
508 RCP<MultiVector> tagged_non_overlap_mv = non_overlap_pair.second;
509 RCP<const Map> non_overlap_map = non_overlap_mv->getMap();
517 fillGIDsFromOverlappedMV(ownedAccess,elementGIDs_,*overlap_map,*overlap_mv);
524 RCP<const Map> overlap_map_neighbor =
525 buildOverlapMapFromElements(neighborAccess);
528 Import imp_neighbor(non_overlap_map,overlap_map_neighbor);
530 Teuchos::RCP<MultiVector> overlap_mv_neighbor =
531 Tpetra::createMultiVector<GO>(overlap_map_neighbor, (size_t)numFields_);
534 overlap_mv_neighbor->doImport(*non_overlap_mv, imp_neighbor,
537 fillGIDsFromOverlappedMV(neighborAccess, elementGIDs_,
538 *overlap_map_neighbor, *overlap_mv_neighbor);
547 panzer::Ordinal64 offset = 0xFFFFFFFFLL;
549 for (
size_t b = 0; b < blockOrder_.size(); ++b) {
550 const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
551 for (std::size_t l = 0; l < myElements.size(); ++l) {
552 std::vector<GO> & localGIDs = elementGIDs_[myElements[l]];
553 for(std::size_t c=0;c<localGIDs.size();c++)
554 localGIDs[c] += offset;
558 Teuchos::ArrayRCP<GO> nvals = non_overlap_mv->get1dViewNonConst();
559 for (
int j = 0; j < nvals.size(); ++j)
568 typedef std::unordered_set<GO> HashTable;
569 HashTable isOwned, remainingOwned;
572 Teuchos::ArrayRCP<const GO> nvals = non_overlap_mv->get1dView();
573 Teuchos::ArrayRCP<const GO> tagged_vals = tagged_non_overlap_mv->get1dView();
574 TEUCHOS_ASSERT(nvals.size()==tagged_vals.size());
575 for (
int j = 0; j < nvals.size(); ++j) {
576 if (nvals[j] != -1) {
577 for(
GO offset=0;offset<tagged_vals[j];offset++)
578 isOwned.insert(nvals[j]+offset);
582 TEUCHOS_ASSERT(tagged_vals[j]==0)
585 remainingOwned = isOwned;
588 for (
size_t b = 0; b < blockOrder_.size(); ++b) {
590 if(fa_fps_[b]==Teuchos::null)
593 const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
595 for (
size_t l = 0; l < myElements.size(); ++l) {
596 const std::vector<GO> & localOrdering = elementGIDs_[myElements[l]];
599 for(std::size_t i=0;i<localOrdering.size();i++) {
601 if(isOwned.find(localOrdering[i])==isOwned.end())
605 std::pair<typename HashTable::iterator,bool> insertResult = hashTable.insert(localOrdering[i]);
606 if(insertResult.second) {
607 owned_.push_back(localOrdering[i]);
608 remainingOwned.erase(localOrdering[i]);
617 for(
typename HashTable::const_iterator itr=remainingOwned.begin();itr!=remainingOwned.end();itr++)
618 owned_.push_back(*itr);
620 if(owned_.size()!=isOwned.size()) {
621 out <<
"I'm about to hang because of unknown numbering failure ... sorry! (line = " << __LINE__ <<
")" << std::endl;
622 TEUCHOS_TEST_FOR_EXCEPTION(owned_.size()!=isOwned.size(),std::logic_error,
623 "DOFManager::buildGlobalUnkonwns: Failure because not all owned unknowns have been accounted for.");
635 "buildGlobalUnknowns::build_ghosted_array");
638 typedef std::unordered_set<GO> HashTable;
640 for (std::size_t i = 0; i < owned_.size(); i++)
641 hashTable.insert(owned_[i]);
646 std::vector<ElementBlockAccess> blockAccessVec;
650 for (std::size_t a = 0; a < blockAccessVec.size(); ++a)
654 for (
size_t b = 0; b < blockOrder_.size(); ++b)
656 if (fa_fps_[b] == Teuchos::null)
658 const std::vector<LO>& myElements =
660 for (
size_t l = 0; l < myElements.size(); ++l)
662 const std::vector<GO>& localOrdering = elementGIDs_[myElements[l]];
665 for (std::size_t i = 0; i < localOrdering.size(); ++i)
667 std::pair<typename HashTable::iterator, bool> insertResult =
668 hashTable.insert(localOrdering[i]);
672 if(insertResult.second)
673 ghosted_.push_back(localOrdering[i]);
680 buildConnectivityRun_ =
true;
683 if(requireOrientations_) {
685 buildUnknownsOrientation();
691 this->buildLocalIdsFromOwnedAndGhostedElements();
695 this->buildLocalIds();
699 template <
typename LO,
typename GO>
700 std::pair<Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> >,
701 Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> > >
703 Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & overlap_mv)
const 707 typedef Tpetra::Map<LO, GO, Node> Map;
709 typedef Tpetra::Export<LO,GO,Node> Export;
710 typedef Tpetra::Import<LO,GO,Node> Import;
713 typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
718 RCP<const Map> overlap_map = tagged_overlap_mv.getMap();
725 RCP<const Map> non_overlap_map;
729 GreedyTieBreak<LO,GO> greedy_tie_break;
730 non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map, greedy_tie_break);
735 HashTieBreak<LO,GO> tie_break;
736 non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map,tie_break);
744 Teuchos::RCP<MultiVector> tagged_non_overlap_mv;
748 tagged_non_overlap_mv = Tpetra::createMultiVector<GO>(non_overlap_map,(size_t)numFields_);
757 RCP<MultiVector> non_overlap_mv;
761 exp = rcp(
new Export(overlap_map,non_overlap_map));
763 #ifdef PANZER_DOFMGR_REQUIRE_CUDA 767 imp = rcp(
new Import(non_overlap_map,overlap_map));
772 tagged_non_overlap_mv->doExport(tagged_overlap_mv,*exp,Tpetra::ABSMAX);
776 non_overlap_mv = rcp(
new MultiVector(*tagged_non_overlap_mv,Teuchos::Copy));
789 typedef typename Tpetra::MultiVector<GO,Node> MV;
790 typedef typename MV::dual_view_type::t_dev KV;
791 typedef typename MV::dual_view_type::t_dev::memory_space DMS;
792 KV values = non_overlap_mv->template getLocalView<DMS>();
793 auto mv_size = values.extent(0);
808 Teuchos::scan<int, GO> (*getComm(), Teuchos::REDUCE_SUM,
static_cast<size_t> (localsum), Teuchos::outArg (scanResult));
809 myOffset = scanResult - localsum;
826 ArrayRCP<ArrayRCP<GO> > editnonoverlap = non_overlap_mv->get2dViewNonConst();
827 for(
size_t i=0; i<non_overlap_mv->getLocalLength(); ++i){
828 for(
int j=0; j<numFields_; ++j){
829 if(editnonoverlap[j][i]!=0){
831 int ndof = Teuchos::as<int>(editnonoverlap[j][i]);
832 editnonoverlap[j][i]=myOffset+which_id;
836 editnonoverlap[j][i]=-1;
854 #ifdef PANZER_DOFMGR_REQUIRE_CUDA 855 overlap_mv.doImport(*non_overlap_mv,*imp,Tpetra::REPLACE);
858 overlap_mv.doImport(*non_overlap_mv,*exp,Tpetra::REPLACE);
865 return std::make_pair(non_overlap_mv,tagged_non_overlap_mv);
868 template <
typename LO,
typename GO>
869 Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> >
874 typedef Tpetra::Map<LO, GO, Node> Map;
875 typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
882 for (
size_t b = 0; b < blockOrder_.size(); ++b) {
883 std::vector<std::tuple< int, panzer::FieldType, RCP<const panzer::FieldPattern> > > faConstruct;
887 for (
size_t i = 0; i < fieldAIDOrder_.size(); ++i) {
888 int looking = fieldAIDOrder_[i];
891 std::vector<int>::const_iterator reu = std::find(blockToAssociatedFP_[b].begin(), blockToAssociatedFP_[b].end(), looking);
892 if(!(reu==blockToAssociatedFP_[b].end())){
893 faConstruct.push_back(std::make_tuple(i, fieldTypes_[fieldAIDOrder_[i]], fieldPatterns_[fieldAIDOrder_[i]]));
898 if(faConstruct.size()>0) {
902 int gidsInBlock = fa_fps_[fa_fps_.size()-1]->numberIds();
903 elementBlockGIDCount_.push_back(gidsInBlock);
906 fa_fps_.push_back(Teuchos::null);
907 elementBlockGIDCount_.push_back(0);
911 RCP<const Map> overlapmap = buildOverlapMapFromElements(ownedAccess);
916 Teuchos::RCP<MultiVector> overlap_mv;
920 overlap_mv = Tpetra::createMultiVector<GO>(overlapmap,(size_t)numFields_);
921 overlap_mv->putScalar(0);
932 std::vector<int> working(overlap_mv->getNumVectors());
933 ArrayRCP<ArrayRCP<GO> > edittwoview = overlap_mv->get2dViewNonConst();
934 for (
size_t b = 0; b < blockOrder_.size(); ++b) {
936 if(fa_fps_[b]==Teuchos::null)
939 const std::vector<LO> &
numFields= fa_fps_[b]->numFieldsPerId();
940 const std::vector<LO> & fieldIds= fa_fps_[b]->fieldIds();
941 const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
942 for (
size_t l = 0; l < myElements.size(); ++l) {
943 LO connSize = connMngr_->getConnectivitySize(myElements[l]);
944 const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
946 for (
int c = 0; c < connSize; ++c) {
947 size_t lid = overlapmap->getLocalElement(elmtConn[c]);
949 for(std::size_t i=0;i<working.size();i++)
952 int whichField = fieldIds[offset];
955 working[whichField]++;
958 for(std::size_t i=0;i<working.size();i++) {
959 auto current = edittwoview[i][lid];
960 edittwoview[i][lid] = (current > working[i]) ? current : working[i];
978 template <
typename LO,
typename GO>
983 while(!found && (
size_t)ind<fieldStringOrder_.size()){
984 if(fieldStringOrder_[ind]==
string)
997 template <
typename LO,
typename GO>
1000 fieldOrder.resize(fieldStringOrder_.size());
1001 for (
size_t i = 0; i < fieldStringOrder_.size(); ++i)
1002 fieldOrder[i]=fieldStringOrder_[i];
1005 template <
typename LO,
typename GO>
1008 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(
field);
1009 if(fitr==fieldNameToAID_.end()) {
1010 std::stringstream ss;
1011 printFieldInformation(ss);
1012 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid field name. DOF information is:\n"+ss.str());
1014 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(block);
1015 if(bitr==blockNameToID_.end()) {
1016 std::stringstream ss;
1017 printFieldInformation(ss);
1018 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid block name. DOF information is:\n"+ss.str());
1020 int fid=fitr->second;
1021 int bid=bitr->second;
1024 for (
size_t i = 0; i < blockToAssociatedFP_[bid].size(); ++i) {
1025 if(blockToAssociatedFP_[bid][i]==fid){
1034 template <
typename LO,
typename GO>
1037 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,
"DOFManager::getBlockFieldNumbers: BuildConnectivity must be run first.");
1039 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1040 if(bitr==blockNameToID_.end())
1041 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid block name");
1042 int bid=bitr->second;
1045 if(fa_fps_[bid]!=Teuchos::null)
1046 return fa_fps_[bid]->fieldIds();
1049 static std::vector<int> empty;
1053 template <
typename LO,
typename GO>
1054 const std::pair<std::vector<int>,std::vector<int> > &
1057 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,
"DOFManager::getGIDFieldOffsets_closure: BuildConnectivity must be run first.");
1058 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1059 if(bitr==blockNameToID_.end())
1060 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::getGIDFieldOffsets_closure: invalid block name.");
1063 if(fa_fps_[bitr->second]!=Teuchos::null)
1064 return fa_fps_[bitr->second]->localOffsets_closure(fieldNum, subcellDim, subcellId);
1066 static std::pair<std::vector<int>,std::vector<int> > empty;
1070 template <
typename LO,
typename GO>
1074 if(indices.size()!=isOwned.size())
1075 isOwned.resize(indices.size(),
false);
1076 typename std::vector<GO>::const_iterator endOf = owned_.end();
1077 for (std::size_t i = 0; i < indices.size(); ++i) {
1078 isOwned[i] = ( std::find(owned_.begin(), owned_.end(), indices[i])!=endOf );
1082 template <
typename LO,
typename GO>
1085 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
1086 "DOFManager::setFieldOrder: setFieldOrder cannot be called after " 1087 "buildGlobalUnknowns has been called");
1088 if(validFieldOrder(fieldOrder)){
1089 fieldStringOrder_=fieldOrder;
1091 for (
size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1092 fieldAIDOrder_[i]=fieldNameToAID_[fieldStringOrder_[i]];
1096 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::setFieldOrder: Invalid Field Ordering!");
1101 template <
typename LO,
typename GO>
1104 if(fieldStringOrder_.size()!=proposed_fieldOrder.size())
1108 for (
size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1110 for (
size_t j = 0; j < proposed_fieldOrder.size(); ++j) {
1111 if(fieldStringOrder_[i]==proposed_fieldOrder[j])
1120 template <
typename LO,
typename GO>
1124 if(num>=(
int)fieldStringOrder_.size())
1125 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"DOFManager::getFieldString: invalid number");
1126 return fieldStringOrder_[num];
1134 template <
typename LO,
typename GO>
1137 orientation_.clear();
1139 std::vector<std::string> elementBlockIds;
1140 connMngr_->getElementBlockIds(elementBlockIds);
1143 std::size_t myElementCount = 0;
1144 for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin(); blockItr!=elementBlockIds.end();++blockItr)
1145 myElementCount += connMngr_->getElementBlock(*blockItr).size();
1148 orientation_.resize(myElementCount);
1151 for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin();
1152 blockItr!=elementBlockIds.end();++blockItr) {
1153 const std::string & blockName = *blockItr;
1156 std::map<std::string,int>::const_iterator fap = blockNameToID_.find(blockName);
1157 if(fap==blockNameToID_.end()) {
1158 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::buildUnknownsOrientation: invalid block name");
1161 int bid=fap->second;
1163 if(fa_fps_[bid]==Teuchos::null)
1170 std::vector<std::pair<int,int> > topEdgeIndices;
1174 std::vector<std::vector<int> > topFaceIndices;
1175 if(ga_fp_->getDimension()==3)
1179 const std::vector<LO> & elmts = connMngr_->getElementBlock(blockName);
1180 for(std::size_t e=0;e<elmts.size();e++) {
1182 std::vector<signed char> & eOrientation = orientation_[elmts[e]];
1187 eOrientation.resize(fieldPattern.
numberIds());
1188 for(std::size_t s=0;s<eOrientation.size();s++)
1189 eOrientation[s] = 1;
1192 LO connSz = connMngr_->getConnectivitySize(elmts[e]);
1193 const GO * connPtr = connMngr_->getConnectivity(elmts[e]);
1194 const std::vector<GO> connectivity(connPtr,connPtr+connSz);
1199 if(ga_fp_->getDimension()==3)
1205 template <
typename LO,
typename GO>
1208 TEUCHOS_TEST_FOR_EXCEPTION(orientation_.size()==0,std::logic_error,
1209 "DOFManager::getElementOrientations: Orientations were not constructed!");
1211 const std::vector<signed char> & local_o = orientation_[localElmtId];
1212 gidsOrientation.resize(local_o.size());
1213 for(std::size_t i=0;i<local_o.size();i++) {
1214 gidsOrientation[i] = double(local_o[i]);
1218 template <
typename LocalOrdinalT,
typename GlobalOrdinalT>
1221 Teuchos::RCP<ConnManager<LocalOrdinalT,GlobalOrdinalT> > connMngr = connMngr_;
1223 connMngr_ = Teuchos::null;
1226 orientation_.clear();
1228 elementGIDs_.clear();
1231 elementBlockGIDCount_.clear();
1236 template <
typename LocalOrdinalT,
typename GlobalOrdinalT>
1239 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1240 if(bitr==blockNameToID_.end())
1241 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"DOFManager::fieldInBlock: invalid block name");
1242 return bitr->second;
1245 template <
typename LocalOrdinalT,
typename GlobalOrdinalT>
1248 os <<
"DOFManager Field Information: " << std::endl;
1251 TEUCHOS_ASSERT(blockOrder_.size()==blockToAssociatedFP_.size());
1253 for(std::size_t i=0;i<blockOrder_.size();i++) {
1254 os <<
" Element Block = " << blockOrder_[i] << std::endl;
1257 const std::vector<int> & fieldIds = blockToAssociatedFP_[i];
1258 for(std::size_t f=0;f<fieldIds.size();f++)
1259 os <<
" \"" << getFieldString(fieldIds[f]) <<
"\" is field ID " << fieldIds[f] << std::endl;
1263 template <
typename LO,
typename GO>
1264 Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> >
1275 std::set<GO> overlapset;
1276 for (
size_t i = 0; i < blockOrder_.size(); ++i) {
1277 const std::vector<LO> & myElements = access.
getElementBlock(blockOrder_[i]);
1278 for (
size_t e = 0; e < myElements.size(); ++e) {
1279 LO connSize = connMngr_->getConnectivitySize(myElements[e]);
1280 const GO * elmtConn = connMngr_->getConnectivity(myElements[e]);
1281 for (
int k = 0; k < connSize; ++k) {
1282 overlapset.insert(elmtConn[k]);
1287 Array<GO> overlapVector;
1288 for (
typename std::set<GO>::const_iterator itr = overlapset.begin(); itr!=overlapset.end(); ++itr) {
1289 overlapVector.push_back(*itr);
1294 return Tpetra::createNonContigMap<LO,GO>(overlapVector,getComm());
1297 template <
typename LO,
typename GO>
1300 std::vector<std::vector< GO > > & elementGIDs,
1301 const Tpetra::Map<LO,GO,panzer::TpetraNodeType> & overlapmap,
1302 const Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & overlap_mv)
const 1304 using Teuchos::ArrayRCP;
1307 ArrayRCP<ArrayRCP<const GO> > twoview = overlap_mv.get2dView();
1311 for (
size_t b = 0; b < blockOrder_.size(); ++b) {
1312 const std::vector<LO> & myElements = access.
getElementBlock(blockOrder_[b]);
1314 if(fa_fps_[b]==Teuchos::null) {
1316 for (
size_t l = 0; l < myElements.size(); ++l) {
1317 LO thisID=myElements[l];
1318 if(elementGIDs.size()<=(size_t)thisID)
1319 elementGIDs.resize(thisID+1);
1324 const std::vector<int> &
numFields= fa_fps_[b]->numFieldsPerId();
1325 const std::vector<int> & fieldIds= fa_fps_[b]->fieldIds();
1328 for (
size_t l = 0; l < myElements.size(); ++l) {
1329 LO connSize = connMngr_->getConnectivitySize(myElements[l]);
1330 const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
1331 std::vector<GO> localOrdering;
1333 for (
int c = 0; c < connSize; ++c) {
1334 size_t lid = overlapmap.getLocalElement(elmtConn[c]);
1335 std::vector<int> dofsPerField(numFields_,0);
1336 for (
int n = 0; n <
numFields[c]; ++n) {
1337 int whichField = fieldIds[offset];
1341 localOrdering.push_back(twoview[whichField][lid]+dofsPerField[whichField]);
1343 dofsPerField[whichField]++;
1346 LO thisID=myElements[l];
1347 if(elementGIDs.size()<=(size_t)thisID){
1348 elementGIDs.resize(thisID+1);
1350 elementGIDs[thisID]=localOrdering;
1355 template <
typename LO,
typename GO>
1358 std::vector<std::vector<LO> > elementLIDs(elementGIDs_.size());
1360 std::vector<GO> ownedAndGhosted;
1361 this->getOwnedAndGhostedIndices(ownedAndGhosted);
1364 std::unordered_map<GO,LO> hashMap;
1365 for(std::size_t i = 0; i < ownedAndGhosted.size(); ++i)
1366 hashMap[ownedAndGhosted[i]] = i;
1368 for (std::size_t i = 0; i < elementGIDs_.size(); ++i) {
1369 const std::vector<GO>& gids = elementGIDs_[i];
1370 std::vector<LO>&
lids = elementLIDs[i];
1371 lids.resize(gids.size());
1372 for (std::size_t g = 0; g < gids.size(); ++g)
1373 lids[g] = hashMap[gids[g]];
1376 this->setLocalIds(elementLIDs);
void buildLocalIdsFromOwnedAndGhostedElements()
const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const
void computeCellEdgeOrientations(const std::vector< std::pair< int, int > > &topEdgeIndices, const std::vector< GlobalOrdinalT > &topology, const FieldPattern &fieldPattern, std::vector< signed char > &orientation)
Sums all entries of a Rank 2 Kokkos View.
int getFieldNum(const std::string &string) const
Get the number used for access to this field.
void getElementGIDs(LO localElementID, std::vector< GO > &gids, const std::string &blockIdHint="") const
get associated GIDs for a given local element
int getNumOwned() const
Get the number of indices owned by this processor.
const std::string & getFieldString(int num) const
Reverse lookup of the field string from a field number.
void fillGIDsFromOverlappedMV(const ElementBlockAccess &access, std::vector< std::vector< GO > > &elementGIDs, const Tpetra::Map< LO, GO, panzer::TpetraNodeType > &overlapmap, const Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &overlap_mv) const
#define PANZER_DOFMGR_FUNC_TIME_MONITOR(a)
void setConnManager(const Teuchos::RCP< ConnManager< LO, GO > > &connMngr, MPI_Comm mpiComm)
Adds a Connection Manager that will be associated with this DOFManager.
const std::vector< int > & getGIDFieldOffsets(const std::string &blockID, int fieldNum) const
Teuchos::RCP< const Tpetra::Map< LO, GO, panzer::TpetraNodeType > > buildOverlapMapFromElements(const ElementBlockAccess &access) const
int getNumGhosted() const
Get the number of indices ghosted for this processor.
FieldType
The type of discretization to use for a field pattern.
const std::vector< LO > & getElementBlock(const std::string &eBlock) const
Teuchos::RCP< ConnManager< LocalOrdinalT, GlobalOrdinalT > > resetIndices()
Reset the indices for this DOF manager.
Kokkos::View< const LO **, PHX::Device > lids
virtual int numberIds() const
std::size_t blockIdToIndex(const std::string &blockId) const
Teuchos::RCP< const FieldPattern > getFieldPattern(const std::string &name) const
Find a field pattern stored for a particular block and field number. This will retrive the pattern ad...
int addField(const std::string &str, const Teuchos::RCP< const FieldPattern > &pattern, const panzer::FieldType &type=panzer::FieldType::CG)
Add a field to the DOF manager.
void setFieldOrder(const std::vector< std::string > &fieldOrder)
void getOwnedAndGhostedIndices(std::vector< GlobalOrdinalT > &indices) const
Get the set of owned and ghosted indices for this processor.
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
void getGhostedIndices(std::vector< GlobalOrdinalT > &indices) const
Get the set of indices ghosted for this processor.
void buildUnknownsOrientation()
void buildGlobalUnknowns()
builds the global unknowns array
void computePatternFaceIndices(const FieldPattern &pattern, std::vector< std::vector< int > > &faceIndices)
bool validFieldOrder(const std::vector< std::string > &proposed_fieldOrder)
void getFieldOrder(std::vector< std::string > &fieldOrder) const
void getOwnedIndices(std::vector< GlobalOrdinalT > &indices) const
Get the set of indices owned by this processor.
void computeCellFaceOrientations(const std::vector< std::pair< int, int > > &topEdgeIndices, const std::vector< GlobalOrdinalT > &topology, const FieldPattern &fieldPattern, std::vector< signed char > &orientation)
std::pair< Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > >, Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > > > buildGlobalUnknowns_GUN(const Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &tagged_overlap_mv, Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &overlap_mv) const
Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > > buildTaggedMultiVector(const ElementBlockAccess &access)
void printFieldInformation(std::ostream &os) const
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
void getElementOrientation(LocalOrdinalT localElmtId, std::vector< double > &gidsOrientation) const
Get a vector containg the orientation of the GIDs relative to the neighbors.
int getNumFields() const
gets the number of fields
bool fieldInBlock(const std::string &field, const std::string &block) const
const Kokkos::View< const int *, PHX::Device > getGIDFieldOffsetsKokkos(const std::string &blockID, int fieldNum) const
const std::pair< std::vector< int >, std::vector< int > > & getGIDFieldOffsets_closure(const std::string &blockId, int fieldNum, int subcellDim, int subcellId) const
Use the field pattern so that you can find a particular field in the GIDs array. This version lets yo...
void computePatternEdgeIndices(const FieldPattern &pattern, std::vector< std::pair< int, int > > &edgeIndices)
void ownedIndices(const std::vector< GlobalOrdinalT > &indices, std::vector< bool > &isOwned) const
int getNumOwnedAndGhosted() const
Get the number of owned and ghosted indices for this processor.