Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_IntrepidFieldPattern.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
44
45#include "Teuchos_Assert.hpp"
46#include "Intrepid2_CellTools.hpp"
47#include "Shards_CellTopology.hpp"
48
49namespace panzer {
50
52 Intrepid2FieldPattern(const Teuchos::RCP< Intrepid2::Basis<PHX::Device,double,double> > &intrepidBasis)
53 : intrepidBasis_(intrepidBasis) {
54 const auto dofOrd = intrepidBasis_->getAllDofOrdinal(); // rank 3 view
55 const auto dofTag = intrepidBasis_->getAllDofTags(); // rank 2 view
56
57 const int
58 iend = dofOrd.extent(0),
59 jend = dofOrd.extent(1);
60
61 subcellIndicies_.resize(iend);
62 for (int i=0;i<iend;++i) {
63 subcellIndicies_[i].resize(jend);
64 for (int j=0;j<jend;++j) {
65 const int ord = dofOrd(i, j, 0);
66 if (ord >= 0) {
67 const int ndofs = dofTag(ord, 3);
68 subcellIndicies_[i][j].resize(ndofs);
69 for (int k=0;k<ndofs;++k)
70 subcellIndicies_[i][j][k] = dofOrd(i, j, k);
71 } else {
72 // if ordinal does not exist empty container.
73 subcellIndicies_[i][j].clear();
74 }
75 }
76 }
77 }
78
79 int
81 getSubcellCount(int dim) const
82 {
83 const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
84 return ct.getSubcellCount(dim);
85 }
86
87 const std::vector<int> &
88 Intrepid2FieldPattern::getSubcellIndices(int dim, int cellIndex) const
89 {
90 if ((dim < static_cast<int>(subcellIndicies_.size() )) and
91 (cellIndex < static_cast<int>(subcellIndicies_[dim].size())))
92 return subcellIndicies_[dim][cellIndex];
93
94 return empty_;
95 }
96
97 void
99 getSubcellClosureIndices(int dim,int cellIndex,std::vector<int> & indices) const
100 {
101 // wipe out previous values
102 indices.clear();
103
104 if (dim == 0) {
105 indices = getSubcellIndices(dim,cellIndex);
106 } else {
107 // construct full topology
108 const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
109
110 std::set<std::pair<unsigned,unsigned> > closure;
111 Intrepid2FieldPattern::buildSubcellClosure(ct,dim,cellIndex,closure);
112
113 // grab basis indices on the closure of the sub cell
114 std::set<std::pair<unsigned,unsigned> >::const_iterator itr;
115 for (itr=closure.begin();itr!=closure.end();++itr) {
116 // grab indices for this sub cell
117 const std::vector<int> & subcellIndices = getSubcellIndices(itr->first,itr->second);
118
119 // concatenate the two vectors
120 indices.insert(indices.end(),subcellIndices.begin(),subcellIndices.end());
121 }
122 }
123 }
124
125 int
127 getDimension() const
128 {
129 return intrepidBasis_->getBaseCellTopology().getDimension();
130 }
131
132 shards::CellTopology
134 getCellTopology() const
135 {
136 return intrepidBasis_->getBaseCellTopology();
137 }
138
139 void
141 getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
142 std::vector<unsigned> & nodes)
143 {
144 if(dim==0) {
145 nodes.push_back(subCell);
146 return;
147 }
148
149 // get all nodes on requested sub cell
150 unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
151 for(unsigned node=0;node<subCellNodeCount;++node)
152 nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));
153
154 // sort them so they are ordered correctly for "includes" call
155 std::sort(nodes.begin(),nodes.end());
156 }
157
158 void
160 findContainedSubcells(const shards::CellTopology & cellTopo,unsigned dim,
161 const std::vector<unsigned> & nodes,
162 std::set<std::pair<unsigned,unsigned> > & subCells)
163 {
164 unsigned subCellCount = cellTopo.getSubcellCount(dim);
165 for(unsigned subCellOrd=0;subCellOrd<subCellCount;++subCellOrd) {
166 // get all nodes in sub cell
167 std::vector<unsigned> subCellNodes;
168 getSubcellNodes(cellTopo,dim,subCellOrd,subCellNodes);
169
170 // if subCellNodes \subset nodes => add (dim,subCellOrd) to subCells
171 bool isSubset = std::includes( nodes.begin(), nodes.end(),
172 subCellNodes.begin(), subCellNodes.end());
173 if(isSubset)
174 subCells.insert(std::make_pair(dim,subCellOrd));
175
176 }
177
178 // stop recursion base case
179 if (dim==0) return;
180
181 // find subcells in next sub dimension
182 findContainedSubcells(cellTopo,dim-1,nodes,subCells);
183 }
184
185 void
187 buildSubcellClosure(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
188 std::set<std::pair<unsigned,unsigned> > & closure)
189 {
190 switch(dim) {
191 case 0:
192 closure.insert(std::make_pair(0,subCell));
193 break;
194 case 1:
195 closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,0)));
196 closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,1)));
197 closure.insert(std::make_pair(1,subCell));
198 break;
199 case 2:
200 {
201 unsigned cnt = (shards::CellTopology(cellTopo.getCellTopologyData(dim,subCell))).getSubcellCount(dim-1);
202 for(unsigned i=0;i<cnt;i++) {
203 int edge = mapCellFaceEdge(cellTopo.getCellTopologyData(),subCell,i);
204 buildSubcellClosure(cellTopo,dim-1,edge,closure);
205 }
206 closure.insert(std::make_pair(2,subCell));
207 }
208 break;
209 default:
210 // beyond a two dimension surface this thing crashes!
211 TEUCHOS_ASSERT(false);
212 };
213 }
214
215 bool
218 {
219 // we no longer use CoordsInterface
220 return true;
221 }
222
228 void
230 getInterpolatoryCoordinates(Kokkos::DynRankView<double,PHX::Device> & coords) const
231 {
232 // this may not be efficient if coords is allocated every time this function is called
233 coords = Kokkos::DynRankView<double,PHX::Device>("coords",intrepidBasis_->getCardinality(),getDimension());
234 intrepidBasis_->getDofCoords(coords);
235 }
236
242 void
244 getInterpolatoryCoordinates(const Kokkos::DynRankView<double,PHX::Device> & cellVertices,
245 Kokkos::DynRankView<double,PHX::Device> & coords) const
246 {
247 TEUCHOS_ASSERT(cellVertices.rank()==3);
248
249 int numCells = cellVertices.extent(0);
250
251 // grab the local coordinates
252 Kokkos::DynRankView<double,PHX::Device> localCoords;
253 getInterpolatoryCoordinates(localCoords);
254
255 // resize the coordinates field container
256 coords = Kokkos::DynRankView<double,PHX::Device>("coords",numCells,localCoords.extent(0),getDimension());
257
258 if(numCells>0) {
259 Intrepid2::CellTools<PHX::Device> cellTools;
260 cellTools.mapToPhysicalFrame(coords,localCoords,cellVertices,intrepidBasis_->getBaseCellTopology());
261 }
262 }
263
264 Teuchos::RCP< Intrepid2::Basis<PHX::Device,double,double> >
267
268}
Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > intrepidBasis_
virtual void getSubcellClosureIndices(int dim, int cellIndex, std::vector< int > &indices) const
static void getSubcellNodes(const shards::CellTopology &cellTopo, unsigned dim, unsigned subCell, std::vector< unsigned > &nodes)
void getInterpolatoryCoordinates(Kokkos::DynRankView< double, PHX::Device > &coords) const
virtual shards::CellTopology getCellTopology() const
Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > getIntrepidBasis() const
Returns the underlying Intrepid2::Basis object.
bool supportsInterpolatoryCoordinates() const
Does this field pattern support interpolatory coordinates?
Intrepid2FieldPattern(const Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > &intrepidBasis)
std::vector< std::vector< std::vector< int > > > subcellIndicies_
static void findContainedSubcells(const shards::CellTopology &cellTopo, unsigned dim, const std::vector< unsigned > &nodes, std::set< std::pair< unsigned, unsigned > > &subCells)
virtual const std::vector< int > & getSubcellIndices(int dim, int cellIndex) const
static void buildSubcellClosure(const shards::CellTopology &cellTopo, unsigned dim, unsigned subCell, std::set< std::pair< unsigned, unsigned > > &closure)