Panzer  Version of the Day
Panzer_Filtered_UniqueGlobalIndexer_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_Filtered_UniqueGlobalIndexer_impl_hpp__
44 #define __Panzer_Filtered_UniqueGlobalIndexer_impl_hpp__
45 
46 #include <unordered_set>
47 
48 #include "PanzerDofMgr_config.hpp"
49 #include "Panzer_NodeType.hpp"
50 
51 #include "Tpetra_Map.hpp"
52 #include "Tpetra_Import.hpp"
53 #include "Tpetra_Export.hpp"
54 #include "Tpetra_Vector.hpp"
55 
56 namespace panzer {
57 
58 template <typename LocalOrdinalT,typename GlobalOrdinalT>
61 { }
62 
64 //
65 // initialize()
66 //
68 template<typename LocalOrdinalT, typename GlobalOrdinalT>
69 void
73  const std::vector<GlobalOrdinalT>& filtered)
74 {
75  typedef GlobalOrdinalT GO;
76  typedef LocalOrdinalT LO;
77  typedef panzer::TpetraNodeType Node;
78  typedef Tpetra::Map<LO, GO, Node> Map;
79  typedef Tpetra::Vector<GO,LO,GO,Node> Vector;
80  typedef Tpetra::Export<LO,GO,Node> Export;
81 
82  using std::size_t;
83  using std::vector;
84  using HashTable = std::unordered_set<GlobalOrdinalT>;
85  using Teuchos::RCP;
86 
87  owned_.clear();
88  ghosted_.clear();
89  base_ = ugi;
90 
91  // From the base global indexer, build the filtered owned indices.
92  vector<GlobalOrdinalT> baseOwned, baseGhosted;
93  base_->getOwnedIndices(baseOwned);
94  base_->getGhostedIndices(baseGhosted);
95 
96  RCP<const Map> ownedMap
97  = Tpetra::createNonContigMap<LO,GO>(baseOwned,getComm());
98  RCP<const Map> ghostedMap
99  = Tpetra::createNonContigMap<LO,GO>(baseGhosted,getComm());
100 
101  Vector ownedFiltered(ownedMap);
102  Vector ghostedFiltered(ghostedMap);
103 
104  ownedFiltered.putScalar(0.0);
105  ghostedFiltered.putScalar(0.0);
106 
107  for(GlobalOrdinalT f : filtered) {
108  bool isOwned = std::find(baseOwned.begin(),baseOwned.end(),f)!=baseOwned.end();
109  bool isGhosted = std::find(baseGhosted.begin(),baseGhosted.end(),f)!=baseGhosted.end();
110 
111  if(isOwned)
112  ownedFiltered.replaceGlobalValue(f,1.0);
113  else if(isGhosted)
114  ghostedFiltered.replaceGlobalValue(f,1.0);
115  // else no one cares...
116  }
117 
118  Export exporter(ghostedMap,ownedMap);
119  ownedFiltered.doExport(ghostedFiltered, exporter, Tpetra::ADD);
120 
121  Teuchos::ArrayRCP<const GlobalOrdinalT> data = ownedFiltered.getData();
122 
123  // Build a hash table for fast searching.
124  HashTable filteredHash;
125  for(int i(0); i < data.size(); ++i) {
126  if(data[i]!=0)
127  filteredHash.insert(baseOwned[i]);
128  }
129  // for (size_t i(0); i < filtered.size(); ++i) {
130  // filteredHash.insert(filtered[i]);
131  // }
132 
133  // Search for indices in the filtered array; add to owned_ if not found, and
134  // add to ghosted_ otherwise.
135  for (size_t i(0); i < baseOwned.size(); ++i)
136  {
137  auto itr = filteredHash.find(baseOwned[i]);
138  if (itr == filteredHash.end())
139  owned_.push_back(baseOwned[i]);
140  else
141  ghosted_.push_back(baseOwned[i]);
142  }
143  ghosted_.insert(ghosted_.end(), baseGhosted.begin(), baseGhosted.end());
144 
145  // Now that we've change the owned_ and ghosted_ vectors, we need to rebuild
146  // the local IDs.
147  this->buildLocalIds();
148 } // end of initialize()
149 
150 template <typename LocalOrdinalT,typename GlobalOrdinalT>
151 void
153 getOwnedAndGhostedNotFilteredIndicator(std::vector<int> & indicator) const
154 {
155  using Teuchos::RCP;
156 
157  typedef GlobalOrdinalT GO;
158  typedef LocalOrdinalT LO;
159  typedef panzer::TpetraNodeType Node;
160  typedef Tpetra::Map<LO, GO, Node> Map;
161  typedef Tpetra::Vector<GO,LO,GO,Node> Vector;
162  typedef Tpetra::Import<LO,GO,Node> Import;
163 
164  std::vector<GlobalOrdinalT> ownedIndices;
165  std::vector<GlobalOrdinalT> ghostedIndices;
166 
167  // build owned and ghosted maps
168  getOwnedIndices(ownedIndices);
169  getOwnedAndGhostedIndices(ghostedIndices);
170 
171  RCP<const Map> ownedMap
172  = Tpetra::createNonContigMap<LO,GO>(ownedIndices,getComm());
173  RCP<const Map> ghostedMap
174  = Tpetra::createNonContigMap<LO,GO>(ghostedIndices,getComm());
175 
176  // allocate the owned vector, mark those GIDs as unfiltered
177  // (they are by definition)
178  Vector ownedActive(ownedMap);
179  ownedActive.putScalar(1);
180 
181  // Initialize all indices to zero
182  Vector ghostedActive(ghostedMap);
183  ghostedActive.putScalar(0);
184 
185  // do communication, marking unfiltered indices as 1 (filtered
186  // indices locally are marked as zero)
187  Import importer(ownedMap,ghostedMap);
188  ghostedActive.doImport(ownedActive,importer,Tpetra::INSERT);
189 
190  Teuchos::ArrayRCP<const GO> data = ghostedActive.getData();
191 
192  // copy communicated data (clear it out first)
193  indicator.clear();
194  indicator.insert(indicator.end(),data.begin(),data.end());
195 }
196 
197 template <typename LocalOrdinalT,typename GlobalOrdinalT>
198 void
200 getFilteredOwnedAndGhostedIndices(std::vector<GlobalOrdinalT> & indices) const
201 {
202  using Teuchos::RCP;
203 
204  // get filtered/unfiltered indicator vector
205  std::vector<int> indicators;
206  getOwnedAndGhostedNotFilteredIndicator(indicators);
207 
208  // build ghosted maps
209  std::vector<GlobalOrdinalT> ghostedIndices;
210  getOwnedAndGhostedIndices(ghostedIndices);
211 
212  // filtered out filtered indices (isn't that a useful comment)
213  for(std::size_t i=0;i<indicators.size();i++) {
214  if(indicators[i]==1)
215  indices.push_back(ghostedIndices[i]);
216  }
217 }
218 
219 template <typename LocalOrdinalT,typename GlobalOrdinalT>
220 void
222 ownedIndices(const std::vector<GlobalOrdinalT> & indices,std::vector<bool> & isOwned) const
223 {
224  //Resizes the isOwned array.
225  if(indices.size()!=isOwned.size())
226  isOwned.resize(indices.size(),false);
227  typename std::vector<GlobalOrdinalT>::const_iterator endOf = owned_.end();
228  for (std::size_t i = 0; i < indices.size(); ++i) {
229  isOwned[i] = ( std::find(owned_.begin(), owned_.end(), indices[i])!=endOf );
230  }
231 }
232 
233 }
234 
235 #endif
void initialize(const Teuchos::RCP< const UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > &ugi, const std::vector< GlobalOrdinalT > &filteredIndices)
size_type size() const
void getOwnedAndGhostedNotFilteredIndicator(std::vector< int > &indicator) const
Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > TpetraNodeType
iterator end() const
virtual void ownedIndices(const std::vector< GlobalOrdinalT > &indices, std::vector< bool > &isOwned) const
iterator begin() const
void getFilteredOwnedAndGhostedIndices(std::vector< GlobalOrdinalT > &indices) const