Panzer  Version of the Day
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 #include "Intrepid2_config.h"
43 #if ! defined( HAVE_INTREPID2_KOKKOS_DYNRANKVIEW )
44 
46 
47 #include "Intrepid2_CellTools.hpp"
48 #include "Shards_CellTopology.hpp"
49 
50 namespace panzer {
51 
53 {
54  const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
55  return ct.getSubcellCount(dim);
56 }
57 
58 const std::vector<int> & Intrepid2FieldPattern::getSubcellIndices(int dim,int cellIndex) const
59 {
60  const std::vector<std::vector<std::vector<int> > > & ordData = intrepidBasis_->getDofOrdinalData();
61 
62  // enough dimemnsions from intrepid?
63  if((std::size_t) dim<ordData.size()) {
64  // enough cells from intrepid?
65  if((std::size_t) cellIndex<ordData[dim].size()) {
66  if(ordData[dim][cellIndex].size()!=1)
67  return ordData[dim][cellIndex];
68  else if(ordData[dim][cellIndex][0]<0)
69  return empty_;
70  else
71  return ordData[dim][cellIndex];
72  }
73  }
74 
75  // not enough information!
76  return empty_;
77 }
78 
79 void Intrepid2FieldPattern::getSubcellClosureIndices(int dim,int cellIndex,std::vector<int> & indices) const
80 {
81  // recursive base case
82  if(dim==0) {
83  indices = getSubcellIndices(dim,cellIndex);
84  return;
85  }
86 
87  indices.clear(); // wipe out previous values...we are using insert!
88 
89  // use topology to build closure of sub cell
90  const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
91  std::set<std::pair<unsigned,unsigned> > closure;
92  Intrepid2FieldPattern::buildSubcellClosure(ct,dim,cellIndex,closure);
93 
94  // grab basis indices on the closure of the sub cell
95  std::set<std::pair<unsigned,unsigned> >::const_iterator itr;
96  for(itr=closure.begin();itr!=closure.end();++itr) {
97  // grab indices for this sub cell
98  const std::vector<int> & subcellIndices = getSubcellIndices(itr->first,itr->second);
99 
100  // concatenate the two vectors
101  indices.insert(indices.end(),subcellIndices.begin(),subcellIndices.end());
102  }
103 }
104 
106 {
107  const shards::CellTopology ct = intrepidBasis_->getBaseCellTopology();
108  return ct.getDimension();
109 }
110 
111 shards::CellTopology Intrepid2FieldPattern::getCellTopology() const
112 {
113  return intrepidBasis_->getBaseCellTopology();
114 }
115 
116 void Intrepid2FieldPattern::getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
117  std::vector<unsigned> & nodes)
118 {
119  if(dim==0) {
120  nodes.push_back(subCell);
121  return;
122  }
123 
124  // get all nodes on requested sub cell
125  unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
126  for(unsigned node=0;node<subCellNodeCount;++node)
127  nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));
128 
129  // sort them so they are ordered correctly for "includes" call
130  std::sort(nodes.begin(),nodes.end());
131 }
132 
133 void Intrepid2FieldPattern::findContainedSubcells(const shards::CellTopology & cellTopo,unsigned dim,
134  const std::vector<unsigned> & nodes,
135  std::set<std::pair<unsigned,unsigned> > & subCells)
136 {
137 
138  unsigned subCellCount = cellTopo.getSubcellCount(dim);
139  for(unsigned subCellOrd=0;subCellOrd<subCellCount;++subCellOrd) {
140  // get all nodes in sub cell
141  std::vector<unsigned> subCellNodes;
142  getSubcellNodes(cellTopo,dim,subCellOrd,subCellNodes);
143 
144  // if subCellNodes \subset nodes => add (dim,subCellOrd) to subCells
145  bool isSubset = std::includes( nodes.begin(), nodes.end(),
146  subCellNodes.begin(), subCellNodes.end());
147  if(isSubset)
148  subCells.insert(std::make_pair(dim,subCellOrd));
149 
150  }
151 
152  // stop recursion base case
153  if(dim==0) return;
154 
155  // find subcells in next sub dimension
156  findContainedSubcells(cellTopo,dim-1,nodes,subCells);
157 }
158 
159 void Intrepid2FieldPattern::buildSubcellClosure(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
160  std::set<std::pair<unsigned,unsigned> > & closure)
161 {
162  switch(dim) {
163  case 0:
164  closure.insert(std::make_pair(0,subCell));
165  break;
166  case 1:
167  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,0)));
168  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,1)));
169  closure.insert(std::make_pair(1,subCell));
170  break;
171  case 2:
172  {
173  unsigned cnt = (shards::CellTopology(cellTopo.getCellTopologyData(dim,subCell))).getSubcellCount(dim-1);
174  for(unsigned i=0;i<cnt;i++) {
175  int edge = mapCellFaceEdge(cellTopo.getCellTopologyData(),subCell,i);
176  buildSubcellClosure(cellTopo,dim-1,edge,closure);
177  }
178  closure.insert(std::make_pair(2,subCell));
179  }
180  break;
181  default:
182  // beyond a two dimension surface this thing crashes!
183  TEUCHOS_ASSERT(false);
184  };
185 }
186 
188 {
189  typedef Intrepid2::DofCoordsInterface<Kokkos::DynRankView<double,PHX::Device> > CoordsInterface;
190 
191  using Teuchos::RCP;
192  using Teuchos::rcp_dynamic_cast;
193 
194  // cast basis object to DofCoordsInterface: throw on failure
195  RCP<CoordsInterface> coordsInterface
196  = rcp_dynamic_cast<CoordsInterface>(intrepidBasis_,false);
197 
198  // if cast succeeds then coordinates are supported! Otherwise they are not!
199  return coordsInterface!=Teuchos::null;
200 }
201 
207 void Intrepid2FieldPattern::getInterpolatoryCoordinates(Kokkos::DynRankView<double,PHX::Device> & coords) const
208 {
209  typedef Intrepid2::DofCoordsInterface<Kokkos::DynRankView<double,PHX::Device> > CoordsInterface;
210 
211  using Teuchos::RCP;
212  using Teuchos::rcp_dynamic_cast;
213 
214  bool throwOnFail = true;
215 
216  // cast basis object to DofCoordsInterface: throw on failure
217  RCP<CoordsInterface> coordsInterface
218  = rcp_dynamic_cast<CoordsInterface>(intrepidBasis_,throwOnFail);
219 
220  // resize coordinates
221  coords = Kokkos::DynRankView<double,PHX::Device>("coords",intrepidBasis_->getCardinality(),getDimension());
222  coordsInterface->getDofCoords(coords);
223 }
224 
230 void Intrepid2FieldPattern::getInterpolatoryCoordinates(const Kokkos::DynRankView<double,PHX::Device> & cellVertices,
231  Kokkos::DynRankView<double,PHX::Device> & coords) const
232 {
233  TEUCHOS_ASSERT(cellVertices.rank()==3);
234 
235  int numCells = cellVertices.dimension(0);
236 
237  // grab the local coordinates
238  Kokkos::DynRankView<double,PHX::Device> localCoords;
239  getInterpolatoryCoordinates(localCoords);
240 
241  // resize the coordinates field container
242  coords = Kokkos::DynRankView<double,PHX::Device>("coords",numCells,localCoords.dimension(0),getDimension());
243 
244  if(numCells>0) {
245  // map to phsyical coordinates
246  Intrepid2::CellTools<double> cellTools;
247  cellTools.mapToPhysicalFrame(coords,localCoords,cellVertices,intrepidBasis_->getBaseCellTopology());
248  }
249 }
250 
251 }
252 #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)