Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_GlobalIndexer_Utilities_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 #include <vector>
44 #include <map>
45 
46 #include "Teuchos_FancyOStream.hpp"
47 #include "Teuchos_ArrayView.hpp"
48 #include "Teuchos_CommHelpers.hpp"
49 
50 #include "Tpetra_Map.hpp"
51 #include "Tpetra_Vector.hpp"
52 #include "Tpetra_Import.hpp"
53 
54 #include "Kokkos_DynRankView.hpp"
55 
56 #include <sstream>
57 #include <cmath>
58 
59 namespace panzer {
60 
61 template <typename ScalarT,typename ArrayT>
62 void updateGhostedDataReducedVector(const std::string & fieldName,const std::string blockId,
63  const GlobalIndexer & ugi,
64  const ArrayT & data,Tpetra::MultiVector<ScalarT,int,panzer::GlobalOrdinal,panzer::TpetraNodeType> & dataVector)
65 {
66  typedef Tpetra::Map<int,panzer::GlobalOrdinal,panzer::TpetraNodeType> Map;
67 
68  TEUCHOS_TEST_FOR_EXCEPTION(!ugi.fieldInBlock(fieldName,blockId),std::runtime_error,
69  "panzer::updateGhostedDataReducedVector: field name = \""+fieldName+"\" is not in element block = \"" +blockId +"\"!");
70 
71  Teuchos::RCP<const Map> dataMap = dataVector.getMap();
72 
73  int fieldNum = ugi.getFieldNum(fieldName);
74  const std::vector<panzer::LocalOrdinal> & elements = ugi.getElementBlock(blockId);
75  const std::vector<int> & fieldOffsets = ugi.getGIDFieldOffsets(blockId,fieldNum);
76 
77  TEUCHOS_TEST_FOR_EXCEPTION(data.extent(0)!=elements.size(),std::runtime_error,
78  "panzer::updateGhostedDataReducedVector: data cell dimension does not match up with block cell count");
79 
80  int rank = data.rank();
81  auto dataVector_host = dataVector.getLocalViewHost(Tpetra::Access::ReadWrite);
82  auto data_host = Kokkos::create_mirror_view(data);
83  Kokkos::deep_copy(data_host,data);
84 
85  if(rank==2) {
86  // loop over elements distributing relevent data to vector
87  std::vector<panzer::GlobalOrdinal> gids;
88  for(std::size_t e=0;e<elements.size();e++) {
89  ugi.getElementGIDs(elements[e],gids);
90 
91  for(std::size_t f=0;f<fieldOffsets.size();f++) {
92  std::size_t localIndex = dataMap->getLocalElement(gids[fieldOffsets[f]]); // hash table lookup
93  dataVector_host(localIndex,0) = data_host(e,f);
94  }
95  }
96  }
97  else if(rank==3) {
98  std::size_t entries = data.extent(2);
99 
100  TEUCHOS_TEST_FOR_EXCEPTION(dataVector.getNumVectors()!=entries,std::runtime_error,
101  "panzer::updateGhostedDataReducedVector: number of columns in data vector inconsistent with data array");
102 
103  // loop over elements distributing relevent data to vector
104  std::vector<panzer::GlobalOrdinal> gids;
105  for(std::size_t e=0;e<elements.size();e++) {
106  ugi.getElementGIDs(elements[e],gids);
107 
108  for(std::size_t f=0;f<fieldOffsets.size();f++) {
109  std::size_t localIndex = dataMap->getLocalElement(gids[fieldOffsets[f]]); // hash table lookup
110  for(std::size_t v=0;v<entries;v++)
111  dataVector_host(localIndex,v) = data_host(e,f,v);
112  }
113  }
114  }
115  else
116  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
117  "panzer::updateGhostedDataReducedVector: data array rank must be 2 or 3");
118 }
119 
120 template <typename ScalarT,typename ArrayT>
122 ArrayToFieldVector::getGhostedDataVector(const std::string & fieldName,const std::map<std::string,ArrayT> & data) const
123 {
124  TEUCHOS_ASSERT(data.size()>0); // there must be at least one "data" item
125 
126  int fieldNum = ugi_->getFieldNum(fieldName);
127  std::vector<std::string> blockIds;
128  ugi_->getElementBlockIds(blockIds);
129 
130  // get rank of first data array, determine column count
131  int rank = data.begin()->second.rank();
132  int numCols = 0;
133  if(rank==2)
134  numCols = 1;
135  else if(rank==3)
136  numCols = data.begin()->second.extent(2);
137  else {
138  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
139  "ArrayToFieldVector::getGhostedDataVector: data array must have rank 2 or 3. This array has rank " << rank << ".");
140  }
141 
142 
143  // first build and fill in final reduced field vector
145 
146  // build field maps as needed
147  Teuchos::RCP<const Map> reducedMap = gh_reducedFieldMaps_[fieldNum];
148  if(gh_reducedFieldMaps_[fieldNum]==Teuchos::null) {
149  reducedMap = panzer::getFieldMap(fieldNum,*gh_reducedFieldVector_);
150  gh_reducedFieldMaps_[fieldNum] = reducedMap;
151  }
152 
154  = Teuchos::rcp(new Tpetra::MultiVector<ScalarT,int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(reducedMap,numCols));
155  for(std::size_t b=0;b<blockIds.size();b++) {
156  std::string block = blockIds[b];
157 
158  // make sure field is in correct block
159  if(!ugi_->fieldInBlock(fieldName,block))
160  continue;
161 
162  // extract data vector
163  typename std::map<std::string,ArrayT>::const_iterator blockItr = data.find(block);
164  TEUCHOS_TEST_FOR_EXCEPTION(blockItr==data.end(),std::runtime_error,
165  "ArrayToFieldVector::getDataVector: can not find block \""+block+"\".");
166 
167  const ArrayT & d = blockItr->second;
168  updateGhostedDataReducedVector<ScalarT,ArrayT>(fieldName,block,*ugi_,d,*finalReducedVec);
169  }
170 
171  // build final (not reduced vector)
173 
174  Teuchos::RCP<const Map> map = gh_fieldMaps_[fieldNum];
175  if(gh_fieldMaps_[fieldNum]==Teuchos::null) {
176  map = panzer::getFieldMap(fieldNum,*gh_fieldVector_);
177  gh_fieldMaps_[fieldNum] = map;
178  }
179 
181  = Teuchos::rcp(new Tpetra::MultiVector<ScalarT,int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(map,numCols));
182 
183  // do import from finalReducedVec
184  Tpetra::Import<int,panzer::GlobalOrdinal,panzer::TpetraNodeType> importer(reducedMap,map);
185  finalVec->doImport(*finalReducedVec,importer,Tpetra::INSERT);
186  PHX::Device::execution_space().fence();
187 
188  return finalVec;
189 }
190 
191 template <typename ScalarT,typename ArrayT>
193 ArrayToFieldVector::getDataVector(const std::string & fieldName,const std::map<std::string,ArrayT> & data) const
194 {
195  // if neccessary build field vector
196  if(fieldVector_==Teuchos::null)
198 
200  = getGhostedDataVector<ScalarT,ArrayT>(fieldName,data);
201 
202  // use lazy construction for each field
203  int fieldNum = ugi_->getFieldNum(fieldName);
204  Teuchos::RCP<const Map> destMap = fieldMaps_[fieldNum];
205  if(fieldMaps_[fieldNum]==Teuchos::null) {
206  destMap = panzer::getFieldMap(fieldNum,*fieldVector_);
207  fieldMaps_[fieldNum] = destMap;
208  }
209 
211  = Teuchos::rcp(new Tpetra::MultiVector<ScalarT,int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(destMap,sourceVec->getNumVectors()));
212 
213  // do import
214  Tpetra::Import<int,panzer::GlobalOrdinal,panzer::TpetraNodeType> importer(sourceVec->getMap(),destMap);
215  destVec->doImport(*sourceVec,importer,Tpetra::INSERT);
216  PHX::Device::execution_space().fence();
217 
218  return destVec;
219 }
220 
221 } // end namspace panzer
std::map< int, Teuchos::RCP< const Map > > fieldMaps_
(unghosted) field vector (as needed)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::map< int, Teuchos::RCP< const Map > > gh_fieldMaps_
Maps for each field (as needed)
Teuchos::RCP< const IntVector > gh_reducedFieldVector_
virtual const std::vector< panzer::LocalOrdinal > & getElementBlock(const std::string &blockId) const =0
std::map< int, Teuchos::RCP< const Map > > gh_reducedFieldMaps_
ghosted field vector
Teuchos::RCP< const IntVector > gh_fieldVector_
ghosted reduced field vector
Teuchos::RCP< const Tpetra::Map< int, panzer::GlobalOrdinal, panzer::TpetraNodeType > > getFieldMap(int fieldNum, const Tpetra::Vector< int, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > &fieldTVector)
virtual bool fieldInBlock(const std::string &field, const std::string &block) const =0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void updateGhostedDataReducedVector(const std::string &fieldName, const std::string blockId, const GlobalIndexer &ugi, const ArrayT &data, Tpetra::MultiVector< ScalarT, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > &dataVector)
Teuchos::RCP< const GlobalIndexer > ugi_
DOF mapping.
Teuchos::RCP< Tpetra::MultiVector< ScalarT, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > > getGhostedDataVector(const std::string &fieldName, const std::map< std::string, ArrayT > &data) const
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
Teuchos::RCP< Tpetra::MultiVector< ScalarT, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > > getDataVector(const std::string &fieldName, const std::map< std::string, ArrayT > &data) const
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.
void buildFieldVector(const Tpetra::Vector< int, int, panzer::GlobalOrdinal, panzer::TpetraNodeType > &source) const
build unghosted field vector from ghosted field vector
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< const IntVector > fieldVector_
Maps for each field (as needed)