Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_STK_Utilities.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 "PanzerAdaptersSTK_config.hpp"
44
45#ifdef PANZER_HAVE_EPETRA_STACK
46
49
50#include "Kokkos_DynRankView.hpp"
51
52#include <stk_mesh/base/FieldBase.hpp>
53
54namespace panzer_stk {
55
56static void gather_in_block(const std::string & blockId, const panzer::GlobalIndexer& dofMngr,
57 const Epetra_Vector & x,const std::vector<std::size_t> & localCellIds,
58 std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > & fc);
59
60static void build_local_ids(const panzer_stk::STK_Interface & mesh,
61 std::map<std::string,Teuchos::RCP<std::vector<std::size_t> > > & localIds);
62
63void write_cell_data(panzer_stk::STK_Interface & mesh,const std::vector<double> & data,const std::string & fieldName)
64{
65 std::vector<std::string> blocks;
66 mesh.getElementBlockNames(blocks);
67
68 // loop over element blocks
69 for(std::size_t eb=0;eb<blocks.size();eb++) {
70 const std::string & blockId = blocks[eb];
72
73 std::vector<stk::mesh::Entity> elements;
74 mesh.getMyElements(blockId,elements);
75
76 // loop over elements in this block
77 for(std::size_t el=0;el<elements.size();el++) {
78 std::size_t localId = mesh.elementLocalId(elements[el]);
79 double * solnData = stk::mesh::field_data(*field,elements[el]);
80 TEUCHOS_ASSERT(solnData!=0); // sanity check
81 solnData[0] = data[localId];
82 }
83 }
84}
85
86void write_solution_data(const panzer::GlobalIndexer& dofMngr,panzer_stk::STK_Interface & mesh,const Epetra_MultiVector & x,const std::string & prefix,const std::string & postfix)
87{
88 write_solution_data(dofMngr,mesh,*x(0),prefix,postfix);
89}
90
91void write_solution_data(const panzer::GlobalIndexer& dofMngr,panzer_stk::STK_Interface & mesh,const Epetra_Vector & x,const std::string & prefix,const std::string & postfix)
92{
93 typedef Kokkos::DynRankView<double,PHX::Device> FieldContainer;
94
95 // get local IDs
96 std::map<std::string,Teuchos::RCP<std::vector<std::size_t> > > localIds;
97 build_local_ids(mesh,localIds);
98
99 // loop over all element blocks
100 for(const auto & itr : localIds) {
101 const auto blockId = itr.first;
102 const auto & localCellIds = *(itr.second);
103
104 std::map<std::string,FieldContainer> data;
105
106 // get all solution data for this block
107 gather_in_block(blockId,dofMngr,x,localCellIds,data);
108
109 // write out to stk mesh
110 for(const auto & dataItr : data)
111 mesh.setSolutionFieldData(prefix+dataItr.first+postfix,blockId,localCellIds,dataItr.second);
112 }
113}
114
115void gather_in_block(const std::string & blockId, const panzer::GlobalIndexer& dofMngr,
116 const Epetra_Vector & x,const std::vector<std::size_t> & localCellIds,
117 std::map<std::string,Kokkos::DynRankView<double,PHX::Device> > & fc)
118{
119 const std::vector<int> & fieldNums = dofMngr.getBlockFieldNumbers(blockId);
120
121 for(std::size_t fieldIndex=0;fieldIndex<fieldNums.size();fieldIndex++) {
122 int fieldNum = fieldNums[fieldIndex];
123 std::string fieldStr = dofMngr.getFieldString(fieldNum);
124
125 // grab the field
126 const std::vector<int> & elmtOffset = dofMngr.getGIDFieldOffsets(blockId,fieldNum);
127 fc[fieldStr] = Kokkos::DynRankView<double,PHX::Device>("fc",localCellIds.size(),elmtOffset.size());
128 auto field = Kokkos::create_mirror_view(fc[fieldStr]);
129
130
131 // gather operation for each cell in workset
132 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
133 std::vector<panzer::GlobalOrdinal> GIDs;
134 std::vector<int> LIDs;
135 std::size_t cellLocalId = localCellIds[worksetCellIndex];
136
137 dofMngr.getElementGIDs(cellLocalId,GIDs);
138
139 // caculate the local IDs for this element
140 LIDs.resize(GIDs.size());
141 for(std::size_t i=0;i<GIDs.size();i++)
142 LIDs[i] = x.Map().LID(GIDs[i]);
143
144 // loop over basis functions and fill the fields
145 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
146 int offset = elmtOffset[basis];
147 int lid = LIDs[offset];
148 field(worksetCellIndex,basis) = x[lid];
149 }
150 }
151 Kokkos::deep_copy(fc[fieldStr], field);
152 }
153}
154
155void build_local_ids(const panzer_stk::STK_Interface & mesh,
156 std::map<std::string,Teuchos::RCP<std::vector<std::size_t> > > & localIds)
157{
158 // defines ordering of blocks
159 std::vector<std::string> blockIds;
160 mesh.getElementBlockNames(blockIds);
161
162 std::vector<std::string>::const_iterator idItr;
163 for(idItr=blockIds.begin();idItr!=blockIds.end();++idItr) {
164 std::string blockId = *idItr;
165
166 localIds[blockId] = Teuchos::rcp(new std::vector<std::size_t>);
167 std::vector<std::size_t> & localBlockIds = *localIds[blockId];
168
169 // grab elements on this block
170 std::vector<stk::mesh::Entity> blockElmts;
171 mesh.getMyElements(blockId,blockElmts);
172
173 std::vector<stk::mesh::Entity>::const_iterator itr;
174 for(itr=blockElmts.begin();itr!=blockElmts.end();++itr)
175 localBlockIds.push_back(mesh.elementLocalId(*itr));
176
177 std::sort(localBlockIds.begin(),localBlockIds.end());
178 }
179}
180
181}
182
183#endif // PANZER_HAVE_EPETRA_STACK
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
int LID(int GID) const
const Epetra_BlockMap & Map() const
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
virtual const std::vector< int > & getBlockFieldNumbers(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.
stk::mesh::Field< double > SolutionFieldType
void setSolutionFieldData(const std::string &fieldName, const std::string &blockId, const std::vector< std::size_t > &localElementIds, const ArrayT &solutionValues, double scaleValue=1.0)
std::size_t elementLocalId(stk::mesh::Entity elmt) const
void getElementBlockNames(std::vector< std::string > &names) const
stk::mesh::Field< double > * getCellField(const std::string &fieldName, const std::string &blockId) const
void getMyElements(std::vector< stk::mesh::Entity > &elements) const