Panzer  Version of the Day
Panzer_DOFManager_impl.hpp
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"
52 #include "Panzer_FieldPattern.hpp"
55 #include "Panzer_ConnManager.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 #define PANZER_DOFMGR_FUNC_TIME_MONITOR(a) \
74  PANZER_FUNC_TIME_MONITOR(a)
75 
76 #ifdef PHX_KOKKOS_DEVICE_TYPE_CUDA
77 #define PANZER_DOFMGR_REQUIRE_CUDA
78 #endif
79 
80 /*
81 #define HAVE_ZOLTAN2
82 #ifdef HAVE_ZOLTAN2
83 #include "Zoltan2_XpetraCrsGraphInput.hpp"
84 #include "Zoltan2_OrderingProblem.hpp"
85 #endif
86 */
87 
88 namespace panzer {
89 
90 namespace {
91 template <typename LocalOrdinal,typename GlobalOrdinal>
92 class HashTieBreak : public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
93  const unsigned int seed_;
94 
95 public:
96  HashTieBreak(const unsigned int seed = (2654435761U))
97  : seed_(seed) { }
98 
99  virtual std::size_t selectedIndex(GlobalOrdinal GID,
100  const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid) const
101  {
102  // this is Epetra's hash/Tpetra's default hash: See Tpetra HashTable
103  int intkey = (int) ((GID & 0x000000007fffffffLL) +
104  ((GID & 0x7fffffff80000000LL) >> 31));
105  return std::size_t((seed_ ^ intkey) % pid_and_lid.size());
106  }
107 };
108 
109 }
110 
111 
112 namespace {
113 template <typename LocalOrdinal,typename GlobalOrdinal>
114 class GreedyTieBreak : public Tpetra::Details::TieBreak<LocalOrdinal,GlobalOrdinal> {
115 
116 public:
117  GreedyTieBreak() { }
118 
119  virtual bool mayHaveSideEffects() const {
120  return true;
121  }
122 
123  virtual std::size_t selectedIndex(GlobalOrdinal /* GID */,
124  const std::vector<std::pair<int,LocalOrdinal> > & pid_and_lid) const
125  {
126  // always choose index of pair with smallest pid
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;
134  minidx = idx;
135  }
136  }
137  return minidx;
138  }
139 };
140 
141 }
142 
143 
144 using Teuchos::RCP;
145 using Teuchos::rcp;
146 using Teuchos::ArrayRCP;
147 using Teuchos::Array;
148 using Teuchos::ArrayView;
149 
150 template <typename LO, typename GO>
152  : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
153 { }
154 
155 template <typename LO, typename GO>
157  : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
158 {
159  setConnManager(connMngr,mpiComm);
160 }
161 
162 template <typename LO, typename GO>
163 void DOFManager<LO,GO>::setConnManager(const Teuchos::RCP<ConnManager<LO,GO> > & connMngr, MPI_Comm mpiComm)
164 {
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;
169  //We acquire the block ordering from the connmanager.
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));
173  //We must also initialize vectors for FP associations.
174  }
175  blockToAssociatedFP_.resize(blockOrder_.size());
176  communicator_ = Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper(mpiComm)));
177 }
178 
179 
180 //Adds a field to be used in creating the Global Numbering
181 //Returns the index for the field pattern
182 template <typename LO, typename GO>
183 int DOFManager<LO,GO>::addField(const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern,
184  const panzer::FieldType& type)
185 {
186  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
187  "DOFManager::addField: addField cannot be called after "
188  "buildGlobalUnknowns has been called");
189 
190  fieldPatterns_.push_back(pattern);
191  fieldTypes_.push_back(type);
192  fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
193 
194  //The default values for IDs are the sequential order they are added in.
195  fieldStringOrder_.push_back(str);
196  fieldAIDOrder_.push_back(numFields_);
197 
198  for(std::size_t i=0;i<blockOrder_.size();i++) {
199  blockToAssociatedFP_[i].push_back(numFields_);
200  }
201 
202  ++numFields_;
203  return numFields_-1;
204 }
205 
206 template <typename LO, typename GO>
207 int DOFManager<LO,GO>::addField(const std::string & blockID, const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern,
208  const panzer::FieldType& type)
209 {
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.")
216  bool found=false;
217  size_t blocknum=0;
218  while(!found && blocknum<blockOrder_.size()){
219  if(blockOrder_[blocknum]==blockID){
220  found=true;
221  break;
222  }
223  blocknum++;
224  }
225  TEUCHOS_TEST_FOR_EXCEPTION(!found,std::logic_error, "DOFManager::addField: Invalid block name.");
226 
227  //This will be different if the FieldPattern is already present.
228  //We need to check for that.
229  found=false;
230  std::map<std::string,int>::const_iterator fpIter = fieldNameToAID_.find(str);
231  if(fpIter!=fieldNameToAID_.end())
232  found=true;
233 
234  if(!found){
235  fieldPatterns_.push_back(pattern);
236  fieldTypes_.push_back(type);
237  fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
238  //The default values for IDs are the sequential order they are added in.
239  fieldStringOrder_.push_back(str);
240  fieldAIDOrder_.push_back(numFields_);
241 
242  //This is going to be associated with blocknum.
243  blockToAssociatedFP_[blocknum].push_back(numFields_);
244  ++numFields_;
245  return numFields_-1;
246  }
247  else{
248  blockToAssociatedFP_[blocknum].push_back(fpIter->second);
249  return numFields_;
250  }
251 }
252 
253 
254 
255  //Returns the fieldpattern of the given name
256 template <typename LO, typename GO>
258 {
259  std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(name);
260  if(fitr==fieldNameToAID_.end())
261  return Teuchos::null;
262 
263  if(fitr->second<int(fieldPatterns_.size()))
264  return fieldPatterns_[fitr->second];
265 
266  return Teuchos::null;
267 }
268 
269 //Returns the fieldpattern of the given name
270 template <typename LO, typename GO>
271 Teuchos::RCP<const FieldPattern> DOFManager<LO,GO>::getFieldPattern(const std::string & blockId, const std::string & fieldName) const
272 {
273  std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(fieldName);
274  if(fitr==fieldNameToAID_.end())
275  return Teuchos::null;
276 
277  bool found=false;
278  size_t blocknum=0;
279  while(!found && blocknum<blockOrder_.size()){
280  if(blockOrder_[blocknum]==blockId){
281  found=true;
282  break;
283  }
284  blocknum++;
285  }
286 
287  if(!found)
288  return Teuchos::null;
289 
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];
295  }
296 
297  return Teuchos::null;
298 }
299 
301 //
302 // getOwnedIndices()
303 //
305 template<typename LO, typename GO>
306 void
309  std::vector<GO>& indices) const
310 {
311  indices = owned_;
312 } // end of getOwnedIndices()
313 
315 //
316 // getGhostedIndices()
317 //
319 template<typename LO, typename GO>
320 void
323  std::vector<GO>& indices) const
324 {
325  indices = ghosted_;
326 } // end of getGhostedIndices()
327 
329 //
330 // getOwnedAndGhostedIndices()
331 //
333 template<typename LO, typename GO>
334 void
337  std::vector<GO>& indices) const
338 {
339  using std::size_t;
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];
345 } // end of getOwnedAndGhostedIndices()
346 
348 //
349 // getNumOwned()
350 //
352 template<typename LO, typename GO>
353 int
355 getNumOwned() const
356 {
357  return owned_.size();
358 } // end of getNumOwned()
359 
361 //
362 // getNumGhosted()
363 //
365 template<typename LO, typename GO>
366 int
369 {
370  return ghosted_.size();
371 } // end of getNumGhosted()
372 
374 //
375 // getNumOwnedAndGhosted()
376 //
378 template<typename LO, typename GO>
379 int
382 {
383  return owned_.size() + ghosted_.size();
384 } // end of getNumOwnedAndGhosted()
385 
386  //gets the number of fields
387 template <typename LO, typename GO>
389 {
390  return numFields_;
391 }
392 
393 template <typename LO, typename GO>
394 const std::vector<int> & DOFManager<LO,GO>::getGIDFieldOffsets(const std::string & blockID, int fieldNum) const
395 {
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);
404 
405  static const std::vector<int> empty;
406  return empty;
407 }
408 
409 template <typename LO, typename GO>
410 const Kokkos::View<const int*,PHX::Device> DOFManager<LO,GO>::
411 getGIDFieldOffsetsKokkos(const std::string & blockID, int fieldNum) const
412 {
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);
421 
422  static const Kokkos::View<int*,PHX::Device> empty("panzer::DOFManager::getGIDFieldOffsetsKokkos() empty",0);
423  return empty;
424 }
425 
426 template <typename LO, typename GO>
427 void DOFManager<LO,GO>::getElementGIDs(LO localElementID, std::vector<GO>& gids, const std::string& /* blockIdHint */) const
428 {
429  gids = elementGIDs_[localElementID];
430 }
431 
432 template <typename LocalOrdinalT,typename GlobalOrdinalT>
434 {
435  /* STEPS.
436  * 1. Build GA_FP and all block's FA_FP's and place into respective data structures.
437  */
438  if(requireOrientations_){
439  fieldPatterns_.push_back(Teuchos::rcp(new NodalFieldPattern(fieldPatterns_[0]->getCellTopology())));
440  fieldTypes_.push_back(FieldType::CG);
441  }
442 
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]));
447 
449 
450  connMngr_->buildConnectivity(*aggFieldPattern);
451 
452  // using new geometric pattern, build global unknowns
453  buildGlobalUnknowns(aggFieldPattern);
454 }
455 
456 //builds the global unknowns array
457 template <typename LO, typename GO>
459 {
460  // some typedefs
461  typedef panzer::TpetraNodeType Node;
462  typedef Tpetra::Map<LO, GO, Node> Map;
463 
464  typedef Tpetra::Import<LO,GO,Node> Import;
465 
466  //the GIDs are of type GO.
467  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
468 
469  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns");
470 
471  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
472  out.setOutputToRootOnly(-1);
473  out.setShowProcRank(true);
474 
475  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
476  "DOFManager::buildGlobalUnknowns: buildGlobalUnknowns cannot be called again "
477  "after buildGlobalUnknowns has been called");
478 
479  // this is a safety check to make sure that nodes are included
480  // in the geometric field pattern when orientations are required
481  if(getOrientationsRequired()) {
482  std::size_t sz = geomPattern->getSubcellIndices(0,0).size();
483 
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!");
487  }
488 
489  /* STEPS.
490  * 1. Build all block's FA_FP's and place into respective data structures.
491  */
492  ga_fp_ = geomPattern;
493 
494  // given a set of elements over each element block build an overlap
495  // map that will provide the required element entities for the
496  // set of elements requested.
497  ElementBlockAccess ownedAccess(true,connMngr_);
498 
499  // INPUT: To the algorithm in the GUN paper
500  RCP<MultiVector> tagged_overlap_mv = buildTaggedMultiVector(ownedAccess);
501  RCP<const Map> overlap_map = tagged_overlap_mv->getMap();
502 
503  RCP<MultiVector> overlap_mv = Tpetra::createMultiVector<GO>(overlap_map,(size_t)numFields_);
504 
505  // call the GUN paper algorithm
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();
510 
511  /* 14. Cross reference local element connectivity and overlap map to
512  * create final GID vectors.
513  */
514 
515  // this bit of code takes the uniquely assigned GIDs and spreads them
516  // out for processing by local element ID
517  fillGIDsFromOverlappedMV(ownedAccess,elementGIDs_,*overlap_map,*overlap_mv);
518 
519  // if neighbor unknowns are required, then make sure they are included
520  // in the elementGIDs_
521  if (useNeighbors_) { // enabling this turns on GID construction for
522  // neighbor processors
523  ElementBlockAccess neighborAccess(false,connMngr_);
524  RCP<const Map> overlap_map_neighbor =
525  buildOverlapMapFromElements(neighborAccess);
526 
527  // Export e(overlap_map_neighbor,non_overlap_map);
528  Import imp_neighbor(non_overlap_map,overlap_map_neighbor);
529 
530  Teuchos::RCP<MultiVector> overlap_mv_neighbor =
531  Tpetra::createMultiVector<GO>(overlap_map_neighbor, (size_t)numFields_);
532 
533  // get all neighbor information
534  overlap_mv_neighbor->doImport(*non_overlap_mv, imp_neighbor,
535  Tpetra::REPLACE);
536 
537  fillGIDsFromOverlappedMV(neighborAccess, elementGIDs_,
538  *overlap_map_neighbor, *overlap_mv_neighbor);
539  }
540 
542  // this is where the code is modified to artificially induce GIDs
543  // over 2 Billion unknowns
545  #if 0
546  {
547  panzer::Ordinal64 offset = 0xFFFFFFFFLL;
548 
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;
555  }
556  }
557 
558  Teuchos::ArrayRCP<GO> nvals = non_overlap_mv->get1dViewNonConst();
559  for (int j = 0; j < nvals.size(); ++j)
560  nvals[j] += offset;
561  }
562  #endif
563 
564  // build owned vector
565  {
566  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_owned_vector");
567 
568  typedef std::unordered_set<GO> HashTable;
569  HashTable isOwned, remainingOwned;
570 
571  // owned_ is made up of owned_ids.: This doesn't work for high order
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);
579  }
580  else {
581  // sanity check
582  TEUCHOS_ASSERT(tagged_vals[j]==0)
583  }
584  }
585  remainingOwned = isOwned;
586 
587  HashTable hashTable; // use to detect if global ID has been added to owned_
588  for (size_t b = 0; b < blockOrder_.size(); ++b) {
589 
590  if(fa_fps_[b]==Teuchos::null)
591  continue;
592 
593  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
594 
595  for (size_t l = 0; l < myElements.size(); ++l) {
596  const std::vector<GO> & localOrdering = elementGIDs_[myElements[l]];
597 
598  // add "novel" global ids that are also owned to owned array
599  for(std::size_t i=0;i<localOrdering.size();i++) {
600  // don't try to add if ID is not owned
601  if(isOwned.find(localOrdering[i])==isOwned.end())
602  continue;
603 
604  // only add a global ID if it has not been added
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]);
609  }
610  }
611  }
612  }
613 
614  // add any other owned GIDs that were not already included.
615  // these are ones that are owned locally but not required by any
616  // element owned by this processor (uggh!)
617  for(typename HashTable::const_iterator itr=remainingOwned.begin();itr!=remainingOwned.end();itr++)
618  owned_.push_back(*itr);
619 
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.");
624  }
625 
626  }
627 
628  // Build the ghosted_ array. The old simple way led to slow Jacobian
629  // assembly; the new way speeds up Jacobian assembly.
630  {
631  // Loop over all the elements and do a greedy ordering of local values over
632  // the elements for building ghosted_. Hopefully this gives a better
633  // layout for an element-ordered assembly.
634  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::" \
635  "buildGlobalUnknowns::build_ghosted_array");
636 
637  // Use a hash table to detect if global IDs have been added to owned_.
638  typedef std::unordered_set<GO> HashTable;
639  HashTable hashTable;
640  for (std::size_t i = 0; i < owned_.size(); i++)
641  hashTable.insert(owned_[i]);
642 
643  // Here we construct an accessor vector, such that we first process
644  // everything in the current element block, optionally followed by
645  // everything in the neighbor element block.
646  std::vector<ElementBlockAccess> blockAccessVec;
647  blockAccessVec.push_back(ElementBlockAccess(true,connMngr_));
648  if(useNeighbors_)
649  blockAccessVec.push_back(ElementBlockAccess(false,connMngr_));
650  for (std::size_t a = 0; a < blockAccessVec.size(); ++a)
651  {
652  // Get the access type (owned or neighbor).
653  const ElementBlockAccess& access = blockAccessVec[a];
654  for (size_t b = 0; b < blockOrder_.size(); ++b)
655  {
656  if (fa_fps_[b] == Teuchos::null)
657  continue;
658  const std::vector<LO>& myElements =
659  access.getElementBlock(blockOrder_[b]);
660  for (size_t l = 0; l < myElements.size(); ++l)
661  {
662  const std::vector<GO>& localOrdering = elementGIDs_[myElements[l]];
663 
664  // Add "novel" global IDs into the ghosted_ vector.
665  for (std::size_t i = 0; i < localOrdering.size(); ++i)
666  {
667  std::pair<typename HashTable::iterator, bool> insertResult =
668  hashTable.insert(localOrdering[i]);
669 
670  // If the insertion succeeds, then this is "novel" to the owned_
671  // and ghosted_ vectors, so include it in ghosted_.
672  if(insertResult.second)
673  ghosted_.push_back(localOrdering[i]);
674  }
675  }
676  }
677  }
678  }
679 
680  buildConnectivityRun_ = true;
681 
682  // build orientations if required
683  if(requireOrientations_) {
684  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_orientation");
685  buildUnknownsOrientation();
686  }
687 
688  // allocate the local IDs
689  if (useNeighbors_) {
690  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_local_ids_from_owned_and_ghosted");
691  this->buildLocalIdsFromOwnedAndGhostedElements();
692  }
693  else {
694  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_local_ids");
695  this->buildLocalIds();
696  }
697 }
698 
699 template <typename LO, typename GO>
700 std::pair<Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> >,
702 DOFManager<LO,GO>::buildGlobalUnknowns_GUN(const Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & tagged_overlap_mv,
703  Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & overlap_mv) const
704 {
705  // some typedefs
706  typedef panzer::TpetraNodeType Node;
707  typedef Tpetra::Map<LO, GO, Node> Map;
708 
709  typedef Tpetra::Export<LO,GO,Node> Export;
710  typedef Tpetra::Import<LO,GO,Node> Import;
711 
712  //the GIDs are of type GO.
713  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
714 
715  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN");
716 
717  // LINE 2: In the GUN paper
718  RCP<const Map> overlap_map = tagged_overlap_mv.getMap();
719 
720  /* 6. Create a OneToOne map from the overlap map.
721  */
722 
723  // LINE 4: In the GUN paper
724 
725  RCP<const Map> non_overlap_map;
726  if(!useTieBreak_) {
727  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_04 createOneToOne");
728 
729  GreedyTieBreak<LO,GO> greedy_tie_break;
730  non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map, greedy_tie_break);
731  }
732  else {
733  // use a hash tie break to get better load balancing from create one to one
734  // Aug. 4, 2016...this is a bad idea and doesn't work
735  HashTieBreak<LO,GO> tie_break;
736  non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map,tie_break);
737  }
738 
739  /* 7. Create a non-overlapped multivector from OneToOne map.
740  */
741 
742  // LINE 5: In the GUN paper
743 
744  Teuchos::RCP<MultiVector> tagged_non_overlap_mv;
745  {
746  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_05 alloc_unique_mv");
747 
748  tagged_non_overlap_mv = Tpetra::createMultiVector<GO>(non_overlap_map,(size_t)numFields_);
749  }
750 
751  /* 8. Create an export between the two maps.
752  */
753 
754  // LINE 6: In the GUN paper
755  RCP<Export> exp;
756  RCP<Import> imp;
757  RCP<MultiVector> non_overlap_mv;
758  {
759  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_06 export");
760 
761  exp = rcp(new Export(overlap_map,non_overlap_map));
762 
763 #ifdef PANZER_DOFMGR_REQUIRE_CUDA
764  // Note: ETP 04/26/16 Temporarily create an importer for all of the
765  // doImport() calls below. This works around mysterious failures when
766  // using the exporter for Cuda builds.
767  imp = rcp(new Import(non_overlap_map,overlap_map));
768 #endif
769 
770  /* 9. Export data using ABSMAX.
771  */
772  tagged_non_overlap_mv->doExport(tagged_overlap_mv,*exp,Tpetra::ABSMAX);
773 
774  // copy the tagged one, so as to preserve the tagged MV so we can overwrite
775  // the non_overlap_mv
776  non_overlap_mv = rcp(new MultiVector(*tagged_non_overlap_mv,Teuchos::Copy));
777  }
778 
779 
780  /* 10. Compute the local sum using Kokkos.
781  */
782 
783  // LINES 7-9: In the GUN paper
784 
785  GO localsum=0;
786  {
787  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_07-09 local_count");
788 
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);
794  Kokkos::parallel_reduce(mv_size,panzer::dof_functors::SumRank2<GO,KV>(values),localsum);
795  }
796 
797  /* 11. Create a map using local sums to generate final GIDs.
798  */
799 
800  // LINE 10: In the GUN paper
801 
802  GO myOffset = -1;
803  {
804  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_10 prefix_sum");
805 
806  // do a prefix sum
807  GO scanResult = 0;
808  Teuchos::scan<int, GO> (*getComm(), Teuchos::REDUCE_SUM, static_cast<size_t> (localsum), Teuchos::outArg (scanResult));
809  myOffset = scanResult - localsum;
810  }
811 
812  // LINE 11 and 12: In the GUN paper, these steps are eliminated because
813  // the non_overlap_mv is reused
814 
815  /* 12. Iterate through the non-overlapping MV and assign GIDs to
816  * the necessary points. (Assign a -1 elsewhere.)
817  */
818 
819  // LINES 13-21: In the GUN paper
820 
821  {
822  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_13-21 gid_assignment");
823 
824  // ArrayView<const GO> owned_ids = gid_map->getNodeElementList();
825  int which_id=0;
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){
830  // editnonoverlap[j][i]=myOffset+which_id;
831  int ndof = Teuchos::as<int>(editnonoverlap[j][i]);
832  editnonoverlap[j][i]=myOffset+which_id;
833  which_id+=ndof;
834  }
835  else{
836  editnonoverlap[j][i]=-1;
837  }
838 
839  }
840  }
841  }
842 
843  // LINE 22: In the GUN paper. Were performed above, and the overlaped_mv is
844  // abused to handle input tagging.
845 
846  /* 13. Import data back to the overlap MV using REPLACE.
847  */
848 
849  // LINE 23: In the GUN paper
850 
851  {
852  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_23 final_import");
853 
854 #ifdef PANZER_DOFMGR_REQUIRE_CUDA
855  overlap_mv.doImport(*non_overlap_mv,*imp,Tpetra::REPLACE);
856 #else
857  // use exporter to save on communication setup costs
858  overlap_mv.doImport(*non_overlap_mv,*exp,Tpetra::REPLACE);
859 #endif
860  }
861 
862  //std::cout << Teuchos::describe(*non_overlap_mv,Teuchos::VERB_EXTREME) << std::endl;
863 
864  // return non_overlap_mv;
865  return std::make_pair(non_overlap_mv,tagged_non_overlap_mv);
866 }
867 
868 template <typename LO, typename GO>
871 {
872  // some typedefs
873  typedef panzer::TpetraNodeType Node;
874  typedef Tpetra::Map<LO, GO, Node> Map;
875  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
876 
877  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector");
878 
879  //We will iterate through all of the blocks, building a FieldAggPattern for
880  //each of them.
881 
882  for (size_t b = 0; b < blockOrder_.size(); ++b) {
883  std::vector<std::tuple< int, panzer::FieldType, RCP<const panzer::FieldPattern> > > faConstruct;
884  //The ID is going to be the AID, and then everything will work.
885  //The ID should not be the AID, it should be the ID it has in the ordering.
886 
887  for (size_t i = 0; i < fieldAIDOrder_.size(); ++i) {
888  int looking = fieldAIDOrder_[i];
889 
890  //Check if in b's fp list
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]]));
894  }
895 
896  }
897 
898  if(faConstruct.size()>0) {
899  fa_fps_.push_back(rcp(new FieldAggPattern(faConstruct, ga_fp_)));
900 
901  // how many global IDs are in this element block?
902  int gidsInBlock = fa_fps_[fa_fps_.size()-1]->numberIds();
903  elementBlockGIDCount_.push_back(gidsInBlock);
904  }
905  else {
906  fa_fps_.push_back(Teuchos::null);
907  elementBlockGIDCount_.push_back(0);
908  }
909  }
910 
911  RCP<const Map> overlapmap = buildOverlapMapFromElements(ownedAccess);
912 
913  // LINE 22: In the GUN paper...the overlap_mv is reused for the tagged multivector.
914  // This is a bit of a practical abuse of the algorithm presented in the paper.
915 
916  Teuchos::RCP<MultiVector> overlap_mv;
917  {
918  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector::allocate_tagged_multivector");
919 
920  overlap_mv = Tpetra::createMultiVector<GO>(overlapmap,(size_t)numFields_);
921  overlap_mv->putScalar(0); // if tpetra is not initialized with zeros
922  }
923 
924  /* 5. Iterate through all local elements again, checking with the FP
925  * information. Mark up the overlap map accordingly.
926  */
927 
928  {
929  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector::fill_tagged_multivector");
930 
931  // temporary working vector to fill each row in tagged array
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) {
935  // there has to be a field pattern assocaited with the block
936  if(fa_fps_[b]==Teuchos::null)
937  continue;
938 
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]);
945  int offset=0;
946  for (int c = 0; c < connSize; ++c) {
947  size_t lid = overlapmap->getLocalElement(elmtConn[c]);
948 
949  for(std::size_t i=0;i<working.size();i++)
950  working[i] = 0;
951  for (int n = 0; n < numFields[c]; ++n) {
952  int whichField = fieldIds[offset];
953  //Row will be lid. column will be whichField.
954  //Shove onto local ordering
955  working[whichField]++;
956  offset++;
957  }
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];
961  }
962 
963  }
964  }
965  }
966 
967  // // verbose output for inspecting overlap_mv
968  // for(int i=0;i<overlap_mv->getLocalLength(); i++) {
969  // for(int j=0;j<overlap_mv->getNumVectors() ; j++)
970  // std::cout << edittwoview[j][i] << " ";
971  // std::cout << std::endl;
972  // }
973  }
974 
975  return overlap_mv;
976 }
977 
978 template <typename LO, typename GO>
979 int DOFManager<LO,GO>::getFieldNum(const std::string & string) const
980 {
981  int ind=0;
982  bool found=false;
983  while(!found && (size_t)ind<fieldStringOrder_.size()){
984  if(fieldStringOrder_[ind]==string)
985  found=true;
986  else
987  ind++;
988  }
989 
990  if(found)
991  return ind;
992 
993  // didn't find anything...return -1
994  return -1;
995 }
996 
997 template <typename LO, typename GO>
998 void DOFManager<LO,GO>::getFieldOrder(std::vector<std::string> & fieldOrder) const
999 {
1000  fieldOrder.resize(fieldStringOrder_.size());
1001  for (size_t i = 0; i < fieldStringOrder_.size(); ++i)
1002  fieldOrder[i]=fieldStringOrder_[i];
1003 }
1004 
1005 template <typename LO, typename GO>
1006 bool DOFManager<LO,GO>::fieldInBlock(const std::string & field, const std::string & block) const
1007 {
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());
1013  }
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());
1019  }
1020  int fid=fitr->second;
1021  int bid=bitr->second;
1022 
1023  bool found=false;
1024  for (size_t i = 0; i < blockToAssociatedFP_[bid].size(); ++i) {
1025  if(blockToAssociatedFP_[bid][i]==fid){
1026  found=true;
1027  break;
1028  }
1029  }
1030 
1031  return found;
1032 }
1033 
1034 template <typename LO, typename GO>
1035 const std::vector<int> & DOFManager<LO,GO>::getBlockFieldNumbers(const std::string & blockId) const
1036 {
1037  TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,"DOFManager::getBlockFieldNumbers: BuildConnectivity must be run first.");
1038 
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;
1043 
1044  // there has to be a field pattern assocaited with the block
1045  if(fa_fps_[bid]!=Teuchos::null)
1046  return fa_fps_[bid]->fieldIds();
1047 
1048  // nothing to return
1049  static std::vector<int> empty;
1050  return empty;
1051 }
1052 
1053 template <typename LO, typename GO>
1054 const std::pair<std::vector<int>,std::vector<int> > &
1055 DOFManager<LO,GO>::getGIDFieldOffsets_closure(const std::string & blockId, int fieldNum, int subcellDim,int subcellId) const
1056 {
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.");
1061 
1062  // there has to be a field pattern assocaited with the block
1063  if(fa_fps_[bitr->second]!=Teuchos::null)
1064  return fa_fps_[bitr->second]->localOffsets_closure(fieldNum, subcellDim, subcellId);
1065 
1066  static std::pair<std::vector<int>,std::vector<int> > empty;
1067  return empty;
1068 }
1069 
1070 template <typename LO, typename GO>
1071 void DOFManager<LO,GO>::ownedIndices(const std::vector<GO> & indices,std::vector<bool> & isOwned) const
1072 {
1073  //Resizes the isOwned array.
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 );
1079  }
1080 }
1081 
1082 template <typename LO, typename GO>
1083 void DOFManager<LO,GO>::setFieldOrder(const std::vector<std::string> & fieldOrder)
1084 {
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;
1090  //We also need to reassign the fieldAIDOrder_.
1091  for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1092  fieldAIDOrder_[i]=fieldNameToAID_[fieldStringOrder_[i]];
1093  }
1094  }
1095  else
1096  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::setFieldOrder: Invalid Field Ordering!");
1097 
1098 }
1099 
1100 
1101 template <typename LO, typename GO>
1102 bool DOFManager<LO,GO>::validFieldOrder(const std::vector<std::string> & proposed_fieldOrder)
1103 {
1104  if(fieldStringOrder_.size()!=proposed_fieldOrder.size())
1105  return false;
1106  //I'm using a not very efficient way of doing this, but there should never
1107  //be that many fields, so it isn't that big of a deal.
1108  for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
1109  bool found=false;
1110  for (size_t j = 0; j < proposed_fieldOrder.size(); ++j) {
1111  if(fieldStringOrder_[i]==proposed_fieldOrder[j])
1112  found=true;
1113  }
1114  if(!found)
1115  return false;
1116  }
1117  return true;
1118 }
1119 
1120 template <typename LO, typename GO>
1121 const std::string & DOFManager<LO,GO>::getFieldString(int num) const
1122 {
1123  //A reverse lookup through fieldStringOrder_.
1124  if(num>=(int)fieldStringOrder_.size())
1125  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "DOFManager::getFieldString: invalid number");
1126  return fieldStringOrder_[num];
1127 }
1128 
1129 //Everything associated with orientation is not yeat built, but this
1130 //is the method as found in Panzer_DOFManager_impl.hpp. There are
1131 //going to need to be some substantial changes to the code as it applies
1132 //to this DOFManager, but the basic ideas and format should be similar.
1133 //
1134 template <typename LO, typename GO>
1136 {
1137  orientation_.clear(); // clean up previous work
1138 
1139  std::vector<std::string> elementBlockIds;
1140  connMngr_->getElementBlockIds(elementBlockIds);
1141 
1142  // figure out how many total elements are owned by this processor (why is this so hard!)
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();
1146 
1147  // allocate for each block
1148  orientation_.resize(myElementCount);
1149 
1150  // loop over all element blocks
1151  for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin();
1152  blockItr!=elementBlockIds.end();++blockItr) {
1153  const std::string & blockName = *blockItr;
1154 
1155  // this block has no unknowns (or elements)
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");
1159  }
1160 
1161  int bid=fap->second;
1162 
1163  if(fa_fps_[bid]==Teuchos::null)
1164  continue;
1165 
1166  // grab field patterns, will be necessary to compute orientations
1167  const FieldPattern & fieldPattern = *fa_fps_[bid];
1168 
1169  //Should be ga_fp_ (geometric aggregate field pattern)
1170  std::vector<std::pair<int,int> > topEdgeIndices;
1171  orientation_helpers::computePatternEdgeIndices(*ga_fp_,topEdgeIndices);
1172 
1173  // grab face orientations if 3D
1174  std::vector<std::vector<int> > topFaceIndices;
1175  if(ga_fp_->getDimension()==3)
1176  orientation_helpers::computePatternFaceIndices(*ga_fp_,topFaceIndices);
1177 
1178  //How many GIDs are associated with a particular element bloc
1179  const std::vector<LO> & elmts = connMngr_->getElementBlock(blockName);
1180  for(std::size_t e=0;e<elmts.size();e++) {
1181  // this is the vector of orientations to fill: initialize it correctly
1182  std::vector<signed char> & eOrientation = orientation_[elmts[e]];
1183 
1184  // This resize seems to be the same as fieldPattern.numberIDs().
1185  // When computer edge orientations is called, that is the assert.
1186  // There should be no reason to make it anymore complicated.
1187  eOrientation.resize(fieldPattern.numberIds());
1188  for(std::size_t s=0;s<eOrientation.size();s++)
1189  eOrientation[s] = 1; // put in 1 by default
1190 
1191  // get geometry ids
1192  LO connSz = connMngr_->getConnectivitySize(elmts[e]);
1193  const GO * connPtr = connMngr_->getConnectivity(elmts[e]);
1194  const std::vector<GO> connectivity(connPtr,connPtr+connSz);
1195 
1196  orientation_helpers::computeCellEdgeOrientations(topEdgeIndices, connectivity, fieldPattern, eOrientation);
1197 
1198  // compute face orientations in 3D
1199  if(ga_fp_->getDimension()==3)
1200  orientation_helpers::computeCellFaceOrientations(topFaceIndices, connectivity, fieldPattern, eOrientation);
1201  }
1202  }
1203 }
1204 
1205 template <typename LO, typename GO>
1206 void DOFManager<LO,GO>::getElementOrientation(LO localElmtId,std::vector<double> & gidsOrientation) const
1207 {
1208  TEUCHOS_TEST_FOR_EXCEPTION(orientation_.size()==0,std::logic_error,
1209  "DOFManager::getElementOrientations: Orientations were not constructed!");
1210 
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]);
1215  }
1216 }
1217 
1218 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1220 {
1222 
1223  connMngr_ = Teuchos::null;
1224 
1225  // wipe out FEI objects
1226  orientation_.clear(); // clean up previous work
1227  fa_fps_.clear();
1228  elementGIDs_.clear();
1229  owned_.clear();
1230  ghosted_.clear();
1231  elementBlockGIDCount_.clear();
1232 
1233  return connMngr;
1234 }
1235 
1236 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1237 std::size_t DOFManager<LocalOrdinalT,GlobalOrdinalT>::blockIdToIndex(const std::string & blockId) const
1238 {
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;
1243 }
1244 
1245 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1247 {
1248  os << "DOFManager Field Information: " << std::endl;
1249 
1250  // sanity check
1251  TEUCHOS_ASSERT(blockOrder_.size()==blockToAssociatedFP_.size());
1252 
1253  for(std::size_t i=0;i<blockOrder_.size();i++) {
1254  os << " Element Block = " << blockOrder_[i] << std::endl;
1255 
1256  // output field information
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;
1260  }
1261 }
1262 
1263 template <typename LO,typename GO>
1267 {
1268  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::builderOverlapMapFromElements");
1269 
1270  /*
1271  * 2. Iterate through all local elements and create the overlapVector
1272  * of concerned elements.
1273  */
1274 
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]);
1283  }
1284  }
1285  }
1286 
1287  Array<GO> overlapVector;
1288  for (typename std::set<GO>::const_iterator itr = overlapset.begin(); itr!=overlapset.end(); ++itr) {
1289  overlapVector.push_back(*itr);
1290  }
1291 
1292  /* 3. Construct an overlap map from this structure.
1293  */
1294  return Tpetra::createNonContigMap<LO,GO>(overlapVector,getComm());
1295 }
1296 
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
1303 {
1304  using Teuchos::ArrayRCP;
1305 
1306  //To generate elementGIDs we need to go through all of the local elements.
1307  ArrayRCP<ArrayRCP<const GO> > twoview = overlap_mv.get2dView();
1308 
1309  //And for each of the things in fa_fp.fieldIds we go to that column. To the the row,
1310  //we move from globalID to localID in the map and use our local value for something.
1311  for (size_t b = 0; b < blockOrder_.size(); ++b) {
1312  const std::vector<LO> & myElements = access.getElementBlock(blockOrder_[b]);
1313 
1314  if(fa_fps_[b]==Teuchos::null) {
1315  // fill elements that are not used with empty vectors
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);
1320  }
1321  continue;
1322  }
1323 
1324  const std::vector<int> & numFields= fa_fps_[b]->numFieldsPerId();
1325  const std::vector<int> & fieldIds= fa_fps_[b]->fieldIds();
1326  //
1327  //
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;
1332  int offset=0;
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];
1338  offset++;
1339  //Row will be lid. column will be whichField.
1340  //Shove onto local ordering
1341  localOrdering.push_back(twoview[whichField][lid]+dofsPerField[whichField]);
1342 
1343  dofsPerField[whichField]++;
1344  }
1345  }
1346  LO thisID=myElements[l];
1347  if(elementGIDs.size()<=(size_t)thisID){
1348  elementGIDs.resize(thisID+1);
1349  }
1350  elementGIDs[thisID]=localOrdering;
1351  }
1352  }
1353 }
1354 
1355 template <typename LO, typename GO>
1357 {
1358  std::vector<std::vector<LO> > elementLIDs(elementGIDs_.size());
1359 
1360  std::vector<GO> ownedAndGhosted;
1361  this->getOwnedAndGhostedIndices(ownedAndGhosted);
1362 
1363  // build global to local hash map (temporary and used only once)
1364  std::unordered_map<GO,LO> hashMap;
1365  for(std::size_t i = 0; i < ownedAndGhosted.size(); ++i)
1366  hashMap[ownedAndGhosted[i]] = i;
1367 
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]];
1374  }
1375 
1376  this->setLocalIds(elementLIDs);
1377 }
1378 
1379 /*
1380 template <typename LO,typename GO>
1381 Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> >
1382 DOFManager<LO,GO>::runLocalRCMReordering(const Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> > & map)
1383 {
1384  typedef panzer::TpetraNodeType Node;
1385  typedef Tpetra::Map<LO, GO, Node> Map;
1386  typedef Tpetra::CrsGraph<LO, GO, Node> Graph;
1387 
1388  Teuchos::RCP<Graph> graph = Teuchos::rcp(new Graph(map,0));
1389 
1390  // build Crs Graph from the mesh
1391  for (size_t b = 0; b < blockOrder_.size(); ++b) {
1392  if(fa_fps_[b]==Teuchos::null)
1393  continue;
1394 
1395  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
1396  for (size_t l = 0; l < myElements.size(); ++l) {
1397  LO connSize = connMngr_->getConnectivitySize(myElements[l]);
1398  const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
1399  for (int c = 0; c < connSize; ++c) {
1400  LO lid = map->getLocalElement(elmtConn[c]);
1401  if(Teuchos::OrdinalTraits<LO>::invalid()!=lid)
1402  continue;
1403 
1404  graph->insertGlobalIndices(elmtConn[c],Teuchos::arrayView(elmtConn,connSize));
1405  }
1406  }
1407  }
1408 
1409  graph->fillComplete();
1410 
1411  std::vector<GO> newOrder(map->getNodeNumElements());
1412  {
1413  // graph is constructed, now run RCM using zoltan2
1414  typedef Zoltan2::XpetraCrsGraphInput<Graph> SparseGraphAdapter;
1415 
1416  Teuchos::ParameterList params;
1417  params.set("order_method", "rcm");
1418  SparseGraphAdapter adapter(graph);
1419 
1420  Zoltan2::OrderingProblem<SparseGraphAdapter> problem(&adapter, &params,MPI_COMM_SELF);
1421  problem.solve();
1422 
1423  // build a new global ording array using permutation
1424  Zoltan2::OrderingSolution<GO,LO> * soln = problem.getSolution();
1425 
1426  size_t dummy;
1427  size_t checkLength = soln->getPermutationSize();
1428  LO * checkPerm = soln->getPermutation(&dummy);
1429 
1430  Teuchos::ArrayView<const GO > oldOrder = map->getNodeElementList();
1431  TEUCHOS_ASSERT(checkLength==oldOrder.size());
1432  TEUCHOS_ASSERT(checkLength==newOrder.size());
1433 
1434  for(size_t i=0;i<checkLength;i++)
1435  newOrder[checkPerm[i]] = oldOrder[i];
1436  }
1437 
1438  return Tpetra::createNonContigMap<LO,GO>(newOrder,communicator_);
1439 }
1440 */
1441 
1442 } /*panzer*/
1443 
1444 #endif
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.
size_type size() const
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
basic_FancyOStream & setShowProcRank(const bool showProcRank)
#define PANZER_DOFMGR_FUNC_TIME_MONITOR(a)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
int numFields
const unsigned int seed_
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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 push_back(const value_type &x)
void printFieldInformation(std::ostream &os) const
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;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
virtual const std::vector< int > & getSubcellIndices(int dim, int cellIndex) const =0
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...
#define TEUCHOS_ASSERT(assertion_test)
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.