Panzer  Version of the Day
Panzer_IntrepidFieldPattern_DynRankView.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 #include "Intrepid2_config.h"
43 #if defined( HAVE_INTREPID2_KOKKOS_DYNRANKVIEW )
45 
46 #include "Teuchos_Assert.hpp"
47 #include "Intrepid2_CellTools.hpp"
48 #include "Shards_CellTopology.hpp"
49 
50 namespace panzer {
51 
52  int
54  getSubcellCount(int dim) const
55  {
56  const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
57  return ct.getSubcellCount(dim);
58  }
59 
60  const std::vector<int> &
61  Intrepid2FieldPattern::getSubcellIndices(int dim, int cellIndex) const
62  {
63  subcellIndices_.clear();
64 
65  const int ord = intrepidBasis_->getDofOrdinal(dim, cellIndex, 0);
66  if (ord >= 0) {
67  const int ndof = intrepidBasis_->getDofTag(ord)(3);
68  const auto tag = intrepidBasis_->getAllDofOrdinal();
69 
70  for (int i=0;i<ndof;++i)
71  subcellIndices_.push_back(tag(dim, cellIndex, i));
72  }
73  return subcellIndices_;
74  }
75 
76  void
78  getSubcellClosureIndices(int dim,int cellIndex,std::vector<int> & indices) const
79  {
80  // wipe out previous values
81  indices.clear();
82 
83  if (dim == 0) {
84  indices = getSubcellIndices(dim,cellIndex);
85  } else {
86  // construct full topology
87  const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
88 
89  std::set<std::pair<unsigned,unsigned> > closure;
90  Intrepid2FieldPattern::buildSubcellClosure(ct,dim,cellIndex,closure);
91 
92  // grab basis indices on the closure of the sub cell
93  std::set<std::pair<unsigned,unsigned> >::const_iterator itr;
94  for (itr=closure.begin();itr!=closure.end();++itr) {
95  // grab indices for this sub cell
96  const std::vector<int> & subcellIndices = getSubcellIndices(itr->first,itr->second);
97 
98  // concatenate the two vectors
99  indices.insert(indices.end(),subcellIndices.begin(),subcellIndices.end());
100  }
101  }
102  }
103 
104  int
106  getDimension() const
107  {
108  return intrepidBasis_->getBaseCellTopology().getDimension();
109  }
110 
111  shards::CellTopology
113  getCellTopology() const
114  {
115  return intrepidBasis_->getBaseCellTopology();
116  }
117 
118  void
120  getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
121  std::vector<unsigned> & nodes)
122  {
123  if(dim==0) {
124  nodes.push_back(subCell);
125  return;
126  }
127 
128  // get all nodes on requested sub cell
129  unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
130  for(unsigned node=0;node<subCellNodeCount;++node)
131  nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));
132 
133  // sort them so they are ordered correctly for "includes" call
134  std::sort(nodes.begin(),nodes.end());
135  }
136 
137  void
139  findContainedSubcells(const shards::CellTopology & cellTopo,unsigned dim,
140  const std::vector<unsigned> & nodes,
141  std::set<std::pair<unsigned,unsigned> > & subCells)
142  {
143  unsigned subCellCount = cellTopo.getSubcellCount(dim);
144  for(unsigned subCellOrd=0;subCellOrd<subCellCount;++subCellOrd) {
145  // get all nodes in sub cell
146  std::vector<unsigned> subCellNodes;
147  getSubcellNodes(cellTopo,dim,subCellOrd,subCellNodes);
148 
149  // if subCellNodes \subset nodes => add (dim,subCellOrd) to subCells
150  bool isSubset = std::includes( nodes.begin(), nodes.end(),
151  subCellNodes.begin(), subCellNodes.end());
152  if(isSubset)
153  subCells.insert(std::make_pair(dim,subCellOrd));
154 
155  }
156 
157  // stop recursion base case
158  if (dim==0) return;
159 
160  // find subcells in next sub dimension
161  findContainedSubcells(cellTopo,dim-1,nodes,subCells);
162  }
163 
164  void
166  buildSubcellClosure(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
167  std::set<std::pair<unsigned,unsigned> > & closure)
168  {
169  switch(dim) {
170  case 0:
171  closure.insert(std::make_pair(0,subCell));
172  break;
173  case 1:
174  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,0)));
175  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,1)));
176  closure.insert(std::make_pair(1,subCell));
177  break;
178  case 2:
179  {
180  unsigned cnt = (shards::CellTopology(cellTopo.getCellTopologyData(dim,subCell))).getSubcellCount(dim-1);
181  for(unsigned i=0;i<cnt;i++) {
182  int edge = mapCellFaceEdge(cellTopo.getCellTopologyData(),subCell,i);
183  buildSubcellClosure(cellTopo,dim-1,edge,closure);
184  }
185  closure.insert(std::make_pair(2,subCell));
186  }
187  break;
188  default:
189  // beyond a two dimension surface this thing crashes!
190  TEUCHOS_ASSERT(false);
191  };
192  }
193 
194  bool
197  {
198  // we no longer use CoordsInterface
199  return true;
200  }
201 
207  void
209  getInterpolatoryCoordinates(Kokkos::DynRankView<double,PHX::Device> & coords) const
210  {
211  // this may not be efficient if coords is allocated every time this function is called
212  coords = Kokkos::DynRankView<double,PHX::Device>("coords",intrepidBasis_->getCardinality(),getDimension());
213  intrepidBasis_->getDofCoords(coords);
214  }
215 
221  void
223  getInterpolatoryCoordinates(const Kokkos::DynRankView<double,PHX::Device> & cellVertices,
224  Kokkos::DynRankView<double,PHX::Device> & coords) const
225  {
226  TEUCHOS_ASSERT(cellVertices.rank()==3);
227 
228  int numCells = cellVertices.dimension(0);
229 
230  // grab the local coordinates
231  Kokkos::DynRankView<double,PHX::Device> localCoords;
232  getInterpolatoryCoordinates(localCoords);
233 
234  // resize the coordinates field container
235  coords = Kokkos::DynRankView<double,PHX::Device>("coords",numCells,localCoords.dimension(0),getDimension());
236 
237  if(numCells>0) {
238  Intrepid2::CellTools<PHX::Device> cellTools;
239  cellTools.mapToPhysicalFrame(coords,localCoords,cellVertices,intrepidBasis_->getBaseCellTopology());
240  }
241  }
242 
243 }
244 #endif
virtual int getSubcellCount(int dim) const
virtual const std::vector< int > & getSubcellIndices(int dim, int cellIndex) const
bool supportsInterpolatoryCoordinates() const
Does this field pattern support interpolatory coordinates?
static void buildSubcellClosure(const shards::CellTopology &cellTopo, unsigned dim, unsigned subCell, std::set< std::pair< unsigned, unsigned > > &closure)
Teuchos::RCP< Intrepid2::Basis< double, Kokkos::DynRankView< double, PHX::Device > > > 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
static void findContainedSubcells(const shards::CellTopology &cellTopo, unsigned dim, const std::vector< unsigned > &nodes, std::set< std::pair< unsigned, unsigned > > &subCells)
#define TEUCHOS_ASSERT(assertion_test)