Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_GlobalIndexer_EpetraUtilities_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#ifndef PANZER_UNIQUEGLOBALINDEXER_EPETRAUTILITIES_IMPL_HPP
43#define PANZER_UNIQUEGLOBALINDEXER_EPETRAUTILITIES_IMPL_HPP
44
45#include <vector>
46#include <map>
47
48#include "Teuchos_FancyOStream.hpp"
49#include "Teuchos_ArrayView.hpp"
50#include "Teuchos_CommHelpers.hpp"
51
52#include "Epetra_Map.h"
53#include "Epetra_IntVector.h"
54#include "Epetra_MultiVector.h"
55#include "Epetra_Import.h"
56#include "Epetra_MpiComm.h"
57
58#include "Epetra_CombineMode.h"
59#include "Epetra_DataAccess.h"
60
61#include <sstream>
62#include <cmath>
63
64namespace panzer {
65
66template <typename ScalarT,typename ArrayT>
67void updateGhostedDataReducedVectorEpetra(const std::string & fieldName,const std::string blockId,
68 const GlobalIndexer & ugi,
69 const ArrayT & data,Epetra_MultiVector & dataVector)
70{
71 typedef Epetra_BlockMap Map;
72
73 TEUCHOS_TEST_FOR_EXCEPTION(!ugi.fieldInBlock(fieldName,blockId),std::runtime_error,
74 "panzer::updateGhostedDataReducedVector: field name = \""+fieldName+"\" is not in element block = \"" +blockId +"\"!");
75
76 Teuchos::RCP<const Map> dataMap = Teuchos::rcpFromRef(dataVector.Map());
77
78 int fieldNum = ugi.getFieldNum(fieldName);
79 const std::vector<panzer::LocalOrdinal> & elements = ugi.getElementBlock(blockId);
80 const std::vector<int> & fieldOffsets = ugi.getGIDFieldOffsets(blockId,fieldNum);
81
82 TEUCHOS_TEST_FOR_EXCEPTION(data.extent(0)!=elements.size(),std::runtime_error,
83 "panzer::updateGhostedDataReducedVector: data cell dimension does not match up with block cell count");
84
85 int rank = data.rank();
86
87 if(rank==2) {
88 // loop over elements distributing relevent data to vector
89 std::vector<panzer::GlobalOrdinal> gids;
90 for(std::size_t e=0;e<elements.size();e++) {
91 ugi.getElementGIDs(elements[e],gids);
92
93 for(std::size_t f=0;f<fieldOffsets.size();f++) {
94 std::size_t localIndex = dataMap->LID(gids[fieldOffsets[f]]); // hash table lookup
95 dataVector.ReplaceMyValue(localIndex,0,data(e,f));
96 }
97 }
98 }
99 else if(rank==3) {
100 std::size_t entries = data.extent(2);
101
102 TEUCHOS_TEST_FOR_EXCEPTION(dataVector.NumVectors()!=static_cast<int>(entries),std::runtime_error,
103 "panzer::updateGhostedDataReducedVector: number of columns in data vector inconsistent with data array");
104
105 // loop over elements distributing relevent data to vector
106 std::vector<panzer::GlobalOrdinal> gids;
107 for(std::size_t e=0;e<elements.size();e++) {
108 ugi.getElementGIDs(elements[e],gids);
109
110 for(std::size_t f=0;f<fieldOffsets.size();f++) {
111 std::size_t localIndex = dataMap->LID(gids[fieldOffsets[f]]); // hash table lookup
112 for(std::size_t v=0;v<entries;v++)
113 dataVector.ReplaceMyValue(localIndex,v,data(e,f,v));
114 }
115 }
116 }
117 else
118 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
119 "panzer::updateGhostedDataReducedVector: data array rank must be 2 or 3");
120}
121
122template <typename ScalarT,typename ArrayT>
123Teuchos::RCP<Epetra_MultiVector>
125 const std::map<std::string,ArrayT> & data) const
126{
127 TEUCHOS_ASSERT(data.size()>0); // there must be at least one "data" item
128
129 int fieldNum = ugi_->getFieldNum(fieldName);
130 std::vector<std::string> blockIds;
131 ugi_->getElementBlockIds(blockIds);
132
133 // get rank of first data array, determine column count
134 int rank = data.begin()->second.rank();
135 int numCols = 0;
136 if(rank==2)
137 numCols = 1;
138 else if(rank==3)
139 numCols = data.begin()->second.extent(2);
140 else {
141 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
142 "ArrayToFieldVectorEpetra::getGhostedDataVector: data array must have rank 2 or 3. This array has rank " << rank << ".");
143 }
144
145
146 // first build and fill in final reduced field vector
148
149 // build field maps as needed
150 Teuchos::RCP<const Map> reducedMap = gh_reducedFieldMaps_[fieldNum];
151 if(gh_reducedFieldMaps_[fieldNum]==Teuchos::null) {
153 gh_reducedFieldMaps_[fieldNum] = reducedMap;
154 }
155
156 Teuchos::RCP<MultiVector> finalReducedVec
157 = Teuchos::rcp(new MultiVector(*reducedMap,numCols));
158 for(std::size_t b=0;b<blockIds.size();b++) {
159 std::string block = blockIds[b];
160
161 // make sure field is in correct block
162 if(!ugi_->fieldInBlock(fieldName,block))
163 continue;
164
165 // extract data vector
166 typename std::map<std::string,ArrayT>::const_iterator blockItr = data.find(block);
167 TEUCHOS_TEST_FOR_EXCEPTION(blockItr==data.end(),std::runtime_error,
168 "ArrayToFieldVectorEpetra::getDataVector: can not find block \""+block+"\".");
169
170 const ArrayT & d = blockItr->second;
171 updateGhostedDataReducedVectorEpetra<ScalarT,ArrayT>(fieldName,block,*ugi_,d,*finalReducedVec);
172 }
173
174 // build final (not reduced vector)
176
177 Teuchos::RCP<const Map> map = gh_fieldMaps_[fieldNum];
178 if(gh_fieldMaps_[fieldNum]==Teuchos::null) {
180 gh_fieldMaps_[fieldNum] = map;
181 }
182
183 Teuchos::RCP<MultiVector> finalVec
184 = Teuchos::rcp(new MultiVector(*map,numCols));
185
186 // do import from finalReducedVec
187 Epetra_Import importer(*map,*reducedMap);
188 finalVec->Import(*finalReducedVec,importer,Insert);
189
190 return finalVec;
191}
192
193template <typename ScalarT,typename ArrayT>
194Teuchos::RCP<Epetra_MultiVector>
195ArrayToFieldVectorEpetra::getDataVector(const std::string & fieldName,
196 const std::map<std::string,ArrayT> & data) const
197{
198 // if neccessary build field vector
199 if(fieldVector_==Teuchos::null)
201
202 Teuchos::RCP<const MultiVector> sourceVec
203 = getGhostedDataVector<ScalarT,ArrayT>(fieldName,data);
204
205 // use lazy construction for each field
206 int fieldNum = ugi_->getFieldNum(fieldName);
207 Teuchos::RCP<const Map> destMap = fieldMaps_[fieldNum];
208 if(fieldMaps_[fieldNum]==Teuchos::null) {
209 destMap = panzer::getFieldMapEpetra(fieldNum,*fieldVector_);
210 fieldMaps_[fieldNum] = destMap;
211 }
212
213 Teuchos::RCP<MultiVector> destVec
214 = Teuchos::rcp(new MultiVector(*destMap,sourceVec->NumVectors()));
215
216 // do import
217 Epetra_Import importer(*destMap,sourceVec->Map());
218 destVec->Import(*sourceVec,importer,Insert);
219
220 return destVec;
221}
222
223} // end namspace panzer
224
225#endif
Insert
const Epetra_BlockMap & Map() const
int NumVectors() const
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Teuchos::RCP< const GlobalIndexer > ugi_
DOF mapping.
std::map< int, Teuchos::RCP< const Map > > gh_reducedFieldMaps_
ghosted field vector
std::map< int, Teuchos::RCP< const Map > > gh_fieldMaps_
Maps for each field (as needed)
void buildFieldVector(const Epetra_IntVector &source) const
build unghosted field vector from ghosted field vector
Teuchos::RCP< const IntVector > gh_fieldVector_
ghosted reduced field vector
Teuchos::RCP< Epetra_MultiVector > getGhostedDataVector(const std::string &fieldName, const std::map< std::string, ArrayT > &data) const
Teuchos::RCP< const IntVector > fieldVector_
Maps for each field (as needed)
std::map< int, Teuchos::RCP< const Map > > fieldMaps_
(unghosted) field vector (as needed)
Teuchos::RCP< Epetra_MultiVector > getDataVector(const std::string &fieldName, const std::map< std::string, ArrayT > &data) const
virtual const std::vector< panzer::LocalOrdinal > & getElementBlock(const std::string &blockId) const =0
virtual const std::vector< int > & getGIDFieldOffsets(const std::string &blockId, int fieldNum) const =0
Use the field pattern so that you can find a particular field in the GIDs array.
virtual void getElementGIDs(panzer::LocalOrdinal localElmtId, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
virtual bool fieldInBlock(const std::string &field, const std::string &block) const =0
void updateGhostedDataReducedVectorEpetra(const std::string &fieldName, const std::string blockId, const GlobalIndexer &ugi, const ArrayT &data, Epetra_MultiVector &dataVector)
Teuchos::RCP< const Epetra_BlockMap > getFieldMapEpetra(int fieldNum, const Epetra_IntVector &fieldVector)