Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_GlobalIndexer_EpetraUtilities.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
43#include "PanzerDofMgr_config.hpp"
44
45#ifdef PANZER_HAVE_EPETRA
46
48
49namespace panzer {
50
51Teuchos::RCP<Epetra_IntVector>
52buildGhostedFieldReducedVectorEpetra(const GlobalIndexer & ugi)
53{
54 typedef Epetra_BlockMap Map;
55 typedef Epetra_IntVector IntVector;
56
57 std::vector<panzer::GlobalOrdinal> indices;
58 std::vector<std::string> blocks;
59
60 ugi.getOwnedAndGhostedIndices(indices);
61 ugi.getElementBlockIds(blocks);
62
63 std::vector<int> fieldNumbers(indices.size(),-1);
64
65 const Teuchos::RCP<const Teuchos::MpiComm<int> > mpiComm = Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(ugi.getComm());
66 Teuchos::RCP<Epetra_MpiComm> comm;
67 if (mpiComm != Teuchos::null)
68 comm = Teuchos::rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()));
69
70 Teuchos::RCP<Map> ghostedMap
71 = Teuchos::rcp(new Map(-1, static_cast<int>(indices.size()), Teuchos::arrayViewFromVector(indices).getRawPtr(),
72 1, Teuchos::OrdinalTraits<int>::zero(), *comm));
73
74 // build a map from local ids to a field number
75 for(std::size_t blk=0;blk<blocks.size();blk++) {
76 std::string blockId = blocks[blk];
77
78 const std::vector<panzer::LocalOrdinal> & elements = ugi.getElementBlock(blockId);
79 const std::vector<int> & fields = ugi.getBlockFieldNumbers(blockId);
80
81 // loop over all elements, and set field number in output array
82 std::vector<panzer::GlobalOrdinal> gids(fields.size());
83 for(std::size_t e=0;e<elements.size();e++) {
84 ugi.getElementGIDs(elements[e],gids);
85
86 for(std::size_t f=0;f<fields.size();f++) {
87 int fieldNum = fields[f];
88 panzer::GlobalOrdinal gid = gids[f];
89 std::size_t lid = ghostedMap->LID(gid); // hash table lookup
90
91 fieldNumbers[lid] = fieldNum;
92 }
93 }
94 }
95
96 // produce a reduced vector containing only fields known by this processor
97 std::vector<panzer::GlobalOrdinal> reducedIndices;
98 std::vector<int> reducedFieldNumbers;
99 for(std::size_t i=0;i<fieldNumbers.size();i++) {
100 if(fieldNumbers[i]>-1) {
101 reducedIndices.push_back(indices[i]);
102 reducedFieldNumbers.push_back(fieldNumbers[i]);
103 }
104 }
105
106 Teuchos::RCP<Map> reducedMap
107 = Teuchos::rcp(new Map(-1, static_cast<int>(reducedIndices.size()), Teuchos::arrayViewFromVector(reducedIndices).getRawPtr(),
108 1, Teuchos::OrdinalTraits<int>::zero(), *comm));
109 return Teuchos::rcp(new IntVector(Copy,*reducedMap,Teuchos::arrayViewFromVector(reducedFieldNumbers).getRawPtr()));
110}
111
112void buildGhostedFieldVectorEpetra(const GlobalIndexer & ugi,
113 std::vector<int> & fieldNumbers,
114 const Teuchos::RCP<const Epetra_IntVector> & reducedVec)
115{
116 typedef Epetra_IntVector IntVector;
117
118 Teuchos::RCP<const IntVector> dest = buildGhostedFieldVectorEpetra(ugi,reducedVec);
119
120 fieldNumbers.resize(dest->MyLength());
121 Teuchos::ArrayView<int> av = Teuchos::arrayViewFromVector(fieldNumbers);
122 dest->ExtractCopy(av.getRawPtr());
123}
124
125Teuchos::RCP<const Epetra_IntVector>
126buildGhostedFieldVectorEpetra(const GlobalIndexer & ugi,
127 const Teuchos::RCP<const Epetra_IntVector> & reducedVec)
128{
129 typedef Epetra_BlockMap Map;
130 typedef Epetra_IntVector IntVector;
131 typedef Epetra_Import Importer;
132
133 // first step: get a reduced field number vector and build a map to
134 // contain the full field number vector
136
137 Teuchos::RCP<Map> destMap;
138 {
139 std::vector<panzer::GlobalOrdinal> indices;
140 ugi.getOwnedAndGhostedIndices(indices);
141
142 const Teuchos::RCP<const Teuchos::MpiComm<int> > mpiComm = Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(ugi.getComm());
143 Teuchos::RCP<Epetra_MpiComm> comm;
144 if (mpiComm != Teuchos::null)
145 comm = Teuchos::rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()));
146
147 destMap = Teuchos::rcp(new Map(-1, static_cast<int>(indices.size()), Teuchos::arrayViewFromVector(indices).getRawPtr(),
148 1, Teuchos::OrdinalTraits<int>::zero(), *comm));
149 }
150
151 Teuchos::RCP<const IntVector> source = reducedVec;
152 if(source==Teuchos::null)
154 Teuchos::RCP<const Map> sourceMap = Teuchos::rcpFromRef(source->Map());
155
156 // second step: perform the global communciation required to fix the
157 // interface conditions (where this processor doesn't know what field
158 // some indices are)
160 Teuchos::RCP<IntVector> dest = Teuchos::rcp(new IntVector(*destMap));
161 Importer importer(*destMap,*sourceMap);
162
163 dest->Import(*source,importer,Insert);
164
165 return dest;
166}
167
168ArrayToFieldVectorEpetra::ArrayToFieldVectorEpetra(const Teuchos::RCP<const GlobalIndexer> & ugi)
169 : ugi_(ugi)
170{
171 gh_reducedFieldVector_ = buildGhostedFieldReducedVectorEpetra(*ugi_);
172 gh_fieldVector_ = buildGhostedFieldVectorEpetra(*ugi_,gh_reducedFieldVector_);
173}
174
175
176void ArrayToFieldVectorEpetra::buildFieldVector(const Epetra_IntVector & source) const
177{
178 // build (unghosted) vector and map
179 std::vector<panzer::GlobalOrdinal> indices;
180 ugi_->getOwnedIndices(indices);
181
182 const Teuchos::RCP<const Teuchos::MpiComm<int> > mpiComm = Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(ugi_->getComm(),true);
183 Teuchos::RCP<Epetra_MpiComm> comm;
184 if (mpiComm != Teuchos::null)
185 comm = Teuchos::rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()));
186
187 Teuchos::RCP<Map> destMap
188 = Teuchos::rcp(new Map(-1, static_cast<int>(indices.size()),
189 Teuchos::arrayViewFromVector(indices).getRawPtr(),
190 // &indices.begin(),
191 1, Teuchos::OrdinalTraits<int>::zero(), *comm));
192
193 Teuchos::RCP<IntVector> localFieldVector = Teuchos::rcp(new IntVector(*destMap));
194
195 Epetra_Import importer(*destMap,source.Map());
196 localFieldVector->Import(source,importer,Insert);
197
198 fieldVector_ = localFieldVector;
199}
200
201Teuchos::RCP<const Epetra_Map>
202ArrayToFieldVectorEpetra::getFieldMap(const std::string & fieldName) const
203{
204 return getFieldMap(ugi_->getFieldNum(fieldName));
205}
206
207Teuchos::RCP<const Epetra_Map>
208ArrayToFieldVectorEpetra::getFieldMap(int fieldNum) const
209{
210 if(fieldMaps_[fieldNum]==Teuchos::null) {
211 // if neccessary build field vector
212 if(fieldVector_==Teuchos::null)
213 buildFieldVector(*gh_fieldVector_);
214
215 fieldMaps_[fieldNum] = panzer::getFieldMapEpetra(fieldNum,*fieldVector_);
216 }
217
218 return Teuchos::rcp_dynamic_cast<const Epetra_Map>(fieldMaps_[fieldNum]);
219}
220
221Teuchos::RCP<const Epetra_BlockMap>
222getFieldMapEpetra(int fieldNum,const Epetra_IntVector & fieldTVector)
223{
224 Teuchos::RCP<const Epetra_BlockMap> origMap = Teuchos::rcpFromRef(fieldTVector.Map());
225 std::vector<int> fieldVector(fieldTVector.MyLength());
226 Teuchos::ArrayView<int> av = Teuchos::arrayViewFromVector(fieldVector);
227 fieldTVector.ExtractCopy(av.getRawPtr());
228
229 std::vector<int> mapVector;
230 for(std::size_t i=0;i<fieldVector.size();i++) {
231 if(fieldVector[i]==fieldNum)
232 mapVector.push_back(origMap->GID(i));
233 }
234
235 Teuchos::RCP<Epetra_BlockMap> finalMap
236 = Teuchos::rcp(new Epetra_BlockMap(-1, static_cast<int>(mapVector.size()), Teuchos::arrayViewFromVector(mapVector).getRawPtr(), 1, Teuchos::OrdinalTraits<int>::zero(), origMap->Comm()));
237
238 return finalMap;
239}
240
241} // end namespace panzer
242
243#endif //end PANZER_HAVE_EPETRA
const Epetra_BlockMap & Map() const
int MyLength() const
int ExtractCopy(int *V) const
ArrayToFieldVectorEpetra()
Maps for each field (as needed)
Teuchos::RCP< const Tpetra::Map< int, panzer::GlobalOrdinal, panzer::TpetraNodeType > > getFieldMap(int fieldNum, const Tpetra::Vector< int, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > &fieldTVector)
Teuchos::RCP< const Epetra_IntVector > buildGhostedFieldVectorEpetra(const GlobalIndexer &ugi, const Teuchos::RCP< const Epetra_IntVector > &reducedVec=Teuchos::null)
Teuchos::RCP< const Epetra_BlockMap > getFieldMapEpetra(int fieldNum, const Epetra_IntVector &fieldVector)
Teuchos::RCP< Epetra_IntVector > buildGhostedFieldReducedVectorEpetra(const GlobalIndexer &ugi)