Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_DOFManager.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43#ifndef PANZER_DOF_MANAGER2_IMPL_HPP
44#define PANZER_DOF_MANAGER2_IMPL_HPP
45
46#include <map>
47#include <set>
48
49#include <mpi.h>
50
51#include "PanzerDofMgr_config.hpp"
57#include "Panzer_DOFManager.hpp"
59#include "Teuchos_GlobalMPISession.hpp"
61
62#include "Teuchos_RCP.hpp"
63#include "Teuchos_Array.hpp"
64#include "Teuchos_ArrayView.hpp"
65
66#include "Tpetra_Map.hpp"
67#include "Tpetra_Export.hpp"
68#include "Tpetra_Vector.hpp"
69#include "Tpetra_MultiVector.hpp"
70
71#include <unordered_set> // a hash table
72
73/*
74#define HAVE_ZOLTAN2
75#ifdef HAVE_ZOLTAN2
76#include "Zoltan2_XpetraCrsGraphInput.hpp"
77#include "Zoltan2_OrderingProblem.hpp"
78#endif
79*/
80
81namespace panzer {
82
84namespace {
85template <typename LocalOrdinal,typename GlobalOrdinal>
86class HashTieBreak : public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
87 const unsigned int seed_;
88
89public:
90 HashTieBreak(const unsigned int seed = (2654435761U))
91 : seed_(seed) { }
92
93 virtual std::size_t selectedIndex(GlobalOrdinal GID,
94 const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid) const
95 {
96 // this is Epetra's hash/Tpetra's default hash: See Tpetra HashTable
97 int intkey = (int) ((GID & 0x000000007fffffffLL) +
98 ((GID & 0x7fffffff80000000LL) >> 31));
99 return std::size_t((seed_ ^ intkey) % pid_and_lid.size());
100 }
101};
102
103}
104
106namespace {
107template <typename LocalOrdinal,typename GlobalOrdinal>
108class GreedyTieBreak : public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
109
110public:
111 GreedyTieBreak() { }
112
113 virtual bool mayHaveSideEffects() const {
114 return true;
115 }
116
117 virtual std::size_t selectedIndex(GlobalOrdinal /* GID */,
118 const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid) const
119 {
120 // always choose index of pair with smallest pid
121 const std::size_t numLids = pid_and_lid.size();
122 std::size_t idx = 0;
123 int minpid = pid_and_lid[0].first;
124 std::size_t minidx = 0;
125 for (idx = 0; idx < numLids; ++idx) {
126 if (pid_and_lid[idx].first < minpid) {
127 minpid = pid_and_lid[idx].first;
128 minidx = idx;
129 }
130 }
131 return minidx;
132 }
133};
134
135}
136
140
141using Teuchos::RCP;
142using Teuchos::rcp;
143using Teuchos::ArrayRCP;
144using Teuchos::Array;
145using Teuchos::ArrayView;
146
149 : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
150{ }
151
153DOFManager::DOFManager(const Teuchos::RCP<ConnManager> & connMngr,MPI_Comm mpiComm)
154 : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
155{
156 setConnManager(connMngr,mpiComm);
157}
158
160void DOFManager::setConnManager(const Teuchos::RCP<ConnManager> & connMngr, MPI_Comm mpiComm)
161{
162 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
163 "DOFManager::setConnManager: setConnManager cannot be called after "
164 "buildGlobalUnknowns has been called");
165 connMngr_ = connMngr;
166 //We acquire the block ordering from the connmanager.
167 connMngr_->getElementBlockIds(blockOrder_);
168 for (size_t i = 0; i < blockOrder_.size(); ++i) {
169 blockNameToID_.insert(std::map<std::string,int>::value_type(blockOrder_[i],i));
170 //We must also initialize vectors for FP associations.
171 }
172 blockToAssociatedFP_.resize(blockOrder_.size());
173 communicator_ = Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper(mpiComm)));
174}
175
177//Adds a field to be used in creating the Global Numbering
178//Returns the index for the field pattern
179int DOFManager::addField(const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern,
180 const panzer::FieldType& type)
181{
182 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
183 "DOFManager::addField: addField cannot be called after "
184 "buildGlobalUnknowns has been called");
185
186 fieldPatterns_.push_back(pattern);
187 fieldTypes_.push_back(type);
188 fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
189
190 //The default values for IDs are the sequential order they are added in.
191 fieldStringOrder_.push_back(str);
192 fieldAIDOrder_.push_back(numFields_);
193
194 for(std::size_t i=0;i<blockOrder_.size();i++) {
195 blockToAssociatedFP_[i].push_back(numFields_);
196 }
197
198 ++numFields_;
199 return numFields_-1;
200}
201
203int DOFManager::addField(const std::string & blockID, const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern,
204 const panzer::FieldType& type)
205{
206 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
207 "DOFManager::addField: addField cannot be called after "
208 "buildGlobalUnknowns has been called");
209 TEUCHOS_TEST_FOR_EXCEPTION((connMngr_==Teuchos::null),std::logic_error,
210 "DOFManager::addField: you must add a ConnManager before"
211 "you can associate a FP with a given block.")
212 bool found=false;
213 size_t blocknum=0;
214 while(!found && blocknum<blockOrder_.size()){
215 if(blockOrder_[blocknum]==blockID){
216 found=true;
217 break;
218 }
219 blocknum++;
220 }
221 TEUCHOS_TEST_FOR_EXCEPTION(!found,std::logic_error, "DOFManager::addField: Invalid block name.");
222
223 //This will be different if the FieldPattern is already present.
224 //We need to check for that.
225 found=false;
226 std::map<std::string,int>::const_iterator fpIter = fieldNameToAID_.find(str);
227 if(fpIter!=fieldNameToAID_.end())
228 found=true;
229
230 if(!found){
231 fieldPatterns_.push_back(pattern);
232 fieldTypes_.push_back(type);
233 fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
234 //The default values for IDs are the sequential order they are added in.
235 fieldStringOrder_.push_back(str);
236 fieldAIDOrder_.push_back(numFields_);
237
238 //This is going to be associated with blocknum.
239 blockToAssociatedFP_[blocknum].push_back(numFields_);
240 ++numFields_;
241 return numFields_-1;
242 }
243 else{
244 blockToAssociatedFP_[blocknum].push_back(fpIter->second);
245 return numFields_;
246 }
247}
248
250Teuchos::RCP<const FieldPattern> DOFManager::getFieldPattern(const std::string & name) const
251{
252 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(name);
253 if(fitr==fieldNameToAID_.end())
254 return Teuchos::null;
255
256 if(fitr->second<int(fieldPatterns_.size()))
257 return fieldPatterns_[fitr->second];
258
259 return Teuchos::null;
260}
261
263Teuchos::RCP<const FieldPattern> DOFManager::getFieldPattern(const std::string & blockId, const std::string & fieldName) const
264{
265 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(fieldName);
266 if(fitr==fieldNameToAID_.end())
267 return Teuchos::null;
268
269 bool found=false;
270 size_t blocknum=0;
271 while(!found && blocknum<blockOrder_.size()){
272 if(blockOrder_[blocknum]==blockId){
273 found=true;
274 break;
275 }
276 blocknum++;
277 }
278
279 if(!found)
280 return Teuchos::null;
281
282 std::vector<int>::const_iterator itr = std::find(blockToAssociatedFP_[blocknum].begin(),
283 blockToAssociatedFP_[blocknum].end(),fitr->second);
284 if(itr!=blockToAssociatedFP_[blocknum].end()) {
285 if(fitr->second<int(fieldPatterns_.size()))
286 return fieldPatterns_[fitr->second];
287 }
288
289 return Teuchos::null;
290}
291
293void DOFManager::getOwnedIndices(std::vector<panzer::GlobalOrdinal>& indices) const
294{
295 indices = owned_;
296}
297
299void DOFManager::getGhostedIndices(std::vector<panzer::GlobalOrdinal>& indices) const
300{
301 indices = ghosted_;
302}
303
305void DOFManager::getOwnedAndGhostedIndices(std::vector<panzer::GlobalOrdinal>& indices) const
306{
307 using std::size_t;
308 indices.resize(owned_.size() + ghosted_.size());
309 for (size_t i(0); i < owned_.size(); ++i)
310 indices[i] = owned_[i];
311 for (size_t i(0); i < ghosted_.size(); ++i)
312 indices[owned_.size() + i] = ghosted_[i];
313}
314
316void DOFManager::getElementGIDsAsInt(panzer::LocalOrdinal localElementID, std::vector<int>& gids, const std::string& /* blockIdHint */) const
317{
318 const auto & element_ids = elementGIDs_[localElementID];
319 gids.resize(element_ids.size());
320 for (std::size_t i=0; i < gids.size(); ++i)
321 gids[i] = element_ids[i];
322}
323
325void DOFManager::getOwnedIndicesAsInt(std::vector<int>& indices) const
326{
327 indices.resize(owned_.size());
328 for (std::size_t i=0; i < owned_.size(); ++i)
329 indices[i] = owned_[i];
330}
331
333void DOFManager::getGhostedIndicesAsInt(std::vector<int>& indices) const
334{
335 indices.resize(ghosted_.size());
336 for (std::size_t i=0; i < ghosted_.size(); ++i)
337 indices[i] = ghosted_[i];
338}
339
341void DOFManager::getOwnedAndGhostedIndicesAsInt(std::vector<int>& indices) const
342{
343 indices.resize(owned_.size() + ghosted_.size());
344 for (std::size_t i=0; i < owned_.size(); ++i)
345 indices[i] = owned_[i];
346 for (std::size_t i=0; i < ghosted_.size(); ++i)
347 indices[owned_.size() + i] = ghosted_[i];
348}
349
352{
353 return owned_.size();
354}
355
358{
359 return ghosted_.size();
360}
361
364{
365 return owned_.size() + ghosted_.size();
366}
367
370{
371 return numFields_;
372}
373
375const std::vector<int> &
376DOFManager::getGIDFieldOffsets(const std::string & blockID, int fieldNum) const
377{
378 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error, "DOFManager::getGIDFieldOffsets: cannot be called before "
379 "buildGlobalUnknowns has been called");
380 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
381 if(bitr==blockNameToID_.end())
382 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
383 int bid=bitr->second;
384 if(fa_fps_[bid]!=Teuchos::null)
385 return fa_fps_[bid]->localOffsets(fieldNum);
386
387 static const std::vector<int> empty;
388 return empty;
389}
390
392const PHX::View<const int*>
393DOFManager::getGIDFieldOffsetsKokkos(const std::string & blockID, int fieldNum) const
394{
395 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error, "DOFManager::getGIDFieldOffsets: cannot be called before "
396 "buildGlobalUnknowns has been called");
397 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
398 if(bitr==blockNameToID_.end()) {
399 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
400 }
401
402 int bid=bitr->second;
403 if(fa_fps_[bid]!=Teuchos::null)
404 return fa_fps_[bid]->localOffsetsKokkos(fieldNum);
405
406 static const PHX::View<int*> empty("panzer::DOFManager::getGIDFieldOffsetsKokkos() empty",0);
407 return empty;
408}
409
411const PHX::ViewOfViews3<1,PHX::View<const int*>>
412DOFManager::getGIDFieldOffsetsKokkos(const std::string & blockID, const std::vector<int> & fieldNums) const
413{
414 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error, "DOFManager::getGIDFieldOffsets: cannot be called before "
415 "buildGlobalUnknowns has been called");
416 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
417 if(bitr==blockNameToID_.end()) {
418 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
419 }
420
421 PHX::ViewOfViews3<1,PHX::View<const int*>> vov("panzer::getGIDFieldOffsetsKokkos vector version",fieldNums.size());
422 vov.disableSafetyCheck(); // Its going to be moved/copied
423
424 int bid=bitr->second;
425
426 for (size_t i=0; i < fieldNums.size(); ++i) {
427 if(fa_fps_[bid]!=Teuchos::null) {
428 vov.addView(fa_fps_[bid]->localOffsetsKokkos(fieldNums[i]),i);
429 }
430 else {
431 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"panzer::DOFManager::getGIDFieldOffsets() - no field pattern exists in block "
432 << blockID << "for field number " << fieldNums[i] << " exists!");
433 // vov.addView(PHX::View<int*>("panzer::DOFManager::getGIDFieldOffsetsKokkos() empty",0),i);
434 }
435 }
436
437 vov.syncHostToDevice();
438
439 return vov;
440}
441
443void DOFManager::getElementGIDs(panzer::LocalOrdinal localElementID, std::vector<panzer::GlobalOrdinal>& gids, const std::string& /* blockIdHint */) const
444{
445 gids = elementGIDs_[localElementID];
446}
447
450{
451 /* STEPS.
452 * 1. Build GA_FP and all block's FA_FP's and place into respective data structures.
453 */
455 fieldPatterns_.push_back(Teuchos::rcp(new NodalFieldPattern(fieldPatterns_[0]->getCellTopology())));
456 fieldTypes_.push_back(FieldType::CG);
457 }
458
459 TEUCHOS_ASSERT(fieldPatterns_.size() == fieldTypes_.size());
460 std::vector<std::pair<FieldType,Teuchos::RCP<const FieldPattern>>> tmp;
461 for (std::size_t i=0; i < fieldPatterns_.size(); ++i)
462 tmp.push_back(std::make_pair(fieldTypes_[i],fieldPatterns_[i]));
463
464 RCP<GeometricAggFieldPattern> aggFieldPattern = Teuchos::rcp(new GeometricAggFieldPattern(tmp));
465
466 connMngr_->buildConnectivity(*aggFieldPattern);
467
468 // using new geometric pattern, build global unknowns
469 buildGlobalUnknowns(aggFieldPattern);
470}
471
473void DOFManager::buildGlobalUnknowns(const Teuchos::RCP<const FieldPattern> & geomPattern)
474{
475 // some typedefs
476 typedef panzer::TpetraNodeType Node;
477 typedef Tpetra::Map<panzer::LocalOrdinal, panzer::GlobalOrdinal, Node> Map;
478
479 typedef Tpetra::Import<panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> Import;
480
481 //the GIDs are of type panzer::GlobalOrdinal.
482 typedef Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> MultiVector;
483
484 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns",buildGlobalUnknowns);
485
486 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
487 out.setOutputToRootOnly(-1);
488 out.setShowProcRank(true);
489
490 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
491 "DOFManager::buildGlobalUnknowns: buildGlobalUnknowns cannot be called again "
492 "after buildGlobalUnknowns has been called");
493
494 // this is a safety check to make sure that nodes are included
495 // in the geometric field pattern when orientations are required
497 std::size_t sz = geomPattern->getSubcellIndices(0,0).size();
498
499 TEUCHOS_TEST_FOR_EXCEPTION(sz==0,std::logic_error,
500 "DOFManager::buildGlobalUnknowns requires a geometric pattern including "
501 "the nodes when orientations are needed!");
502 }
503
504 /* STEPS.
505 * 1. Build all block's FA_FP's and place into respective data structures.
506 */
507 ga_fp_ = geomPattern;
508
509 // given a set of elements over each element block build an overlap
510 // map that will provide the required element entities for the
511 // set of elements requested.
512 ElementBlockAccess ownedAccess(true,connMngr_);
513
514 // INPUT: To the algorithm in the GUN paper
515 RCP<MultiVector> tagged_overlap_mv = buildTaggedMultiVector(ownedAccess);
516 RCP<const Map> overlap_map = tagged_overlap_mv->getMap();
517
518 RCP<MultiVector> overlap_mv = Tpetra::createMultiVector<panzer::GlobalOrdinal>(overlap_map,(size_t)numFields_);
519
520 // call the GUN paper algorithm
521 auto non_overlap_pair = buildGlobalUnknowns_GUN(*tagged_overlap_mv,*overlap_mv);
522 RCP<MultiVector> non_overlap_mv = non_overlap_pair.first;
523 RCP<MultiVector> tagged_non_overlap_mv = non_overlap_pair.second;
524 RCP<const Map> non_overlap_map = non_overlap_mv->getMap();
525
526 /* 14. Cross reference local element connectivity and overlap map to
527 * create final GID vectors.
528 */
529
530 // this bit of code takes the uniquely assigned GIDs and spreads them
531 // out for processing by local element ID
532 fillGIDsFromOverlappedMV(ownedAccess,elementGIDs_,*overlap_map,*overlap_mv);
533
534 // if neighbor unknowns are required, then make sure they are included
535 // in the elementGIDs_
536 if (useNeighbors_) { // enabling this turns on GID construction for
537 // neighbor processors
538 ElementBlockAccess neighborAccess(false,connMngr_);
539 RCP<const Map> overlap_map_neighbor =
540 buildOverlapMapFromElements(neighborAccess);
541
542 // Export e(overlap_map_neighbor,non_overlap_map);
543 Import imp_neighbor(non_overlap_map,overlap_map_neighbor);
544
545 Teuchos::RCP<MultiVector> overlap_mv_neighbor =
546 Tpetra::createMultiVector<panzer::GlobalOrdinal>(overlap_map_neighbor, (size_t)numFields_);
547
548 // get all neighbor information
549 overlap_mv_neighbor->doImport(*non_overlap_mv, imp_neighbor,
550 Tpetra::REPLACE);
551
553 *overlap_map_neighbor, *overlap_mv_neighbor);
554 }
555
557 // this is where the code is modified to artificially induce GIDs
558 // over 2 Billion unknowns
560 #if 0
561 {
562 panzer::GlobalOrdinal offset = 0xFFFFFFFFLL;
563
564 for (size_t b = 0; b < blockOrder_.size(); ++b) {
565 const std::vector<panzer::LocalOrdinal> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
566 for (std::size_t l = 0; l < myElements.size(); ++l) {
567 std::vector<panzer::GlobalOrdinal> & localGIDs = elementGIDs_[myElements[l]];
568 for(std::size_t c=0;c<localGIDs.size();c++)
569 localGIDs[c] += offset;
570 }
571 }
572
573 Teuchos::ArrayRCP<panzer::GlobalOrdinal> nvals = non_overlap_mv->get1dViewNonConst();
574 for (int j = 0; j < nvals.size(); ++j)
575 nvals[j] += offset;
576 }
577 #endif
578
579 // build owned vector
580 {
581 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns::build_owned_vector",build_owned_vector);
582
583 typedef std::unordered_set<panzer::GlobalOrdinal> HashTable;
584 HashTable isOwned, remainingOwned;
585
586 // owned_ is made up of owned_ids.: This doesn't work for high order
587 auto nvals = non_overlap_mv->getLocalViewHost(Tpetra::Access::ReadOnly);
588 auto tagged_vals = tagged_non_overlap_mv->getLocalViewHost(Tpetra::Access::ReadOnly);
589 TEUCHOS_ASSERT(nvals.size()==tagged_vals.size());
590 for (size_t i = 0; i < nvals.extent(1); ++i) {
591 for (size_t j = 0; j < nvals.extent(0); ++j) {
592 if (nvals(j,i) != -1) {
593 for(panzer::GlobalOrdinal offset=0;offset<tagged_vals(j,i);++offset)
594 isOwned.insert(nvals(j,i)+offset);
595 }
596 else {
597 // sanity check
598 TEUCHOS_ASSERT(tagged_vals(j,i)==0);
599 }
600 }
601 }
602 remainingOwned = isOwned;
603
604 HashTable hashTable; // use to detect if global ID has been added to owned_
605 for (size_t b = 0; b < blockOrder_.size(); ++b) {
606
607 if(fa_fps_[b]==Teuchos::null)
608 continue;
609
610 const std::vector<panzer::LocalOrdinal> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
611
612 for (size_t l = 0; l < myElements.size(); ++l) {
613 const std::vector<panzer::GlobalOrdinal> & localOrdering = elementGIDs_[myElements[l]];
614
615 // add "novel" global ids that are also owned to owned array
616 for(std::size_t i=0;i<localOrdering.size();i++) {
617 // don't try to add if ID is not owned
618 if(isOwned.find(localOrdering[i])==isOwned.end())
619 continue;
620
621 // only add a global ID if it has not been added
622 std::pair<typename HashTable::iterator,bool> insertResult = hashTable.insert(localOrdering[i]);
623 if(insertResult.second) {
624 owned_.push_back(localOrdering[i]);
625 remainingOwned.erase(localOrdering[i]);
626 }
627 }
628 }
629 }
630
631 // add any other owned GIDs that were not already included.
632 // these are ones that are owned locally but not required by any
633 // element owned by this processor (uggh!)
634 for(typename HashTable::const_iterator itr=remainingOwned.begin();itr!=remainingOwned.end();itr++)
635 owned_.push_back(*itr);
636
637 if(owned_.size()!=isOwned.size()) {
638 out << "I'm about to hang because of unknown numbering failure ... sorry! (line = " << __LINE__ << ")" << std::endl;
639 TEUCHOS_TEST_FOR_EXCEPTION(owned_.size()!=isOwned.size(),std::logic_error,
640 "DOFManager::buildGlobalUnkonwns: Failure because not all owned unknowns have been accounted for.");
641 }
642
643 }
644
645 // Build the ghosted_ array. The old simple way led to slow Jacobian
646 // assembly; the new way speeds up Jacobian assembly.
647 {
648 // Loop over all the elements and do a greedy ordering of local values over
649 // the elements for building ghosted_. Hopefully this gives a better
650 // layout for an element-ordered assembly.
651 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns::build_ghosted_array",BGA);
652
653 // Use a hash table to detect if global IDs have been added to owned_.
654 typedef std::unordered_set<panzer::GlobalOrdinal> HashTable;
655 HashTable hashTable;
656 for (std::size_t i = 0; i < owned_.size(); i++)
657 hashTable.insert(owned_[i]);
658
659 // Here we construct an accessor vector, such that we first process
660 // everything in the current element block, optionally followed by
661 // everything in the neighbor element block.
662 std::vector<ElementBlockAccess> blockAccessVec;
663 blockAccessVec.push_back(ElementBlockAccess(true,connMngr_));
664 if(useNeighbors_)
665 blockAccessVec.push_back(ElementBlockAccess(false,connMngr_));
666 for (std::size_t a = 0; a < blockAccessVec.size(); ++a)
667 {
668 // Get the access type (owned or neighbor).
669 const ElementBlockAccess& access = blockAccessVec[a];
670 for (size_t b = 0; b < blockOrder_.size(); ++b)
671 {
672 if (fa_fps_[b] == Teuchos::null)
673 continue;
674 const std::vector<panzer::LocalOrdinal>& myElements =
675 access.getElementBlock(blockOrder_[b]);
676 for (size_t l = 0; l < myElements.size(); ++l)
677 {
678 const std::vector<panzer::GlobalOrdinal>& localOrdering = elementGIDs_[myElements[l]];
679
680 // Add "novel" global IDs into the ghosted_ vector.
681 for (std::size_t i = 0; i < localOrdering.size(); ++i)
682 {
683 std::pair<typename HashTable::iterator, bool> insertResult =
684 hashTable.insert(localOrdering[i]);
685
686 // If the insertion succeeds, then this is "novel" to the owned_
687 // and ghosted_ vectors, so include it in ghosted_.
688 if(insertResult.second)
689 ghosted_.push_back(localOrdering[i]);
690 }
691 }
692 }
693 }
694 }
695
697
698 // build orientations if required
700 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns::build_orientation",BO);
702 }
703
704 // allocate the local IDs
705 if (useNeighbors_) {
706 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns::build_local_ids_from_owned_and_ghosted",BLOFOG);
708 }
709 else {
710 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns::build_local_ids",BLI);
711 this->buildLocalIds();
712 }
713}
714
716std::pair<Teuchos::RCP<Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >,
717 Teuchos::RCP<Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> > >
718DOFManager::buildGlobalUnknowns_GUN(const Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> & tagged_overlap_mv,
719 Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> & overlap_mv) const
720{
721 // some typedefs
722 typedef panzer::TpetraNodeType Node;
723 typedef Tpetra::Map<panzer::LocalOrdinal, panzer::GlobalOrdinal, Node> Map;
724
725 typedef Tpetra::Export<panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> Export;
726 typedef Tpetra::Import<panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> Import;
727
728 //the GIDs are of type panzer::GlobalOrdinal.
729 typedef Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> MultiVector;
730
731 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN",BGU_GUN);
732
733 // LINE 2: In the GUN paper
734 RCP<const Map> overlap_map = tagged_overlap_mv.getMap();
735
736 /* 6. Create a OneToOne map from the overlap map.
737 */
738
739 // LINE 4: In the GUN paper
740
741 RCP<const Map> non_overlap_map;
742 if(!useTieBreak_) {
743 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_04 createOneToOne",GUN04);
744
745 GreedyTieBreak<panzer::LocalOrdinal,panzer::GlobalOrdinal> greedy_tie_break;
746 non_overlap_map = Tpetra::createOneToOne<panzer::LocalOrdinal,panzer::GlobalOrdinal,Node>(overlap_map, greedy_tie_break);
747 }
748 else {
749 // use a hash tie break to get better load balancing from create one to one
750 // Aug. 4, 2016...this is a bad idea and doesn't work
751 HashTieBreak<panzer::LocalOrdinal,panzer::GlobalOrdinal> tie_break;
752 non_overlap_map = Tpetra::createOneToOne<panzer::LocalOrdinal,panzer::GlobalOrdinal,Node>(overlap_map,tie_break);
753 }
754
755 /* 7. Create a non-overlapped multivector from OneToOne map.
756 */
757
758 // LINE 5: In the GUN paper
759
760 Teuchos::RCP<MultiVector> tagged_non_overlap_mv;
761 {
762 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_05 alloc_unique_mv",GUN05);
763
764 tagged_non_overlap_mv = Tpetra::createMultiVector<panzer::GlobalOrdinal>(non_overlap_map,(size_t)numFields_);
765 }
766
767 /* 8. Create an export between the two maps.
768 */
769
770 // LINE 6: In the GUN paper
771 RCP<Export> exp;
772 RCP<Import> imp;
773 RCP<MultiVector> non_overlap_mv;
774 {
775 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_06 export",GUN06);
776
777 exp = rcp(new Export(overlap_map,non_overlap_map));
778
779 /* 9. Export data using ABSMAX.
780 */
781 tagged_non_overlap_mv->doExport(tagged_overlap_mv,*exp,Tpetra::ABSMAX);
782
783 // copy the tagged one, so as to preserve the tagged MV so we can overwrite
784 // the non_overlap_mv
785 non_overlap_mv = rcp(new MultiVector(*tagged_non_overlap_mv,Teuchos::Copy));
786 }
787
788 /* 10. Compute the local sum using Kokkos.
789 */
790
791 // LINES 7-9: In the GUN paper
792
793 panzer::GlobalOrdinal localsum=0;
794 {
795 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_07-09 local_count",GUN07_09);
796 auto values = non_overlap_mv->getLocalViewDevice(Tpetra::Access::ReadOnly);
797 auto mv_size = values.extent(0);
798 Kokkos::parallel_reduce(mv_size,panzer::dof_functors::SumRank2<panzer::GlobalOrdinal,decltype(values)>(values),localsum);
799 }
800
801 /* 11. Create a map using local sums to generate final GIDs.
802 */
803
804 // LINE 10: In the GUN paper
805
806 panzer::GlobalOrdinal myOffset = -1;
807 {
808 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_10 prefix_sum",GUN_10);
809
810 // do a prefix sum
811 panzer::GlobalOrdinal scanResult = 0;
812 Teuchos::scan<int, panzer::GlobalOrdinal> (*getComm(), Teuchos::REDUCE_SUM, static_cast<panzer::GlobalOrdinal> (localsum), Teuchos::outArg (scanResult));
813 myOffset = scanResult - localsum;
814 }
815
816 // LINE 11 and 12: In the GUN paper, these steps are eliminated because
817 // the non_overlap_mv is reused
818
819 /* 12. Iterate through the non-overlapping MV and assign GIDs to
820 * the necessary points. (Assign a -1 elsewhere.)
821 */
822
823 // LINES 13-21: In the GUN paper
824
825 {
826 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_13-21 gid_assignment",GUN13_21);
827 int which_id=0;
828 auto editnonoverlap = non_overlap_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
829 for(size_t i=0; i<non_overlap_mv->getLocalLength(); ++i){
830 for(int j=0; j<numFields_; ++j){
831 if(editnonoverlap(i,j)!=0){
832 // editnonoverlap[j][i]=myOffset+which_id;
833 int ndof = Teuchos::as<int>(editnonoverlap(i,j));
834 editnonoverlap(i,j)=myOffset+which_id;
835 which_id+=ndof;
836 }
837 else{
838 editnonoverlap(i,j)=-1;
839 }
840
841 }
842 }
843 }
844
845 // LINE 22: In the GUN paper. Were performed above, and the overlaped_mv is
846 // abused to handle input tagging.
847
848 /* 13. Import data back to the overlap MV using REPLACE.
849 */
850
851 // LINE 23: In the GUN paper
852
853 {
854 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildGlobalUnknowns_GUN::line_23 final_import",GUN23);
855
856 // use exporter to save on communication setup costs
857 overlap_mv.doImport(*non_overlap_mv,*exp,Tpetra::REPLACE);
858 }
859
860 //std::cout << Teuchos::describe(*non_overlap_mv,Teuchos::VERB_EXTREME) << std::endl;
861
862 // return non_overlap_mv;
863 return std::make_pair(non_overlap_mv,tagged_non_overlap_mv);
864}
865
867Teuchos::RCP<Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >
869{
870 // some typedefs
871 typedef panzer::TpetraNodeType Node;
872 typedef Tpetra::Map<panzer::LocalOrdinal, panzer::GlobalOrdinal, Node> Map;
873 typedef Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,Node> MultiVector;
874
875 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildTaggedMultiVector",BTMV);
876
877 //We will iterate through all of the blocks, building a FieldAggPattern for
878 //each of them.
879
880 for (size_t b = 0; b < blockOrder_.size(); ++b) {
881 std::vector<std::tuple< int, panzer::FieldType, RCP<const panzer::FieldPattern> > > faConstruct;
882 //The ID is going to be the AID, and then everything will work.
883 //The ID should not be the AID, it should be the ID it has in the ordering.
884
885 for (size_t i = 0; i < fieldAIDOrder_.size(); ++i) {
886 int looking = fieldAIDOrder_[i];
887
888 //Check if in b's fp list
889 std::vector<int>::const_iterator reu = std::find(blockToAssociatedFP_[b].begin(), blockToAssociatedFP_[b].end(), looking);
890 if(!(reu==blockToAssociatedFP_[b].end())){
891 faConstruct.push_back(std::make_tuple(i, fieldTypes_[fieldAIDOrder_[i]], fieldPatterns_[fieldAIDOrder_[i]]));
892 }
893
894 }
895
896 if(faConstruct.size()>0) {
897 fa_fps_.push_back(rcp(new FieldAggPattern(faConstruct, ga_fp_)));
898
899 // how many global IDs are in this element block?
900 int gidsInBlock = fa_fps_[fa_fps_.size()-1]->numberIds();
901 elementBlockGIDCount_.push_back(gidsInBlock);
902 }
903 else {
904 fa_fps_.push_back(Teuchos::null);
905 elementBlockGIDCount_.push_back(0);
906 }
907 }
908
909 RCP<const Map> overlapmap = buildOverlapMapFromElements(ownedAccess);
910
911 // LINE 22: In the GUN paper...the overlap_mv is reused for the tagged multivector.
912 // This is a bit of a practical abuse of the algorithm presented in the paper.
913
914 Teuchos::RCP<MultiVector> overlap_mv;
915 {
916 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildTaggedMultiVector::allocate_tagged_multivector",ATMV);
917
918 overlap_mv = Tpetra::createMultiVector<panzer::GlobalOrdinal>(overlapmap,(size_t)numFields_);
919 overlap_mv->putScalar(0); // if tpetra is not initialized with zeros
920 }
921
922 /* 5. Iterate through all local elements again, checking with the FP
923 * information. Mark up the overlap map accordingly.
924 */
925
926 {
927 PANZER_FUNC_TIME_MONITOR_DIFF("panzer::DOFManager::buildTaggedMultiVector::fill_tagged_multivector",FTMV);
928
929 // temporary working vector to fill each row in tagged array
930 std::vector<int> working(overlap_mv->getNumVectors());
931 auto edittwoview_host = overlap_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
932 for (size_t b = 0; b < blockOrder_.size(); ++b) {
933 // there has to be a field pattern associated with the block
934 if(fa_fps_[b]==Teuchos::null)
935 continue;
936
937 const std::vector<panzer::LocalOrdinal> & numFields= fa_fps_[b]->numFieldsPerId();
938 const std::vector<panzer::LocalOrdinal> & fieldIds= fa_fps_[b]->fieldIds();
939 const std::vector<panzer::LocalOrdinal> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
940 for (size_t l = 0; l < myElements.size(); ++l) {
941 auto connSize = connMngr_->getConnectivitySize(myElements[l]);
942 const auto * elmtConn = connMngr_->getConnectivity(myElements[l]);
943 int offset=0;
944 for (int c = 0; c < connSize; ++c) {
945 size_t lid = overlapmap->getLocalElement(elmtConn[c]);
946
947 for(std::size_t i=0;i<working.size();i++)
948 working[i] = 0;
949 for (int n = 0; n < numFields[c]; ++n) {
950 int whichField = fieldIds[offset];
951 //Row will be lid. column will be whichField.
952 //Shove onto local ordering
953 working[whichField]++;
954 offset++;
955 }
956 for(std::size_t i=0;i<working.size();i++) {
957 auto current = edittwoview_host(lid,i);
958 edittwoview_host(lid,i) = (current > working[i]) ? current : working[i];
959 }
960
961 }
962 }
963 }
964 // // verbose output for inspecting overlap_mv
965 // for(int i=0;i<overlap_mv->getLocalLength(); i++) {
966 // for(int j=0;j<overlap_mv->getNumVectors() ; j++)
967 // std::cout << edittwoview[j][i] << " ";
968 // std::cout << std::endl;
969 // }
970 }
971
972 return overlap_mv;
973}
974
976int DOFManager::getFieldNum(const std::string & string) const
977{
978 int ind=0;
979 bool found=false;
980 while(!found && (size_t)ind<fieldStringOrder_.size()){
981 if(fieldStringOrder_[ind]==string)
982 found=true;
983 else
984 ind++;
985 }
986
987 if(found)
988 return ind;
989
990 // didn't find anything...return -1
991 return -1;
992}
993
995void DOFManager::getFieldOrder(std::vector<std::string> & fieldOrder) const
996{
997 fieldOrder.resize(fieldStringOrder_.size());
998 for (size_t i = 0; i < fieldStringOrder_.size(); ++i)
999 fieldOrder[i]=fieldStringOrder_[i];
1000}
1001
1003bool DOFManager::fieldInBlock(const std::string & field, const std::string & block) const
1004{
1005 std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(field);
1006 if(fitr==fieldNameToAID_.end()) {
1007 std::stringstream ss;
1009 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid field name. DOF information is:\n"+ss.str());
1010 }
1011 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(block);
1012 if(bitr==blockNameToID_.end()) {
1013 std::stringstream ss;
1015 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name. DOF information is:\n"+ss.str());
1016 }
1017 int fid=fitr->second;
1018 int bid=bitr->second;
1019
1020 bool found=false;
1021 for (size_t i = 0; i < blockToAssociatedFP_[bid].size(); ++i) {
1022 if(blockToAssociatedFP_[bid][i]==fid){
1023 found=true;
1024 break;
1025 }
1026 }
1027
1028 return found;
1029}
1030
1032const std::vector<int> & DOFManager::getBlockFieldNumbers(const std::string & blockId) const
1033{
1034 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,"DOFManager::getBlockFieldNumbers: BuildConnectivity must be run first.");
1035
1036 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1037 if(bitr==blockNameToID_.end())
1038 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
1039 int bid=bitr->second;
1040
1041 // there has to be a field pattern assocaited with the block
1042 if(fa_fps_[bid]!=Teuchos::null)
1043 return fa_fps_[bid]->fieldIds();
1044
1045 // nothing to return
1046 static std::vector<int> empty;
1047 return empty;
1048}
1049
1051const std::pair<std::vector<int>,std::vector<int> > &
1052DOFManager::getGIDFieldOffsets_closure(const std::string & blockId, int fieldNum, int subcellDim,int subcellId) const
1053{
1054 TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,"DOFManager::getGIDFieldOffsets_closure: BuildConnectivity must be run first.");
1055 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1056 if(bitr==blockNameToID_.end())
1057 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "DOFManager::getGIDFieldOffsets_closure: invalid block name.");
1058
1059 // there has to be a field pattern assocaited with the block
1060 if(fa_fps_[bitr->second]!=Teuchos::null)
1061 return fa_fps_[bitr->second]->localOffsets_closure(fieldNum, subcellDim, subcellId);
1062
1063 static std::pair<std::vector<int>,std::vector<int> > empty;
1064 return empty;
1065}
1066
1068void DOFManager::ownedIndices(const std::vector<panzer::GlobalOrdinal> & indices,std::vector<bool> & isOwned) const
1069{
1070 //Resizes the isOwned array.
1071 if(indices.size()!=isOwned.size())
1072 isOwned.resize(indices.size(),false);
1073 typename std::vector<panzer::GlobalOrdinal>::const_iterator endOf = owned_.end();
1074 for (std::size_t i = 0; i < indices.size(); ++i) {
1075 isOwned[i] = ( std::find(owned_.begin(), owned_.end(), indices[i])!=endOf );
1076 }
1077}
1078
1080void DOFManager::setFieldOrder(const std::vector<std::string> & fieldOrder)
1081{
1082 TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
1083 "DOFManager::setFieldOrder: setFieldOrder cannot be called after "
1084 "buildGlobalUnknowns has been called");
1085 if(validFieldOrder(fieldOrder)){
1086 fieldStringOrder_=fieldOrder;
1087 //We also need to reassign the fieldAIDOrder_.
1088 for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1090 }
1091 }
1092 else {
1093 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::setFieldOrder: Invalid Field Ordering!");
1094 }
1095}
1096
1098bool DOFManager::validFieldOrder(const std::vector<std::string> & proposed_fieldOrder)
1099{
1100 if(fieldStringOrder_.size()!=proposed_fieldOrder.size())
1101 return false;
1102 //I'm using a not very efficient way of doing this, but there should never
1103 //be that many fields, so it isn't that big of a deal.
1104 for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1105 bool found=false;
1106 for (size_t j = 0; j < proposed_fieldOrder.size(); ++j) {
1107 if(fieldStringOrder_[i]==proposed_fieldOrder[j])
1108 found=true;
1109 }
1110 if(!found)
1111 return false;
1112 }
1113 return true;
1114}
1115
1117const std::string & DOFManager::getFieldString(int num) const
1118{
1119 //A reverse lookup through fieldStringOrder_.
1120 if(num>=(int)fieldStringOrder_.size())
1121 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "DOFManager::getFieldString: invalid number");
1122 return fieldStringOrder_[num];
1123}
1124
1126// Everything associated with orientation is not yet built, but this
1127// is the method as found in Panzer_DOFManager_impl.hpp. There are
1128// going to need to be some substantial changes to the code as it applies
1129// to this DOFManager, but the basic ideas and format should be similar.
1130//
1132{
1133 orientation_.clear(); // clean up previous work
1134
1135 std::vector<std::string> elementBlockIds;
1136 connMngr_->getElementBlockIds(elementBlockIds);
1137
1138 // figure out how many total elements are owned by this processor (why is this so hard!)
1139 std::size_t myElementCount = 0;
1140 for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin(); blockItr!=elementBlockIds.end();++blockItr)
1141 myElementCount += connMngr_->getElementBlock(*blockItr).size();
1142
1143 // allocate for each block
1144 orientation_.resize(myElementCount);
1145
1146 // loop over all element blocks
1147 for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin();
1148 blockItr!=elementBlockIds.end();++blockItr) {
1149 const std::string & blockName = *blockItr;
1150
1151 // this block has no unknowns (or elements)
1152 std::map<std::string,int>::const_iterator fap = blockNameToID_.find(blockName);
1153 if(fap==blockNameToID_.end()) {
1154 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::buildUnknownsOrientation: invalid block name");
1155 }
1156
1157 int bid=fap->second;
1158
1159 if(fa_fps_[bid]==Teuchos::null)
1160 continue;
1161
1162 // grab field patterns, will be necessary to compute orientations
1163 const FieldPattern & fieldPattern = *fa_fps_[bid];
1164
1165 //Should be ga_fp_ (geometric aggregate field pattern)
1166 std::vector<std::pair<int,int> > topEdgeIndices;
1168
1169 // grab face orientations if 3D
1170 std::vector<std::vector<int> > topFaceIndices;
1171 if(ga_fp_->getDimension()==3)
1173
1174 //How many GIDs are associated with a particular element bloc
1175 const std::vector<panzer::LocalOrdinal> & elmts = connMngr_->getElementBlock(blockName);
1176 for(std::size_t e=0;e<elmts.size();e++) {
1177 // this is the vector of orientations to fill: initialize it correctly
1178 std::vector<signed char> & eOrientation = orientation_[elmts[e]];
1179
1180 // This resize seems to be the same as fieldPattern.numberIDs().
1181 // When computer edge orientations is called, that is the assert.
1182 // There should be no reason to make it anymore complicated.
1183 eOrientation.resize(fieldPattern.numberIds());
1184 for(std::size_t s=0;s<eOrientation.size();s++)
1185 eOrientation[s] = 1; // put in 1 by default
1186
1187 // get geometry ids
1188 auto connSz = connMngr_->getConnectivitySize(elmts[e]);
1189 const panzer::GlobalOrdinal * connPtr = connMngr_->getConnectivity(elmts[e]);
1190 const std::vector<panzer::GlobalOrdinal> connectivity(connPtr,connPtr+connSz);
1191
1192 orientation_helpers::computeCellEdgeOrientations(topEdgeIndices, connectivity, fieldPattern, eOrientation);
1193
1194 // compute face orientations in 3D
1195 if(ga_fp_->getDimension()==3)
1196 orientation_helpers::computeCellFaceOrientations(topFaceIndices, connectivity, fieldPattern, eOrientation);
1197 }
1198 }
1199}
1200
1202void DOFManager::getElementOrientation(panzer::LocalOrdinal localElmtId,std::vector<double> & gidsOrientation) const
1203{
1204 TEUCHOS_TEST_FOR_EXCEPTION(orientation_.size()==0,std::logic_error,
1205 "DOFManager::getElementOrientations: Orientations were not constructed!");
1206
1207 const std::vector<signed char> & local_o = orientation_[localElmtId];
1208 gidsOrientation.resize(local_o.size());
1209 for(std::size_t i=0;i<local_o.size();i++) {
1210 gidsOrientation[i] = double(local_o[i]);
1211 }
1212}
1213
1214Teuchos::RCP<ConnManager> DOFManager::resetIndices()
1215{
1216 Teuchos::RCP<ConnManager> connMngr = connMngr_;
1217
1218 connMngr_ = Teuchos::null;
1219
1220 // wipe out FEI objects
1221 orientation_.clear(); // clean up previous work
1222 fa_fps_.clear();
1223 elementGIDs_.clear();
1224 owned_.clear();
1225 ghosted_.clear();
1226 elementBlockGIDCount_.clear();
1227
1228 return connMngr;
1229}
1230
1232std::size_t DOFManager::blockIdToIndex(const std::string & blockId) const
1233{
1234 std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1235 if(bitr==blockNameToID_.end())
1236 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
1237 return bitr->second;
1238}
1239
1241void DOFManager::printFieldInformation(std::ostream & os) const
1242{
1243 os << "DOFManager Field Information: " << std::endl;
1244
1245 // sanity check
1246 TEUCHOS_ASSERT(blockOrder_.size()==blockToAssociatedFP_.size());
1247
1248 for(std::size_t i=0;i<blockOrder_.size();i++) {
1249 os << " Element Block = " << blockOrder_[i] << std::endl;
1250
1251 // output field information
1252 const std::vector<int> & fieldIds = blockToAssociatedFP_[i];
1253 for(std::size_t f=0;f<fieldIds.size();f++)
1254 os << " \"" << getFieldString(fieldIds[f]) << "\" is field ID " << fieldIds[f] << std::endl;
1255 }
1256}
1257
1259Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >
1261{
1262 PANZER_FUNC_TIME_MONITOR("panzer::DOFManager::builderOverlapMapFromElements");
1263
1264 /*
1265 * 2. Iterate through all local elements and create the overlapVector
1266 * of concerned elements.
1267 */
1268
1269 std::set<panzer::GlobalOrdinal> overlapset;
1270 for (size_t i = 0; i < blockOrder_.size(); ++i) {
1271 const std::vector<panzer::LocalOrdinal> & myElements = access.getElementBlock(blockOrder_[i]);
1272 for (size_t e = 0; e < myElements.size(); ++e) {
1273 auto connSize = connMngr_->getConnectivitySize(myElements[e]);
1274 const panzer::GlobalOrdinal * elmtConn = connMngr_->getConnectivity(myElements[e]);
1275 for (int k = 0; k < connSize; ++k) {
1276 overlapset.insert(elmtConn[k]);
1277 }
1278 }
1279 }
1280
1281 Array<panzer::GlobalOrdinal> overlapVector;
1282 for (typename std::set<panzer::GlobalOrdinal>::const_iterator itr = overlapset.begin(); itr!=overlapset.end(); ++itr) {
1283 overlapVector.push_back(*itr);
1284 }
1285
1286 /* 3. Construct an overlap map from this structure.
1287 */
1288 return Tpetra::createNonContigMapWithNode<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>(overlapVector,getComm());
1289}
1290
1294 std::vector<std::vector< panzer::GlobalOrdinal > > & elementGIDs,
1295 const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> & overlapmap,
1296 const Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> & const_overlap_mv) const
1297{
1298 using Teuchos::ArrayRCP;
1299
1300 //To generate elementGIDs we need to go through all of the local elements.
1301 auto overlap_mv = const_cast<Tpetra::MultiVector<panzer::GlobalOrdinal,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>&>(const_overlap_mv);
1302 const auto twoview_host = overlap_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
1303
1304 //And for each of the things in fa_fp.fieldIds we go to that column. To the the row,
1305 //we move from globalID to localID in the map and use our local value for something.
1306 for (size_t b = 0; b < blockOrder_.size(); ++b) {
1307 const std::vector<panzer::LocalOrdinal> & myElements = access.getElementBlock(blockOrder_[b]);
1308
1309 if(fa_fps_[b]==Teuchos::null) {
1310 // fill elements that are not used with empty vectors
1311 for (size_t l = 0; l < myElements.size(); ++l) {
1312 panzer::LocalOrdinal thisID=myElements[l];
1313 if(elementGIDs.size()<=(size_t)thisID)
1314 elementGIDs.resize(thisID+1);
1315 }
1316 continue;
1317 }
1318
1319 const std::vector<int> & numFields= fa_fps_[b]->numFieldsPerId();
1320 const std::vector<int> & fieldIds= fa_fps_[b]->fieldIds();
1321 //
1322 //
1323 for (size_t l = 0; l < myElements.size(); ++l) {
1324 auto connSize = connMngr_->getConnectivitySize(myElements[l]);
1325 const panzer::GlobalOrdinal * elmtConn = connMngr_->getConnectivity(myElements[l]);
1326 std::vector<panzer::GlobalOrdinal> localOrdering;
1327 int offset=0;
1328 for (int c = 0; c < connSize; ++c) {
1329 size_t lid = overlapmap.getLocalElement(elmtConn[c]);
1330 std::vector<int> dofsPerField(numFields_,0);
1331 for (int n = 0; n < numFields[c]; ++n) {
1332 int whichField = fieldIds[offset];
1333 offset++;
1334 //Row will be lid. column will be whichField.
1335 //Shove onto local ordering
1336 localOrdering.push_back(twoview_host(lid,whichField)+dofsPerField[whichField]);
1337
1338 dofsPerField[whichField]++;
1339 }
1340 }
1341 panzer::LocalOrdinal thisID=myElements[l];
1342 if(elementGIDs.size()<=(size_t)thisID){
1343 elementGIDs.resize(thisID+1);
1344 }
1345 elementGIDs[thisID]=localOrdering;
1346 }
1347 }
1348}
1349
1351{
1352 std::vector<std::vector<panzer::LocalOrdinal> > elementLIDs(elementGIDs_.size());
1353
1354 std::vector<panzer::GlobalOrdinal> ownedAndGhosted;
1355 this->getOwnedAndGhostedIndices(ownedAndGhosted);
1356
1357 // build global to local hash map (temporary and used only once)
1358 std::unordered_map<panzer::GlobalOrdinal,panzer::LocalOrdinal> hashMap;
1359 for(std::size_t i = 0; i < ownedAndGhosted.size(); ++i)
1360 hashMap[ownedAndGhosted[i]] = i;
1361
1362 for (std::size_t i = 0; i < elementGIDs_.size(); ++i) {
1363 const std::vector<panzer::GlobalOrdinal>& gids = elementGIDs_[i];
1364 std::vector<panzer::LocalOrdinal>& lids = elementLIDs[i];
1365 lids.resize(gids.size());
1366 for (std::size_t g = 0; g < gids.size(); ++g)
1367 lids[g] = hashMap[gids[g]];
1368 }
1369
1370 this->setLocalIds(elementLIDs);
1371}
1372
1373/*
1374template <typename panzer::LocalOrdinal,typename panzer::GlobalOrdinal>
1375Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> >
1376DOFManager<panzer::LocalOrdinal,panzer::GlobalOrdinal>::runLocalRCMReordering(const Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> > & map)
1377{
1378 typedef panzer::TpetraNodeType Node;
1379 typedef Tpetra::Map<panzer::LocalOrdinal, panzer::GlobalOrdinal, Node> Map;
1380 typedef Tpetra::CrsGraph<panzer::LocalOrdinal, panzer::GlobalOrdinal, Node> Graph;
1381
1382 Teuchos::RCP<Graph> graph = Teuchos::rcp(new Graph(map,0));
1383
1384 // build Crs Graph from the mesh
1385 for (size_t b = 0; b < blockOrder_.size(); ++b) {
1386 if(fa_fps_[b]==Teuchos::null)
1387 continue;
1388
1389 const std::vector<panzer::LocalOrdinal> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
1390 for (size_t l = 0; l < myElements.size(); ++l) {
1391 panzer::LocalOrdinal connSize = connMngr_->getConnectivitySize(myElements[l]);
1392 const panzer::GlobalOrdinal * elmtConn = connMngr_->getConnectivity(myElements[l]);
1393 for (int c = 0; c < connSize; ++c) {
1394 panzer::LocalOrdinal lid = map->getLocalElement(elmtConn[c]);
1395 if(Teuchos::OrdinalTraits<panzer::LocalOrdinal>::invalid()!=lid)
1396 continue;
1397
1398 graph->insertGlobalIndices(elmtConn[c],Teuchos::arrayView(elmtConn,connSize));
1399 }
1400 }
1401 }
1402
1403 graph->fillComplete();
1404
1405 std::vector<panzer::GlobalOrdinal> newOrder(map->getLocalNumElements());
1406 {
1407 // graph is constructed, now run RCM using zoltan2
1408 typedef Zoltan2::XpetraCrsGraphInput<Graph> SparseGraphAdapter;
1409
1410 Teuchos::ParameterList params;
1411 params.set("order_method", "rcm");
1412 SparseGraphAdapter adapter(graph);
1413
1414 Zoltan2::OrderingProblem<SparseGraphAdapter> problem(&adapter, &params,MPI_COMM_SELF);
1415 problem.solve();
1416
1417 // build a new global ording array using permutation
1418 Zoltan2::OrderingSolution<panzer::GlobalOrdinal,panzer::LocalOrdinal> * soln = problem.getSolution();
1419
1420 size_t dummy;
1421 size_t checkLength = soln->getPermutationSize();
1422 panzer::LocalOrdinal * checkPerm = soln->getPermutation(&dummy);
1423
1424 Teuchos::ArrayView<const panzer::GlobalOrdinal > oldOrder = map->getLocalElementList();
1425 TEUCHOS_ASSERT(checkLength==oldOrder.size());
1426 TEUCHOS_ASSERT(checkLength==newOrder.size());
1427
1428 for(size_t i=0;i<checkLength;i++)
1429 newOrder[checkPerm[i]] = oldOrder[i];
1430 }
1431
1432 return Tpetra::createNonContigMap<panzer::LocalOrdinal,panzer::GlobalOrdinal>(newOrder,communicator_);
1433}
1434*/
1435
1436} /*panzer*/
1437
1438#endif
int numFields
const unsigned int seed_
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.
Kokkos::View< const LO **, Kokkos::LayoutRight, PHX::Device > lids
const std::vector< panzer::LocalOrdinal > & getElementBlock(const std::string &eBlock) const
void ownedIndices(const std::vector< panzer::GlobalOrdinal > &indices, std::vector< bool > &isOwned) const
std::vector< int > elementBlockGIDCount_
bool validFieldOrder(const std::vector< std::string > &proposed_fieldOrder)
Teuchos::RCP< Teuchos::Comm< int > > getComm() const
const PHX::View< const int * > getGIDFieldOffsetsKokkos(const std::string &blockID, int fieldNum) const
int getNumGhosted() const
Get the number of indices ghosted for this processor.
std::vector< std::string > blockOrder_
std::vector< panzer::GlobalOrdinal > ghosted_
std::vector< FieldType > fieldTypes_
std::map< std::string, int > blockNameToID_
std::vector< int > fieldAIDOrder_
void getOwnedIndicesAsInt(std::vector< int > &indices) const
Get the set of indices owned by this processor.
void getGhostedIndices(std::vector< panzer::GlobalOrdinal > &indices) const
Get the set of indices ghosted for this processor.
void getFieldOrder(std::vector< std::string > &fieldOrder) const
bool getOrientationsRequired() const
void setConnManager(const Teuchos::RCP< ConnManager > &connMngr, MPI_Comm mpiComm)
Adds a Connection Manager that will be associated with this DOFManager.
std::vector< panzer::GlobalOrdinal > owned_
void fillGIDsFromOverlappedMV(const ElementBlockAccess &access, std::vector< std::vector< panzer::GlobalOrdinal > > &elementGIDs, const Tpetra::Map< panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > &overlapmap, const Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > &overlap_mv) const
void getElementOrientation(panzer::LocalOrdinal localElmtId, std::vector< double > &gidsOrientation) const
Get a vector containg the orientation of the GIDs relative to the neighbors.
Teuchos::RCP< const panzer::FieldPattern > ga_fp_
void buildLocalIdsFromOwnedAndGhostedElements()
Teuchos::RCP< ConnManager > connMngr_
void getOwnedAndGhostedIndicesAsInt(std::vector< int > &indices) const
Get the set of owned and ghosted indices for this processor.
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...
std::pair< Teuchos::RCP< Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > >, Teuchos::RCP< Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > > > buildGlobalUnknowns_GUN(const Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > &tagged_overlap_mv, Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > &overlap_mv) const
std::vector< std::vector< int > > blockToAssociatedFP_
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 setFieldOrder(const std::vector< std::string > &fieldOrder)
Teuchos::RCP< Tpetra::MultiVector< panzer::GlobalOrdinal, panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > > buildTaggedMultiVector(const ElementBlockAccess &access)
std::vector< std::vector< panzer::GlobalOrdinal > > elementGIDs_
std::vector< std::vector< signed char > > orientation_
Teuchos::RCP< ConnManager > resetIndices()
Reset the indices for this DOF manager.
std::size_t blockIdToIndex(const std::string &blockId) const
void getElementGIDsAsInt(panzer::LocalOrdinal localElementID, std::vector< int > &gids, const std::string &blockIdHint="") const
Get the global IDs for a particular element. This function overwrites the gids variable.
void buildGlobalUnknowns()
builds the global unknowns array
Teuchos::RCP< Teuchos::Comm< int > > communicator_
std::vector< Teuchos::RCP< panzer::FieldAggPattern > > fa_fps_
void getOwnedIndices(std::vector< panzer::GlobalOrdinal > &indices) const
Get the set of indices owned by this processor.
int getNumOwnedAndGhosted() const
Get the number of owned and ghosted indices for this processor.
void getElementGIDs(panzer::LocalOrdinal localElementID, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const
get associated GIDs for a given local element
std::map< std::string, int > fieldNameToAID_
std::vector< std::string > fieldStringOrder_
void printFieldInformation(std::ostream &os) const
void getGhostedIndicesAsInt(std::vector< int > &indices) const
Get the set of indices ghosted for this processor.
int getFieldNum(const std::string &string) const
Get the number used for access to this field.
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.
std::vector< Teuchos::RCP< const FieldPattern > > fieldPatterns_
bool fieldInBlock(const std::string &field, const std::string &block) const
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.
int getNumFields() const
gets the number of fields
void getOwnedAndGhostedIndices(std::vector< panzer::GlobalOrdinal > &indices) const
Get the set of owned and ghosted indices for this processor.
Teuchos::RCP< const Tpetra::Map< panzer::LocalOrdinal, panzer::GlobalOrdinal, panzer::TpetraNodeType > > buildOverlapMapFromElements(const ElementBlockAccess &access) const
const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const
const std::vector< int > & getGIDFieldOffsets(const std::string &blockID, int fieldNum) const
virtual int numberIds() const
void setLocalIds(const std::vector< std::vector< panzer::LocalOrdinal > > &localIDs)
void computeCellEdgeOrientations(const std::vector< std::pair< int, int > > &topEdgeIndices, const std::vector< panzer::GlobalOrdinal > &topology, const FieldPattern &fieldPattern, std::vector< signed char > &orientation)
void computePatternEdgeIndices(const FieldPattern &pattern, std::vector< std::pair< int, int > > &edgeIndices)
void computeCellFaceOrientations(const std::vector< std::vector< int > > &topFaceIndices, const std::vector< panzer::GlobalOrdinal > &topology, const FieldPattern &fieldPattern, std::vector< signed char > &orientation)
void computePatternFaceIndices(const FieldPattern &pattern, std::vector< std::vector< int > > &faceIndices)
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
FieldType
The type of discretization to use for a field pattern.
Sums all entries of a Rank 2 Kokkos View.