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 
48 #include "mpi.h"
49 
50 #include "PanzerDofMgr_config.hpp"
51 #include "Panzer_FieldPattern.hpp"
54 #include "Panzer_ConnManager.hpp"
58 #include "Teuchos_GlobalMPISession.hpp"
60 
61 #include "Teuchos_RCP.hpp"
62 #include "Teuchos_Array.hpp"
63 #include "Teuchos_ArrayView.hpp"
64 
65 #include "Tpetra_DefaultPlatform.hpp"
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 using Teuchos::RCP;
112 using Teuchos::rcp;
113 using Teuchos::ArrayRCP;
114 using Teuchos::Array;
115 using Teuchos::ArrayView;
116 
117 template <typename LO, typename GO>
119  : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
120 { }
121 
122 template <typename LO, typename GO>
123 DOFManager<LO,GO>::DOFManager(const Teuchos::RCP<ConnManager<LO,GO> > & connMngr,MPI_Comm mpiComm)
124  : numFields_(0),buildConnectivityRun_(false),requireOrientations_(false), useTieBreak_(false), useNeighbors_(false)
125 {
126  setConnManager(connMngr,mpiComm);
127 }
128 
129 template <typename LO, typename GO>
130 void DOFManager<LO,GO>::setConnManager(const Teuchos::RCP<ConnManager<LO,GO> > & connMngr, MPI_Comm mpiComm)
131 {
132  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
133  "DOFManager::setConnManager: setConnManager cannot be called after "
134  "buildGlobalUnknowns has been called");
135  connMngr_ = connMngr;
136  //We acquire the block ordering from the connmanager.
137  connMngr_->getElementBlockIds(blockOrder_);
138  for (size_t i = 0; i < blockOrder_.size(); ++i) {
139  blockNameToID_.insert(std::map<std::string,int>::value_type(blockOrder_[i],i));
140  //We must also initialize vectors for FP associations.
141  }
142  blockToAssociatedFP_.resize(blockOrder_.size());
143  communicator_ = Teuchos::rcp(new Teuchos::MpiComm<int>(Teuchos::opaqueWrapper(mpiComm)));
144 }
145 
146 
147 //Adds a field to be used in creating the Global Numbering
148 //Returns the index for the field pattern
149 template <typename LO, typename GO>
150 int DOFManager<LO,GO>::addField(const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern)
151 {
152  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
153  "DOFManager::addField: addField cannot be called after "
154  "buildGlobalUnknowns has been called");
155 
156  fieldPatterns_.push_back(pattern);
157  fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
158 
159  //The default values for IDs are the sequential order they are added in.
160  fieldStringOrder_.push_back(str);
161  fieldAIDOrder_.push_back(numFields_);
162 
163  for(std::size_t i=0;i<blockOrder_.size();i++) {
164  blockToAssociatedFP_[i].push_back(numFields_);
165  }
166 
167  ++numFields_;
168  return numFields_-1;
169 }
170 
171 template <typename LO, typename GO>
172 int DOFManager<LO,GO>::addField(const std::string & blockID, const std::string & str, const Teuchos::RCP<const FieldPattern> & pattern)
173 {
174  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
175  "DOFManager::addField: addField cannot be called after "
176  "buildGlobalUnknowns has been called");
177  TEUCHOS_TEST_FOR_EXCEPTION((connMngr_==Teuchos::null),std::logic_error,
178  "DOFManager::addField: you must add a ConnManager before"
179  "you can associate a FP with a given block.")
180  bool found=false;
181  size_t blocknum=0;
182  while(!found && blocknum<blockOrder_.size()){
183  if(blockOrder_[blocknum]==blockID){
184  found=true;
185  break;
186  }
187  blocknum++;
188  }
189  TEUCHOS_TEST_FOR_EXCEPTION(!found,std::logic_error, "DOFManager::addField: Invalid block name.");
190 
191  //This will be different if the FieldPattern is already present.
192  //We need to check for that.
193  found=false;
194  std::map<std::string,int>::const_iterator fpIter = fieldNameToAID_.find(str);
195  if(fpIter!=fieldNameToAID_.end())
196  found=true;
197 
198  if(!found){
199  fieldPatterns_.push_back(pattern);
200  fieldNameToAID_.insert(std::map<std::string,int>::value_type(str, numFields_));
201  //The default values for IDs are the sequential order they are added in.
202  fieldStringOrder_.push_back(str);
203  fieldAIDOrder_.push_back(numFields_);
204 
205  //This is going to be associated with blocknum.
206  blockToAssociatedFP_[blocknum].push_back(numFields_);
207  ++numFields_;
208  return numFields_-1;
209  }
210  else{
211  blockToAssociatedFP_[blocknum].push_back(fpIter->second);
212  return numFields_;
213  }
214 }
215 
216 
217 
218  //Returns the fieldpattern of the given name
219 template <typename LO, typename GO>
220 Teuchos::RCP<const FieldPattern> DOFManager<LO,GO>::getFieldPattern(const std::string & name) const
221 {
222  std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(name);
223  if(fitr==fieldNameToAID_.end())
224  return Teuchos::null;
225 
226  if(fitr->second<int(fieldPatterns_.size()))
227  return fieldPatterns_[fitr->second];
228 
229  return Teuchos::null;
230 }
231 
232 //Returns the fieldpattern of the given name
233 template <typename LO, typename GO>
234 Teuchos::RCP<const FieldPattern> DOFManager<LO,GO>::getFieldPattern(const std::string & blockId, const std::string & fieldName) const
235 {
236  std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(fieldName);
237  if(fitr==fieldNameToAID_.end())
238  return Teuchos::null;
239 
240  bool found=false;
241  size_t blocknum=0;
242  while(!found && blocknum<blockOrder_.size()){
243  if(blockOrder_[blocknum]==blockId){
244  found=true;
245  break;
246  }
247  blocknum++;
248  }
249 
250  if(!found)
251  return Teuchos::null;
252 
253  std::vector<int>::const_iterator itr = std::find(blockToAssociatedFP_[blocknum].begin(),
254  blockToAssociatedFP_[blocknum].end(),fitr->second);
255  if(itr!=blockToAssociatedFP_[blocknum].end()) {
256  if(fitr->second<int(fieldPatterns_.size()))
257  return fieldPatterns_[fitr->second];
258  }
259 
260  return Teuchos::null;
261 }
262 
263 template <typename LO, typename GO>
264 void DOFManager<LO, GO>::getOwnedIndices(std::vector<GO>& indices) const
265 {
266  indices = owned_;
267 }
268 
269 template <typename LO, typename GO>
270 void DOFManager<LO, GO>::getOwnedAndGhostedIndices(std::vector<GO>& indices)
271  const
272 {
273  indices.resize(owned_.size() + ghosted_.size());
274  for (size_t i = 0; i < owned_.size(); ++i)
275  indices[i] = owned_[i];
276  for (size_t i = 0; i < ghosted_.size(); ++i)
277  indices[owned_.size() + i] = ghosted_[i];
278 }
279 
280  //gets the number of fields
281 template <typename LO, typename GO>
283 {
284  return numFields_;
285 }
286 
287 template <typename LO, typename GO>
288 const std::vector<int> & DOFManager<LO,GO>::getGIDFieldOffsets(const std::string & blockID, int fieldNum) const
289 {
290  TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error, "DOFManager::getGIDFieldOffsets: cannot be called before "
291  "buildGlobalUnknowns has been called");
292  std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockID);
293  if(bitr==blockNameToID_.end())
294  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
295  int bid=bitr->second;
296  if(fa_fps_[bid]!=Teuchos::null)
297  return fa_fps_[bid]->localOffsets(fieldNum);
298 
299  static const std::vector<int> empty;
300  return empty;
301 }
302 
303 template <typename LO, typename GO>
304 void DOFManager<LO,GO>::getElementGIDs(LO localElementID, std::vector<GO> & gids, const std::string & blockIdHint) const
305 {
306  gids = elementGIDs_[localElementID];
307 }
308 
309 template <typename LocalOrdinalT,typename GlobalOrdinalT>
311 {
312  /* STEPS.
313  * 1. Build GA_FP and all block's FA_FP's and place into respective data structures.
314  */
315  if(requireOrientations_){
316  fieldPatterns_.push_back(Teuchos::rcp(new NodalFieldPattern(fieldPatterns_[0]->getCellTopology())));
317  }
318  RCP<GeometricAggFieldPattern> aggFieldPattern = Teuchos::rcp(new GeometricAggFieldPattern);;
319  aggFieldPattern = Teuchos::rcp(new GeometricAggFieldPattern(fieldPatterns_));
320 
321  connMngr_->buildConnectivity(*aggFieldPattern);
322 
323  // using new geometric pattern, build global unknowns
324  buildGlobalUnknowns(aggFieldPattern);
325 }
326 
327 //builds the global unknowns array
328 template <typename LO, typename GO>
329 void DOFManager<LO,GO>::buildGlobalUnknowns(const Teuchos::RCP<const FieldPattern> & geomPattern)
330 {
331  // some typedefs
332  typedef panzer::TpetraNodeType Node;
333  typedef Tpetra::Map<LO, GO, Node> Map;
334 
335  typedef Tpetra::Import<LO,GO,Node> Import;
336 
337  //the GIDs are of type GO.
338  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
339 
340  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns");
341 
342  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
343  out.setOutputToRootOnly(-1);
344  out.setShowProcRank(true);
345 
346  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
347  "DOFManager::buildGlobalUnknowns: buildGlobalUnknowns cannot be called again "
348  "after buildGlobalUnknowns has been called");
349 
350  // this is a safety check to make sure that nodes are included
351  // in the geometric field pattern when orientations are required
352  if(getOrientationsRequired()) {
353  std::size_t sz = geomPattern->getSubcellIndices(0,0).size();
354 
355  TEUCHOS_TEST_FOR_EXCEPTION(sz==0,std::logic_error,
356  "DOFManager::buildGlobalUnknowns requires a geometric pattern including "
357  "the nodes when orientations are needed!");
358  }
359 
360  /* STEPS.
361  * 1. Build all block's FA_FP's and place into respective data structures.
362  */
363  ga_fp_ = geomPattern;
364 
365  // given a set of elements over each element block build an overlap
366  // map that will provide the required element entities for the
367  // set of elements requested.
368  ElementBlockAccess ownedAccess(true,connMngr_);
369 
370  // INPUT: To the algorithm in the GUN paper
371  RCP<MultiVector> tagged_overlap_mv = buildTaggedMultiVector(ownedAccess);
372  RCP<const Map> overlap_map = tagged_overlap_mv->getMap();
373 
374  RCP<MultiVector> overlap_mv = Tpetra::createMultiVector<GO>(overlap_map,(size_t)numFields_);
375 
376  // call the GUN paper algorithm
377  auto non_overlap_pair = buildGlobalUnknowns_GUN(*tagged_overlap_mv,*overlap_mv);
378  RCP<MultiVector> non_overlap_mv = non_overlap_pair.first;
379  RCP<MultiVector> tagged_non_overlap_mv = non_overlap_pair.second;
380  RCP<const Map> non_overlap_map = non_overlap_mv->getMap();
381 
382  /* 14. Cross reference local element connectivity and overlap map to
383  * create final GID vectors.
384  */
385 
386  // this bit of code takes the uniquely assigned GIDs and spreads them
387  // out for processing by local element ID
388  fillGIDsFromOverlappedMV(ownedAccess,elementGIDs_,*overlap_map,*overlap_mv);
389 
390  // if neighbor unknowns are required, then make sure they are included
391  // in the elementGIDs_
392  if (useNeighbors_) { // enabling this turns on GID construction for
393  // neighbor processors
394  ElementBlockAccess neighborAccess(false,connMngr_);
395  RCP<const Map> overlap_map_neighbor =
396  buildOverlapMapFromElements(neighborAccess);
397 
398  // Export e(overlap_map_neighbor,non_overlap_map);
399  Import imp_neighbor(non_overlap_map,overlap_map_neighbor);
400 
401  Teuchos::RCP<MultiVector> overlap_mv_neighbor =
402  Tpetra::createMultiVector<GO>(overlap_map_neighbor, (size_t)numFields_);
403 
404  // get all neighbor information
405  overlap_mv_neighbor->doImport(*non_overlap_mv, imp_neighbor,
406  Tpetra::REPLACE);
407 
408  fillGIDsFromOverlappedMV(neighborAccess, elementGIDs_,
409  *overlap_map_neighbor, *overlap_mv_neighbor);
410  }
411 
413  // this is where the code is modified to artificially induce GIDs
414  // over 2 Billion unknowns
416  #if 0
417  {
418  panzer::Ordinal64 offset = 0xFFFFFFFFLL;
419 
420  for (size_t b = 0; b < blockOrder_.size(); ++b) {
421  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
422  for (std::size_t l = 0; l < myElements.size(); ++l) {
423  std::vector<GO> & localGIDs = elementGIDs_[myElements[l]];
424  for(std::size_t c=0;c<localGIDs.size();c++)
425  localGIDs[c] += offset;
426  }
427  }
428 
429  Teuchos::ArrayRCP<GO> nvals = non_overlap_mv->get1dViewNonConst();
430  for (int j = 0; j < nvals.size(); ++j)
431  nvals[j] += offset;
432  }
433  #endif
434 
435  // build owned vector
436  {
437  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_owned_vector");
438 
439  typedef std::unordered_set<GO> HashTable;
440  HashTable isOwned, remainingOwned;
441 
442  // owned_ is made up of owned_ids.: This doesn't work for high order
443  Teuchos::ArrayRCP<const GO> nvals = non_overlap_mv->get1dView();
444  Teuchos::ArrayRCP<const GO> tagged_vals = tagged_non_overlap_mv->get1dView();
445  TEUCHOS_ASSERT(nvals.size()==tagged_vals.size());
446  for (int j = 0; j < nvals.size(); ++j) {
447  if (nvals[j] != -1) {
448  for(GO offset=0;offset<tagged_vals[j];offset++)
449  isOwned.insert(nvals[j]+offset);
450  }
451  else {
452  // sanity check
453  TEUCHOS_ASSERT(tagged_vals[j]==0)
454  }
455  }
456  remainingOwned = isOwned;
457 
458  HashTable hashTable; // use to detect if global ID has been added to owned_
459  for (size_t b = 0; b < blockOrder_.size(); ++b) {
460 
461  if(fa_fps_[b]==Teuchos::null)
462  continue;
463 
464  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
465 
466  for (size_t l = 0; l < myElements.size(); ++l) {
467  const std::vector<GO> & localOrdering = elementGIDs_[myElements[l]];
468 
469  // add "novel" global ids that are also owned to owned array
470  for(std::size_t i=0;i<localOrdering.size();i++) {
471  // don't try to add if ID is not owned
472  if(isOwned.find(localOrdering[i])==isOwned.end())
473  continue;
474 
475  // only add a global ID if it has not been added
476  std::pair<typename HashTable::iterator,bool> insertResult = hashTable.insert(localOrdering[i]);
477  if(insertResult.second) {
478  owned_.push_back(localOrdering[i]);
479  remainingOwned.erase(localOrdering[i]);
480  }
481  }
482  }
483  }
484 
485  // add any other owned GIDs that were not already included.
486  // these are ones that are owned locally but not required by any
487  // element owned by this processor (uggh!)
488  for(typename HashTable::const_iterator itr=remainingOwned.begin();itr!=remainingOwned.end();itr++)
489  owned_.push_back(*itr);
490 
491  if(owned_.size()!=isOwned.size()) {
492  out << "I'm about to hang because of unknown numbering failure ... sorry! (line = " << __LINE__ << ")" << std::endl;
493  TEUCHOS_TEST_FOR_EXCEPTION(owned_.size()!=isOwned.size(),std::logic_error,
494  "DOFManager::buildGlobalUnkonwns: Failure because not all owned unknowns have been accounted for.");
495  }
496 
497  }
498 
499  // Build the ghosted_ array. The old simple way led to slow Jacobian
500  // assembly; the new way speeds up Jacobian assembly.
501  {
502  // Loop over all the elements and do a greedy ordering of local values over
503  // the elements for building ghosted_. Hopefully this gives a better
504  // layout for an element-ordered assembly.
505  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::" \
506  "buildGlobalUnknowns::build_ghosted_array");
507 
508  // Use a hash table to detect if global IDs have been added to owned_.
509  typedef std::unordered_set<GO> HashTable;
510  HashTable hashTable;
511  for (std::size_t i = 0; i < owned_.size(); i++)
512  hashTable.insert(owned_[i]);
513 
514  // Here we construct an accessor vector, such that we first process
515  // everything in the current element block, optionally followed by
516  // everything in the neighbor element block.
517  std::vector<ElementBlockAccess> blockAccessVec;
518  blockAccessVec.push_back(ElementBlockAccess(true,connMngr_));
519  if(useNeighbors_)
520  blockAccessVec.push_back(ElementBlockAccess(false,connMngr_));
521  for (std::size_t a = 0; a < blockAccessVec.size(); ++a)
522  {
523  // Get the access type (owned or neighbor).
524  const ElementBlockAccess& access = blockAccessVec[a];
525  for (size_t b = 0; b < blockOrder_.size(); ++b)
526  {
527  if (fa_fps_[b] == Teuchos::null)
528  continue;
529  const std::vector<LO>& myElements =
530  access.getElementBlock(blockOrder_[b]);
531  for (size_t l = 0; l < myElements.size(); ++l)
532  {
533  const std::vector<GO>& localOrdering = elementGIDs_[myElements[l]];
534 
535  // Add "novel" global IDs into the ghosted_ vector.
536  for (std::size_t i = 0; i < localOrdering.size(); ++i)
537  {
538  std::pair<typename HashTable::iterator, bool> insertResult =
539  hashTable.insert(localOrdering[i]);
540 
541  // If the insertion succeeds, then this is "novel" to the owned_
542  // and ghosted_ vectors, so include it in ghosted_.
543  if(insertResult.second)
544  ghosted_.push_back(localOrdering[i]);
545  }
546  }
547  }
548  }
549  }
550 
551  buildConnectivityRun_ = true;
552 
553  // build orientations if required
554  if(requireOrientations_) {
555  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_orientation");
556  buildUnknownsOrientation();
557  }
558 
559  // allocate the local IDs
560  if (useNeighbors_) {
561  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_local_ids_from_owned_and_ghosted");
562  this->buildLocalIdsFromOwnedAndGhostedElements();
563  }
564  else {
565  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns::build_local_ids");
566  this->buildLocalIds();
567  }
568 }
569 
570 template <typename LO, typename GO>
571 std::pair<Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> >,
572  Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> > >
573 DOFManager<LO,GO>::buildGlobalUnknowns_GUN(const Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & tagged_overlap_mv,
574  Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & overlap_mv) const
575 {
576  // some typedefs
577  typedef panzer::TpetraNodeType Node;
578  typedef Tpetra::Map<LO, GO, Node> Map;
579 
580  typedef Tpetra::Export<LO,GO,Node> Export;
581  typedef Tpetra::Import<LO,GO,Node> Import;
582 
583  //the GIDs are of type GO.
584  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
585 
586  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN");
587 
588  // LINE 2: In the GUN paper
589  RCP<const Map> overlap_map = tagged_overlap_mv.getMap();
590 
591  /* 6. Create a OneToOne map from the overlap map.
592  */
593 
594  // LINE 4: In the GUN paper
595 
596  RCP<const Map> non_overlap_map;
597  if(!useTieBreak_) {
598  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_04 createOneToOne");
599 
600  non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map);
601  }
602  else {
603  // use a hash tie break to get better load balancing from create one to one
604  // Aug. 4, 2016...this is a bad idea and doesn't work
605  HashTieBreak<LO,GO> tie_break;
606  non_overlap_map = Tpetra::createOneToOne<LO,GO,Node>(overlap_map,tie_break);
607  }
608 
609  /* 7. Create a non-overlapped multivector from OneToOne map.
610  */
611 
612  // LINE 5: In the GUN paper
613 
614  Teuchos::RCP<MultiVector> tagged_non_overlap_mv;
615  {
616  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_05 alloc_unique_mv");
617 
618  tagged_non_overlap_mv = Tpetra::createMultiVector<GO>(non_overlap_map,(size_t)numFields_);
619  }
620 
621  /* 8. Create an export between the two maps.
622  */
623 
624  // LINE 6: In the GUN paper
625  RCP<Export> exp;
626  RCP<Import> imp;
627  RCP<MultiVector> non_overlap_mv;
628  {
629  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_06 export");
630 
631  exp = rcp(new Export(overlap_map,non_overlap_map));
632 
633 #ifdef PANZER_DOFMGR_REQUIRE_CUDA
634  // Note: ETP 04/26/16 Temporarily create an importer for all of the
635  // doImport() calls below. This works around mysterious failures when
636  // using the exporter for Cuda builds.
637  imp = rcp(new Import(non_overlap_map,overlap_map));
638 #endif
639 
640  /* 9. Export data using ABSMAX.
641  */
642  tagged_non_overlap_mv->doExport(tagged_overlap_mv,*exp,Tpetra::ABSMAX);
643 
644  // copy the tagged one, so as to preserve the tagged MV so we can overwrite
645  // the non_overlap_mv
646  non_overlap_mv = rcp(new MultiVector(*tagged_non_overlap_mv,Teuchos::Copy));
647  }
648 
649 
650  /* 10. Compute the local sum using Kokkos.
651  */
652 
653  // LINES 7-9: In the GUN paper
654 
655  GO localsum=0;
656  {
657  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_07-09 local_count");
658 
659  typedef typename Tpetra::MultiVector<GO,Node> MV;
660  typedef typename MV::dual_view_type::t_dev KV;
661  typedef typename MV::dual_view_type::t_dev::memory_space DMS;
662  KV values = non_overlap_mv->template getLocalView<DMS>();
663  auto mv_size = values.dimension_0();
664  Kokkos::parallel_reduce(mv_size,panzer::dof_functors::SumRank2<GO,KV>(values),localsum);
665  }
666 
667  /* 11. Create a map using local sums to generate final GIDs.
668  */
669 
670  // LINE 10: In the GUN paper
671 
672  GO myOffset = -1;
673  {
674  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_10 prefix_sum");
675 
676  // do a prefix sum
677  GO scanResult = 0;
678  Teuchos::scan<int, GO> (*getComm(), Teuchos::REDUCE_SUM, static_cast<size_t> (localsum), Teuchos::outArg (scanResult));
679  myOffset = scanResult - localsum;
680  }
681 
682  // LINE 11 and 12: In the GUN paper, these steps are eliminated because
683  // the non_overlap_mv is reused
684 
685  /* 12. Iterate through the non-overlapping MV and assign GIDs to
686  * the necessary points. (Assign a -1 elsewhere.)
687  */
688 
689  // LINES 13-21: In the GUN paper
690 
691  {
692  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_13-21 gid_assignment");
693 
694  // ArrayView<const GO> owned_ids = gid_map->getNodeElementList();
695  int which_id=0;
696  ArrayRCP<ArrayRCP<GO> > editnonoverlap = non_overlap_mv->get2dViewNonConst();
697  for(size_t i=0; i<non_overlap_mv->getLocalLength(); ++i){
698  for(int j=0; j<numFields_; ++j){
699  if(editnonoverlap[j][i]!=0){
700  // editnonoverlap[j][i]=myOffset+which_id;
701  int ndof = Teuchos::as<int>(editnonoverlap[j][i]);
702  editnonoverlap[j][i]=myOffset+which_id;
703  which_id+=ndof;
704  }
705  else{
706  editnonoverlap[j][i]=-1;
707  }
708 
709  }
710  }
711  }
712 
713  // LINE 22: In the GUN paper. Were performed above, and the overlaped_mv is
714  // abused to handle input tagging.
715 
716  /* 13. Import data back to the overlap MV using REPLACE.
717  */
718 
719  // LINE 23: In the GUN paper
720 
721  {
722  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildGlobalUnknowns_GUN::line_23 final_import");
723 
724 #ifdef PANZER_DOFMGR_REQUIRE_CUDA
725  overlap_mv.doImport(*non_overlap_mv,*imp,Tpetra::REPLACE);
726 #else
727  // use exporter to save on communication setup costs
728  overlap_mv.doImport(*non_overlap_mv,*exp,Tpetra::REPLACE);
729 #endif
730  }
731 
732  //std::cout << Teuchos::describe(*non_overlap_mv,Teuchos::VERB_EXTREME) << std::endl;
733 
734  // return non_overlap_mv;
735  return std::make_pair(non_overlap_mv,tagged_non_overlap_mv);
736 }
737 
738 template <typename LO, typename GO>
739 Teuchos::RCP<Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> >
741 {
742  // some typedefs
743  typedef panzer::TpetraNodeType Node;
744  typedef Tpetra::Map<LO, GO, Node> Map;
745  typedef Tpetra::MultiVector<GO,LO,GO,Node> MultiVector;
746 
747  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector");
748 
749  //We will iterate through all of the blocks, building a FieldAggPattern for
750  //each of them.
751 
752  for (size_t b = 0; b < blockOrder_.size(); ++b) {
753  std::vector<std::pair< int, RCP<const panzer::FieldPattern> > > faConstruct;
754  //The ID is going to be the AID, and then everything will work.
755  //The ID should not be the AID, it should be the ID it has in the ordering.
756 
757  for (size_t i = 0; i < fieldAIDOrder_.size(); ++i) {
758  int looking = fieldAIDOrder_[i];
759 
760  //Check if in b's fp list
761  std::vector<int>::const_iterator reu = std::find(blockToAssociatedFP_[b].begin(), blockToAssociatedFP_[b].end(), looking);
762  if(!(reu==blockToAssociatedFP_[b].end())){
763  faConstruct.push_back(std::make_pair(i, fieldPatterns_[fieldAIDOrder_[i]]));
764  }
765 
766  }
767 
768  if(faConstruct.size()>0) {
769  fa_fps_.push_back(rcp(new FieldAggPattern(faConstruct, ga_fp_)));
770 
771  // how many global IDs are in this element block?
772  int gidsInBlock = fa_fps_[fa_fps_.size()-1]->numberIds();
773  elementBlockGIDCount_.push_back(gidsInBlock);
774  }
775  else {
776  fa_fps_.push_back(Teuchos::null);
777  elementBlockGIDCount_.push_back(0);
778  }
779  }
780 
781  RCP<const Map> overlapmap = buildOverlapMapFromElements(ownedAccess);
782 
783  // LINE 22: In the GUN paper...the overlap_mv is reused for the tagged multivector.
784  // This is a bit of a practical abuse of the algorithm presented in the paper.
785 
786  Teuchos::RCP<MultiVector> overlap_mv;
787  {
788  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector::allocate_tagged_multivector");
789 
790  overlap_mv = Tpetra::createMultiVector<GO>(overlapmap,(size_t)numFields_);
791  overlap_mv->putScalar(0); // if tpetra is not initialized with zeros
792  }
793 
794  /* 5. Iterate through all local elements again, checking with the FP
795  * information. Mark up the overlap map accordingly.
796  */
797 
798  {
799  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::buildTaggedMultiVector::fill_tagged_multivector");
800 
801  // temporary working vector to fill each row in tagged array
802  std::vector<int> working(overlap_mv->getNumVectors());
803  ArrayRCP<ArrayRCP<GO> > edittwoview = overlap_mv->get2dViewNonConst();
804  for (size_t b = 0; b < blockOrder_.size(); ++b) {
805  // there has to be a field pattern assocaited with the block
806  if(fa_fps_[b]==Teuchos::null)
807  continue;
808 
809  const std::vector<LO> & numFields= fa_fps_[b]->numFieldsPerId();
810  const std::vector<LO> & fieldIds= fa_fps_[b]->fieldIds();
811  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
812  for (size_t l = 0; l < myElements.size(); ++l) {
813  LO connSize = connMngr_->getConnectivitySize(myElements[l]);
814  const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
815  int offset=0;
816  for (int c = 0; c < connSize; ++c) {
817  size_t lid = overlapmap->getLocalElement(elmtConn[c]);
818 
819  for(std::size_t i=0;i<working.size();i++)
820  working[i] = 0;
821  for (int n = 0; n < numFields[c]; ++n) {
822  int whichField = fieldIds[offset];
823  //Row will be lid. column will be whichField.
824  //Shove onto local ordering
825  working[whichField]++;
826  offset++;
827  }
828  for(std::size_t i=0;i<working.size();i++) {
829  auto current = edittwoview[i][lid];
830  edittwoview[i][lid] = (current > working[i]) ? current : working[i];
831  }
832 
833  }
834  }
835  }
836 
837  // // verbose output for inspecting overlap_mv
838  // for(int i=0;i<overlap_mv->getLocalLength(); i++) {
839  // for(int j=0;j<overlap_mv->getNumVectors() ; j++)
840  // std::cout << edittwoview[j][i] << " ";
841  // std::cout << std::endl;
842  // }
843  }
844 
845  return overlap_mv;
846 }
847 
848 template <typename LO, typename GO>
849 int DOFManager<LO,GO>::getFieldNum(const std::string & string) const
850 {
851  int ind=0;
852  bool found=false;
853  while(!found && (size_t)ind<fieldStringOrder_.size()){
854  if(fieldStringOrder_[ind]==string)
855  found=true;
856  else
857  ind++;
858  }
859 
860  if(found)
861  return ind;
862 
863  // didn't find anything...return -1
864  return -1;
865 }
866 
867 template <typename LO, typename GO>
868 void DOFManager<LO,GO>::getFieldOrder(std::vector<std::string> & fieldOrder) const
869 {
870  fieldOrder.resize(fieldStringOrder_.size());
871  for (size_t i = 0; i < fieldStringOrder_.size(); ++i)
872  fieldOrder[i]=fieldStringOrder_[i];
873 }
874 
875 template <typename LO, typename GO>
876 bool DOFManager<LO,GO>::fieldInBlock(const std::string & field, const std::string & block) const
877 {
878  std::map<std::string,int>::const_iterator fitr = fieldNameToAID_.find(field);
879  if(fitr==fieldNameToAID_.end()) {
880  std::stringstream ss;
881  printFieldInformation(ss);
882  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid field name. DOF information is:\n"+ss.str());
883  }
884  std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(block);
885  if(bitr==blockNameToID_.end()) {
886  std::stringstream ss;
887  printFieldInformation(ss);
888  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name. DOF information is:\n"+ss.str());
889  }
890  int fid=fitr->second;
891  int bid=bitr->second;
892 
893  bool found=false;
894  for (size_t i = 0; i < blockToAssociatedFP_[bid].size(); ++i) {
895  if(blockToAssociatedFP_[bid][i]==fid){
896  found=true;
897  break;
898  }
899  }
900 
901  return found;
902 }
903 
904 template <typename LO, typename GO>
905 const std::vector<int> & DOFManager<LO,GO>::getBlockFieldNumbers(const std::string & blockId) const
906 {
907  TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,"DOFManager::getBlockFieldNumbers: BuildConnectivity must be run first.");
908 
909  std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
910  if(bitr==blockNameToID_.end())
911  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
912  int bid=bitr->second;
913 
914  // there has to be a field pattern assocaited with the block
915  if(fa_fps_[bid]!=Teuchos::null)
916  return fa_fps_[bid]->fieldIds();
917 
918  // nothing to return
919  static std::vector<int> empty;
920  return empty;
921 }
922 
923 template <typename LO, typename GO>
924 const std::pair<std::vector<int>,std::vector<int> > &
925 DOFManager<LO,GO>::getGIDFieldOffsets_closure(const std::string & blockId, int fieldNum, int subcellDim,int subcellId) const
926 {
927  TEUCHOS_TEST_FOR_EXCEPTION(!buildConnectivityRun_,std::logic_error,"DOFManager::getGIDFieldOffsets_closure: BuildConnectivity must be run first.");
928  std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
929  if(bitr==blockNameToID_.end())
930  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, "DOFManager::getGIDFieldOffsets_closure: invalid block name.");
931 
932  // there has to be a field pattern assocaited with the block
933  if(fa_fps_[bitr->second]!=Teuchos::null)
934  return fa_fps_[bitr->second]->localOffsets_closure(fieldNum, subcellDim, subcellId);
935 
936  static std::pair<std::vector<int>,std::vector<int> > empty;
937  return empty;
938 }
939 
940 template <typename LO, typename GO>
941 void DOFManager<LO,GO>::ownedIndices(const std::vector<GO> & indices,std::vector<bool> & isOwned) const
942 {
943  //Resizes the isOwned array.
944  if(indices.size()!=isOwned.size())
945  isOwned.resize(indices.size(),false);
946  typename std::vector<GO>::const_iterator endOf = owned_.end();
947  for (std::size_t i = 0; i < indices.size(); ++i) {
948  isOwned[i] = ( std::find(owned_.begin(), owned_.end(), indices[i])!=endOf );
949  }
950 }
951 
952 template <typename LO, typename GO>
953 void DOFManager<LO,GO>::setFieldOrder(const std::vector<std::string> & fieldOrder)
954 {
955  TEUCHOS_TEST_FOR_EXCEPTION(buildConnectivityRun_,std::logic_error,
956  "DOFManager::setFieldOrder: setFieldOrder cannot be called after "
957  "buildGlobalUnknowns has been called");
958  if(validFieldOrder(fieldOrder)){
959  fieldStringOrder_=fieldOrder;
960  //We also need to reassign the fieldAIDOrder_.
961  for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
962  fieldAIDOrder_[i]=fieldNameToAID_[fieldStringOrder_[i]];
963  }
964  }
965  else
966  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::setFieldOrder: Invalid Field Ordering!");
967 
968 }
969 
970 
971 template <typename LO, typename GO>
972 bool DOFManager<LO,GO>::validFieldOrder(const std::vector<std::string> & proposed_fieldOrder)
973 {
974  if(fieldStringOrder_.size()!=proposed_fieldOrder.size())
975  return false;
976  //I'm using a not very efficient way of doing this, but there should never
977  //be that many fields, so it isn't that big of a deal.
978  for (size_t i = 0; i < fieldStringOrder_.size(); ++i) {
979  bool found=false;
980  for (size_t j = 0; j < proposed_fieldOrder.size(); ++j) {
981  if(fieldStringOrder_[i]==proposed_fieldOrder[j])
982  found=true;
983  }
984  if(!found)
985  return false;
986  }
987  return true;
988 }
989 
990 template <typename LO, typename GO>
991 const std::string & DOFManager<LO,GO>::getFieldString(int num) const
992 {
993  //A reverse lookup through fieldStringOrder_.
994  if(num>=(int)fieldStringOrder_.size())
995  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "DOFManager::getFieldString: invalid number");
996  return fieldStringOrder_[num];
997 }
998 
999 //Everything associated with orientation is not yeat built, but this
1000 //is the method as found in Panzer_DOFManager_impl.hpp. There are
1001 //going to need to be some substantial changes to the code as it applies
1002 //to this DOFManager, but the basic ideas and format should be similar.
1003 //
1004 template <typename LO, typename GO>
1006 {
1007  orientation_.clear(); // clean up previous work
1008 
1009  std::vector<std::string> elementBlockIds;
1010  connMngr_->getElementBlockIds(elementBlockIds);
1011 
1012  // figure out how many total elements are owned by this processor (why is this so hard!)
1013  std::size_t myElementCount = 0;
1014  for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin(); blockItr!=elementBlockIds.end();++blockItr)
1015  myElementCount += connMngr_->getElementBlock(*blockItr).size();
1016 
1017  // allocate for each block
1018  orientation_.resize(myElementCount);
1019 
1020  // loop over all element blocks
1021  for(std::vector<std::string>::const_iterator blockItr=elementBlockIds.begin();
1022  blockItr!=elementBlockIds.end();++blockItr) {
1023  const std::string & blockName = *blockItr;
1024 
1025  // this block has no unknowns (or elements)
1026  std::map<std::string,int>::const_iterator fap = blockNameToID_.find(blockName);
1027  if(fap==blockNameToID_.end()) {
1028  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::buildUnknownsOrientation: invalid block name");
1029  }
1030 
1031  int bid=fap->second;
1032 
1033  if(fa_fps_[bid]==Teuchos::null)
1034  continue;
1035 
1036  // grab field patterns, will be necessary to compute orientations
1037  const FieldPattern & fieldPattern = *fa_fps_[bid];
1038 
1039  //Should be ga_fp_ (geometric aggregate field pattern)
1040  std::vector<std::pair<int,int> > topEdgeIndices;
1041  orientation_helpers::computePatternEdgeIndices(*ga_fp_,topEdgeIndices);
1042 
1043  // grab face orientations if 3D
1044  std::vector<std::vector<int> > topFaceIndices;
1045  if(ga_fp_->getDimension()==3)
1046  orientation_helpers::computePatternFaceIndices(*ga_fp_,topFaceIndices);
1047 
1048  //How many GIDs are associated with a particular element bloc
1049  const std::vector<LO> & elmts = connMngr_->getElementBlock(blockName);
1050  for(std::size_t e=0;e<elmts.size();e++) {
1051  // this is the vector of orientations to fill: initialize it correctly
1052  std::vector<char> & eOrientation = orientation_[elmts[e]];
1053 
1054  // This resize seems to be the same as fieldPattern.numberIDs().
1055  // When computer edge orientations is called, that is the assert.
1056  // There should be no reason to make it anymore complicated.
1057  eOrientation.resize(fieldPattern.numberIds());
1058  for(std::size_t s=0;s<eOrientation.size();s++)
1059  eOrientation[s] = 1; // put in 1 by default
1060 
1061  // get geometry ids
1062  LO connSz = connMngr_->getConnectivitySize(elmts[e]);
1063  const GO * connPtr = connMngr_->getConnectivity(elmts[e]);
1064  const std::vector<GO> connectivity(connPtr,connPtr+connSz);
1065 
1066  orientation_helpers::computeCellEdgeOrientations(topEdgeIndices, connectivity, fieldPattern, eOrientation);
1067 
1068  // compute face orientations in 3D
1069  if(ga_fp_->getDimension()==3)
1070  orientation_helpers::computeCellFaceOrientations(topFaceIndices, connectivity, fieldPattern, eOrientation);
1071  }
1072  }
1073 }
1074 
1075 template <typename LO, typename GO>
1076 void DOFManager<LO,GO>::getElementOrientation(LO localElmtId,std::vector<double> & gidsOrientation) const
1077 {
1078  TEUCHOS_TEST_FOR_EXCEPTION(orientation_.size()==0,std::logic_error,
1079  "DOFManager::getElementOrientations: Orientations were not constructed!");
1080 
1081  const std::vector<char> & local_o = orientation_[localElmtId];
1082  gidsOrientation.resize(local_o.size());
1083  for(std::size_t i=0;i<local_o.size();i++) {
1084  gidsOrientation[i] = double(local_o[i]);
1085  }
1086 }
1087 
1088 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1089 Teuchos::RCP<ConnManager<LocalOrdinalT,GlobalOrdinalT> > DOFManager<LocalOrdinalT,GlobalOrdinalT>::resetIndices()
1090 {
1091  Teuchos::RCP<ConnManager<LocalOrdinalT,GlobalOrdinalT> > connMngr = connMngr_;
1092 
1093  connMngr_ = Teuchos::null;
1094 
1095  // wipe out FEI objects
1096  orientation_.clear(); // clean up previous work
1097  fa_fps_.clear();
1098  elementGIDs_.clear();
1099  owned_.clear();
1100  ghosted_.clear();
1101  elementBlockGIDCount_.clear();
1102 
1103  return connMngr;
1104 }
1105 
1106 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1107 std::size_t DOFManager<LocalOrdinalT,GlobalOrdinalT>::blockIdToIndex(const std::string & blockId) const
1108 {
1109  std::map<std::string,int>::const_iterator bitr = blockNameToID_.find(blockId);
1110  if(bitr==blockNameToID_.end())
1111  TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"DOFManager::fieldInBlock: invalid block name");
1112  return bitr->second;
1113 }
1114 
1115 template <typename LocalOrdinalT,typename GlobalOrdinalT>
1117 {
1118  os << "DOFManager Field Information: " << std::endl;
1119 
1120  // sanity check
1121  TEUCHOS_ASSERT(blockOrder_.size()==blockToAssociatedFP_.size());
1122 
1123  for(std::size_t i=0;i<blockOrder_.size();i++) {
1124  os << " Element Block = " << blockOrder_[i] << std::endl;
1125 
1126  // output field information
1127  const std::vector<int> & fieldIds = blockToAssociatedFP_[i];
1128  for(std::size_t f=0;f<fieldIds.size();f++)
1129  os << " \"" << getFieldString(fieldIds[f]) << "\" is field ID " << fieldIds[f] << std::endl;
1130  }
1131 }
1132 
1133 template <typename LO,typename GO>
1134 Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> >
1137 {
1138  PANZER_DOFMGR_FUNC_TIME_MONITOR("panzer::DOFManager::builderOverlapMapFromElements");
1139 
1140  /*
1141  * 2. Iterate through all local elements and create the overlapVector
1142  * of concerned elements.
1143  */
1144 
1145  std::set<GO> overlapset;
1146  for (size_t i = 0; i < blockOrder_.size(); ++i) {
1147  const std::vector<LO> & myElements = access.getElementBlock(blockOrder_[i]);
1148  for (size_t e = 0; e < myElements.size(); ++e) {
1149  LO connSize = connMngr_->getConnectivitySize(myElements[e]);
1150  const GO * elmtConn = connMngr_->getConnectivity(myElements[e]);
1151  for (int k = 0; k < connSize; ++k) {
1152  overlapset.insert(elmtConn[k]);
1153  }
1154  }
1155  }
1156 
1157  Array<GO> overlapVector;
1158  for (typename std::set<GO>::const_iterator itr = overlapset.begin(); itr!=overlapset.end(); ++itr) {
1159  overlapVector.push_back(*itr);
1160  }
1161 
1162  /* 3. Construct an overlap map from this structure.
1163  */
1164  return Tpetra::createNonContigMap<LO,GO>(overlapVector,getComm());
1165 }
1166 
1167 template <typename LO,typename GO>
1170  std::vector<std::vector< GO > > & elementGIDs,
1171  const Tpetra::Map<LO,GO,panzer::TpetraNodeType> & overlapmap,
1172  const Tpetra::MultiVector<GO,LO,GO,panzer::TpetraNodeType> & overlap_mv) const
1173 {
1174  using Teuchos::ArrayRCP;
1175 
1176  //To generate elementGIDs we need to go through all of the local elements.
1177  ArrayRCP<ArrayRCP<const GO> > twoview = overlap_mv.get2dView();
1178 
1179  //And for each of the things in fa_fp.fieldIds we go to that column. To the the row,
1180  //we move from globalID to localID in the map and use our local value for something.
1181  for (size_t b = 0; b < blockOrder_.size(); ++b) {
1182  const std::vector<LO> & myElements = access.getElementBlock(blockOrder_[b]);
1183 
1184  if(fa_fps_[b]==Teuchos::null) {
1185  // fill elements that are not used with empty vectors
1186  for (size_t l = 0; l < myElements.size(); ++l) {
1187  LO thisID=myElements[l];
1188  if(elementGIDs.size()<=(size_t)thisID)
1189  elementGIDs.resize(thisID+1);
1190  }
1191  continue;
1192  }
1193 
1194  const std::vector<int> & numFields= fa_fps_[b]->numFieldsPerId();
1195  const std::vector<int> & fieldIds= fa_fps_[b]->fieldIds();
1196  //
1197  //
1198  for (size_t l = 0; l < myElements.size(); ++l) {
1199  LO connSize = connMngr_->getConnectivitySize(myElements[l]);
1200  const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
1201  std::vector<GO> localOrdering;
1202  int offset=0;
1203  for (int c = 0; c < connSize; ++c) {
1204  size_t lid = overlapmap.getLocalElement(elmtConn[c]);
1205  std::vector<int> dofsPerField(numFields_,0);
1206  for (int n = 0; n < numFields[c]; ++n) {
1207  int whichField = fieldIds[offset];
1208  offset++;
1209  //Row will be lid. column will be whichField.
1210  //Shove onto local ordering
1211  localOrdering.push_back(twoview[whichField][lid]+dofsPerField[whichField]);
1212 
1213  dofsPerField[whichField]++;
1214  }
1215  }
1216  LO thisID=myElements[l];
1217  if(elementGIDs.size()<=(size_t)thisID){
1218  elementGIDs.resize(thisID+1);
1219  }
1220  elementGIDs[thisID]=localOrdering;
1221  }
1222  }
1223 }
1224 
1225 template <typename LO, typename GO>
1227 {
1228  std::vector<std::vector<LO> > elementLIDs(elementGIDs_.size());
1229 
1230  std::vector<GO> ownedAndGhosted;
1231  this->getOwnedAndGhostedIndices(ownedAndGhosted);
1232 
1233  // build global to local hash map (temporary and used only once)
1234  std::unordered_map<GO,LO> hashMap;
1235  for(std::size_t i = 0; i < ownedAndGhosted.size(); ++i)
1236  hashMap[ownedAndGhosted[i]] = i;
1237 
1238  for (std::size_t i = 0; i < elementGIDs_.size(); ++i) {
1239  const std::vector<GO>& gids = elementGIDs_[i];
1240  std::vector<LO>& lids = elementLIDs[i];
1241  lids.resize(gids.size());
1242  for (std::size_t g = 0; g < gids.size(); ++g)
1243  lids[g] = hashMap[gids[g]];
1244  }
1245 
1246  this->setLocalIds(elementLIDs);
1247 }
1248 
1249 /*
1250 template <typename LO,typename GO>
1251 Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> >
1252 DOFManager<LO,GO>::runLocalRCMReordering(const Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType> > & map)
1253 {
1254  typedef panzer::TpetraNodeType Node;
1255  typedef Tpetra::Map<LO, GO, Node> Map;
1256  typedef Tpetra::CrsGraph<LO, GO, Node> Graph;
1257 
1258  Teuchos::RCP<Graph> graph = Teuchos::rcp(new Graph(map,0));
1259 
1260  // build Crs Graph from the mesh
1261  for (size_t b = 0; b < blockOrder_.size(); ++b) {
1262  if(fa_fps_[b]==Teuchos::null)
1263  continue;
1264 
1265  const std::vector<LO> & myElements = connMngr_->getElementBlock(blockOrder_[b]);
1266  for (size_t l = 0; l < myElements.size(); ++l) {
1267  LO connSize = connMngr_->getConnectivitySize(myElements[l]);
1268  const GO * elmtConn = connMngr_->getConnectivity(myElements[l]);
1269  for (int c = 0; c < connSize; ++c) {
1270  LO lid = map->getLocalElement(elmtConn[c]);
1271  if(Teuchos::OrdinalTraits<LO>::invalid()!=lid)
1272  continue;
1273 
1274  graph->insertGlobalIndices(elmtConn[c],Teuchos::arrayView(elmtConn,connSize));
1275  }
1276  }
1277  }
1278 
1279  graph->fillComplete();
1280 
1281  std::vector<GO> newOrder(map->getNodeNumElements());
1282  {
1283  // graph is constructed, now run RCM using zoltan2
1284  typedef Zoltan2::XpetraCrsGraphInput<Graph> SparseGraphAdapter;
1285 
1286  Teuchos::ParameterList params;
1287  params.set("order_method", "rcm");
1288  SparseGraphAdapter adapter(graph);
1289 
1290  Zoltan2::OrderingProblem<SparseGraphAdapter> problem(&adapter, &params,MPI_COMM_SELF);
1291  problem.solve();
1292 
1293  // build a new global ording array using permutation
1294  Zoltan2::OrderingSolution<GO,LO> * soln = problem.getSolution();
1295 
1296  size_t dummy;
1297  size_t checkLength = soln->getPermutationSize();
1298  LO * checkPerm = soln->getPermutation(&dummy);
1299 
1300  Teuchos::ArrayView<const GO > oldOrder = map->getNodeElementList();
1301  TEUCHOS_ASSERT(checkLength==oldOrder.size());
1302  TEUCHOS_ASSERT(checkLength==newOrder.size());
1303 
1304  for(size_t i=0;i<checkLength;i++)
1305  newOrder[checkPerm[i]] = oldOrder[i];
1306  }
1307 
1308  return Tpetra::createNonContigMap<LO,GO>(newOrder,communicator_);
1309 }
1310 */
1311 
1312 } /*panzer*/
1313 
1314 #endif
const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const
void computeCellFaceOrientations(const std::vector< std::pair< int, int > > &topEdgeIndices, const std::vector< GlobalOrdinalT > &topology, const FieldPattern &fieldPattern, std::vector< 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
const std::string & getFieldString(int num) const
Reverse lookup of the field string from a field number.
void fillGIDsFromOverlappedMV(const ElementBlockAccess &access, std::vector< std::vector< GO > > &elementGIDs, const Tpetra::Map< LO, GO, panzer::TpetraNodeType > &overlapmap, const Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &overlap_mv) const
#define PANZER_DOFMGR_FUNC_TIME_MONITOR(a)
void setConnManager(const Teuchos::RCP< ConnManager< LO, GO > > &connMngr, MPI_Comm mpiComm)
Adds a Connection Manager that will be associated with this DOFManager.
const std::vector< int > & getGIDFieldOffsets(const std::string &blockID, int fieldNum) const
Teuchos::RCP< const Tpetra::Map< LO, GO, panzer::TpetraNodeType > > buildOverlapMapFromElements(const ElementBlockAccess &access) const
PHX::MDField< ScalarT > vector
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
void computeCellEdgeOrientations(const std::vector< std::pair< int, int > > &topEdgeIndices, const std::vector< GlobalOrdinalT > &topology, const FieldPattern &fieldPattern, std::vector< char > &orientation)
std::size_t blockIdToIndex(const std::string &blockId) const
Teuchos::RCP< const FieldPattern > getFieldPattern(const std::string &name) const
Find a field pattern stored for a particular block and field number. This will retrive the pattern ad...
void setFieldOrder(const std::vector< std::string > &fieldOrder)
void getOwnedAndGhostedIndices(std::vector< GlobalOrdinalT > &indices) const
std::vector< ScalarT > values
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
int addField(const std::string &str, const Teuchos::RCP< const FieldPattern > &pattern)
Add a field to the DOF manager.
void buildGlobalUnknowns()
builds the global unknowns array
PHX::MDField< const ScalarT, Cell, IP > field
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
std::pair< Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > >, Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > > > buildGlobalUnknowns_GUN(const Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &tagged_overlap_mv, Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > &overlap_mv) const
Teuchos::RCP< Tpetra::MultiVector< GO, LO, GO, panzer::TpetraNodeType > > buildTaggedMultiVector(const ElementBlockAccess &access)
void printFieldInformation(std::ostream &os) const
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 std::pair< std::vector< int >, std::vector< int > > & getGIDFieldOffsets_closure(const std::string &blockId, int fieldNum, int subcellDim, int subcellId) const
Use the field pattern so that you can find a particular field in the GIDs array. This version lets yo...
void computePatternEdgeIndices(const FieldPattern &pattern, std::vector< std::pair< int, int > > &edgeIndices)
void ownedIndices(const std::vector< GlobalOrdinalT > &indices, std::vector< bool > &isOwned) const