Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
Map_Xpetra_def.hpp
1#ifndef MAP_XPETRA_DEF_hpp
2#define MAP_XPETRA_DEF_hpp
3#include "Map_Xpetra_decl.hpp"
4
13
14namespace FEDD {
15template < class LO, class GO, class NO>
16Map_Xpetra<LO,GO,NO>::Map_Xpetra():
17map_()
18{
19
20}
21
22template <class LO, class GO, class NO>
23Map_Xpetra<LO,GO,NO>::Map_Xpetra( XpetraMapConstPtrConst_Type& xpetraMapPtrIn ):
24map_( xpetraMapPtrIn )
25{
26
27}
28
29template < class LO, class GO, class NO>
30Map_Xpetra<LO,GO,NO>::Map_Xpetra( const Map_Type& mapIn ):
31map_()
32{
33 Xpetra::UnderlyingLib ulib;
34 if (!mapIn.getUnderlyingLib().compare("Epetra"))
35 ulib = Xpetra::UseEpetra;
36 else if (!mapIn.getUnderlyingLib().compare("Tpetra"))
37 ulib = Xpetra::UseTpetra;
38
39 map_ = Xpetra::MapFactory<LO,GO,NO>::Build( ulib, mapIn.getGlobalNumElements(), mapIn.getNodeElementList(), mapIn.getIndexBase(), mapIn.getComm() );
40}
41
42template < class LO, class GO, class NO>
43Map_Xpetra<LO,GO,NO>::Map_Xpetra(std::string lib,
44 GO numGlobalElements,
45 const Teuchos::ArrayView<const GO> &elementList,
46 GO indexBase,
47 const CommConstPtr_Type &comm):
48map_()
49{
50 Xpetra::UnderlyingLib ulib;
51 if (!lib.compare("Epetra"))
52 ulib = Xpetra::UseEpetra;
53 else if (!lib.compare("Tpetra"))
54 ulib = Xpetra::UseTpetra;
55 map_ = Xpetra::MapFactory<LO,GO,NO>::Build( ulib, numGlobalElements, elementList, indexBase, comm);
56}
57
58template < class LO, class GO, class NO>
59Map_Xpetra<LO,GO,NO>::Map_Xpetra(std::string lib,
60 GO numGlobalElements,
61 LO numLocalElements,
62 GO indexBase,
63 const CommConstPtr_Type &comm):
64map_()
65{
66 Xpetra::UnderlyingLib ulib;
67 if (!lib.compare("Epetra"))
68 ulib = Xpetra::UseEpetra;
69 else if (!lib.compare("Tpetra"))
70 ulib = Xpetra::UseTpetra;
71
72 map_ = Xpetra::MapFactory<LO,GO,NO>::Build( ulib, numGlobalElements, numLocalElements, indexBase, comm);
73}
74
75template < class LO, class GO, class NO>
76Map_Xpetra<LO,GO,NO>::~Map_Xpetra(){
77
78}
79
80template < class LO, class GO, class NO>
81std::string Map_Xpetra<LO,GO,NO>::getUnderlyingLib( ) const{
82 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
83 std::string uLib;
84 if (map_->lib() == Xpetra::UseEpetra)
85 uLib = "Epetra";
86 else if (map_->lib() == Xpetra::UseTpetra)
87 uLib = "Tpetra";
88 else
89 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Unknown underlying lib");
90 return uLib;
91}
92
93template < class LO, class GO, class NO>
94typename Map_Xpetra<LO,GO,NO>::MapPtr_Type Map_Xpetra<LO,GO,NO>::buildVecFieldMap(UN numDofs, std::string ordering) const{
95
96 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
97 TEUCHOS_TEST_FOR_EXCEPTION(ordering.compare("NodeWise"), std::logic_error,"Select a valid ordering: NodeWise");
98 Teuchos::ArrayView<const GO> elementList = map_->getLocalElementList();
99 Teuchos::Array<GO> elementListField( numDofs * elementList.size() );
100 for (UN i=0; i<elementList.size(); i++) {
101 for (UN dof=0; dof<numDofs; dof++)
102 elementListField[ numDofs * i + dof ] = numDofs * elementList[i] + dof;
103 }
104 typedef Teuchos::OrdinalTraits<GO> GOOT;
105 return Teuchos::rcp(new Map_Type( this->getUnderlyingLib(), GOOT::invalid(), elementListField(), map_->getIndexBase(), map_->getComm() ) );
106}
107
108template < class LO, class GO, class NO>
109typename Map_Xpetra<LO,GO,NO>::XpetraMapConstPtr_Type Map_Xpetra<LO,GO,NO>::getXpetraMap() const{
110
111 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"getXpetraMap(): map_ is null.");
112
113 return map_;
114}
115template < class LO, class GO, class NO>
116typename Map_Xpetra<LO,GO,NO>::ThyraVSBConstPtr_Type Map_Xpetra<LO,GO,NO>::getThyraVectorSpaceBase() const{
117
118 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"getThyraVectorSpaceBase(): map_ is null.");
119 return Xpetra::ThyraUtils<default_sc,LO,GO,NO>::toThyra( map_ );
120}
121
122
123template < class LO, class GO, class NO>
124GO Map_Xpetra<LO,GO,NO>::getGlobalElement(LO id) const{
125 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
126 return map_->getGlobalElement(id);
127}
128
129template < class LO, class GO, class NO>
130LO Map_Xpetra<LO,GO,NO>::getLocalElement(GO id) const{
131 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
132 return map_->getLocalElement(id);
133}
134
135template < class LO, class GO, class NO>
136LO Map_Xpetra<LO,GO,NO>::getNodeNumElements() const{
137 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
138 return map_->getLocalNumElements();
139}
140
141template < class LO, class GO, class NO>
142GO Map_Xpetra<LO,GO,NO>::getGlobalNumElements() const{
143 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
144 return map_->getGlobalNumElements();
145}
146
147template < class LO, class GO, class NO>
148Teuchos::ArrayView< const GO > Map_Xpetra<LO,GO,NO>::getNodeElementList() const{
149 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
150 return map_->getLocalElementList();
151}
152
153template < class LO, class GO, class NO>
154GO Map_Xpetra<LO,GO,NO>::getMaxAllGlobalIndex() const{
155 return map_->getMaxAllGlobalIndex();
156}
157
158template < class LO, class GO, class NO>
159LO Map_Xpetra<LO,GO,NO>::getMaxLocalIndex() const{
160 return map_->getMaxLocalIndex();
161}
162
163template < class LO, class GO, class NO>
164void Map_Xpetra<LO,GO,NO>::print(Teuchos::EVerbosityLevel verbLevel) const{
165
166 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
167 map_->describe(*out,verbLevel);
168}
169
170template < class LO, class GO, class NO>
171typename Map_Xpetra<LO,GO,NO>::CommConstPtr_Type Map_Xpetra<LO,GO,NO>::getComm() const{
172
173 return map_->getComm();
174}
175
176template < class LO, class GO, class NO>
177typename Map_Xpetra<LO,GO,NO>::CommPtr_Type Map_Xpetra<LO,GO,NO>::getCommNonConst() {
178 CommConstPtr_Type commConst = map_->getComm();
179 return Teuchos::rcp_const_cast<Comm_Type>(commConst);
180}
181
182template <class LO,class GO,class NO>
183Teuchos::RCP<Map_Xpetra<LO,GO,NO> > Map_Xpetra<LO,GO,NO>::buildUniqueMap( int numFreeProcs ) const
184{
185 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
186 if (numFreeProcs==0) {
187
188 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > myIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(map_);
189 myIndices->putScalar(map_->getComm()->getRank()+1);
190
191 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > linearMap = Xpetra::MapFactory<LO,GO,NO>::Build(map_->lib(),map_->getMaxAllGlobalIndex()+1,0,map_->getComm());
192 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > globalIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(linearMap);
193
194 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer = Xpetra::ImportFactory<LO,GO,NO>::Build(map_,linearMap);
195 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer2 = Xpetra::ImportFactory<LO,GO,NO>::Build(linearMap,map_);
196 globalIndices->doImport(*myIndices,*importer,Xpetra::INSERT);
197 myIndices->putScalar(0);
198 myIndices->doImport(*globalIndices,*importer2,Xpetra::ADD);
199
200 Teuchos::Array<GO> uniqueVector;
201 for (unsigned i=0; i<myIndices->getLocalLength(); i++) {
202 if (myIndices->getData(0)[i] == map_->getComm()->getRank()+1) {
203 uniqueVector.push_back(map_->getGlobalElement(i));
204 }
205 }
206 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapXpetra = Xpetra::MapFactory<LO,GO,NO>::Build(map_->lib(),-1,uniqueVector(),0,map_->getComm());
207 Teuchos::RCP<Map_Xpetra<LO,GO,NO> > map = Teuchos::rcp( new Map_Xpetra<LO,GO,NO>( mapXpetra ) );
208 return map;
209 }
210 else{
211 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > myIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(map_);
212 myIndices->putScalar(map_->getComm()->getRank()+1);
213 GO maxGID = map_->getMaxAllGlobalIndex();
214 int numAvailableRanks = map_->getComm()->getSize() - numFreeProcs;
215 int numElementsForAvailRank = (int) ( ( (maxGID+1) / numAvailableRanks ) + 100.*std::numeric_limits<double>::epsilon() );
216
217 int remainingElement = maxGID+1 - numAvailableRanks * numElementsForAvailRank;
218 bool hasOneMoreElement = false;
219
220 if ( remainingElement > map_->getComm()->getRank() ) {
221 numElementsForAvailRank++;
222 hasOneMoreElement = true;
223 }
224
225 if ( map_->getComm()->getRank() + 1 > map_->getComm()->getSize() - numFreeProcs){
226 numElementsForAvailRank = 0;
227 }
228
229 Teuchos::Array<GO> myElements( numElementsForAvailRank );
230 GO offset = numElementsForAvailRank * map_->getComm()->getRank();
231 if (!hasOneMoreElement) {
232 offset += remainingElement;
233 }
234 for (int i=0; i<myElements.size(); i++) {
235 myElements[i] = i + offset;
236 }
237
238 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > linearMapAvailRanks = Xpetra::MapFactory<LO,GO,NO>::Build( map_->lib(),Teuchos::OrdinalTraits<GO>::invalid(), myElements(), 0, map_->getComm() );
239
240 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > globalIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(linearMapAvailRanks);
241
242 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer = Xpetra::ImportFactory<LO,GO,NO>::Build( map_, linearMapAvailRanks );
243 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer2 = Xpetra::ImportFactory<LO,GO,NO>::Build( linearMapAvailRanks, map_ );
244
245 globalIndices->doImport(*myIndices,*importer,Xpetra::INSERT);
246
247 myIndices->putScalar(0);
248 myIndices->doImport(*globalIndices,*importer2,Xpetra::ADD);
249
250 Teuchos::Array<GO> uniqueVector;
251 for (unsigned i=0; i<myIndices->getLocalLength(); i++) {
252 if (myIndices->getData(0)[i] == map_->getComm()->getRank()+1) {
253 uniqueVector.push_back(map_->getGlobalElement(i));
254 }
255 }
256 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapXpetra = Xpetra::MapFactory<LO,GO,NO>::Build(map_->lib(),-1,uniqueVector(),0,map_->getComm());
257 Teuchos::RCP<Map_Xpetra<LO,GO,NO> > map = Teuchos::rcp( new Map_Xpetra<LO,GO,NO>( mapXpetra ) );
258 return map;
259 }
260
261}
262
263//merge with above function
264template <class LO,class GO,class NO>
265Teuchos::RCP<Map_Xpetra<LO,GO,NO> > Map_Xpetra<LO,GO,NO>::buildUniqueMap( tuple_intint_Type rankRange ) const
266{
267 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
268 int rank = map_->getComm()->getRank();
269 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > myIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(map_);
270 myIndices->putScalar(rank + 1);
271 GO maxGID = map_->getMaxAllGlobalIndex();
272
273 int numAvailableRanks = std::get<1>(rankRange) - std::get<0>(rankRange) + 1;
274 int numElementsForAvailRank = (int) ( ( (maxGID+1) / numAvailableRanks ) + 100.*std::numeric_limits<double>::epsilon() );
275
276 int remainingElement = maxGID+1 - numAvailableRanks * numElementsForAvailRank;
277 bool hasOneMoreElement = false;
278
279 if ( remainingElement > map_->getComm()->getRank() - std::get<0>(rankRange) ) {
280 numElementsForAvailRank++;
281 hasOneMoreElement = true;
282 }
283 if ( map_->getComm()->getRank() < std::get<0>(rankRange) || map_->getComm()->getRank() > std::get<1>(rankRange) )
284 numElementsForAvailRank = 0;
285
286
287 Teuchos::Array<GO> myElements( numElementsForAvailRank );
288 GO offset = numElementsForAvailRank * ( map_->getComm()->getRank() - std::get<0>(rankRange) );
289
290 if (!hasOneMoreElement)
291 offset += remainingElement;
292
293 for (int i=0; i<myElements.size(); i++)
294 myElements[i] = i + offset;
295
296
297 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > linearMapAvailRanks = Xpetra::MapFactory<LO,GO,NO>::Build( map_->lib(),Teuchos::OrdinalTraits<GO>::invalid(), myElements(), 0, map_->getComm() );
298
299 Teuchos::RCP<Xpetra::Vector<GO,LO,GO,NO> > globalIndices = Xpetra::VectorFactory<GO,LO,GO,NO>::Build(linearMapAvailRanks);
300
301 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer = Xpetra::ImportFactory<LO,GO,NO>::Build( map_, linearMapAvailRanks );
302 Teuchos::RCP<Xpetra::Import<LO,GO,NO> > importer2 = Xpetra::ImportFactory<LO,GO,NO>::Build( linearMapAvailRanks, map_ );
303
304 globalIndices->doImport(*myIndices,*importer,Xpetra::INSERT);
305
306 myIndices->putScalar(0);
307 myIndices->doImport(*globalIndices,*importer2,Xpetra::ADD);
308
309 Teuchos::Array<GO> uniqueVector;
310 for (unsigned i=0; i<myIndices->getLocalLength(); i++) {
311 if (myIndices->getData(0)[i] == map_->getComm()->getRank()+1) {
312 uniqueVector.push_back(map_->getGlobalElement(i));
313 }
314 }
315 Teuchos::RCP<Xpetra::Map<LO,GO,NO> > mapXpetra = Xpetra::MapFactory<LO,GO,NO>::Build(map_->lib(),-1,uniqueVector(),0,map_->getComm());
316 Teuchos::RCP<Map_Xpetra<LO,GO,NO> > map = Teuchos::rcp( new Map_Xpetra<LO,GO,NO>( mapXpetra ) );
317
318 return map;
319
320}
321
322template <class LO,class GO,class NO>
323GO Map_Xpetra<LO,GO,NO>::getIndexBase() const{
324 TEUCHOS_TEST_FOR_EXCEPTION(map_.is_null(),std::runtime_error,"map is null.");
325 return map_->getIndexBase();
326}
327}
328#endif
FEDD::Map_Xpetra::buildUniqueMap
Teuchos::RCP< Map_Xpetra< LO, GO, NO > > buildUniqueMap(int numFreeProcs=0) const
Definition Map_Xpetra_def.hpp:183
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5