Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
49 namespace 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
128  {
129  return intrepidBasis_->getBaseCellTopology().getDimension();
130  }
131 
132  shards::CellTopology
135  {
136  // TODO BWR Probably should change the name of this call...
137  // TODO BWR However this is a virtual function
138  return intrepidBasis_->getBaseCellTopology();
139  }
140 
141  void
143  getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
144  std::vector<unsigned> & nodes)
145  {
146  if(dim==0) {
147  nodes.push_back(subCell);
148  return;
149  }
150 
151  // get all nodes on requested sub cell
152  unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
153  for(unsigned node=0;node<subCellNodeCount;++node)
154  nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));
155 
156  // sort them so they are ordered correctly for "includes" call
157  std::sort(nodes.begin(),nodes.end());
158  }
159 
160  void
162  findContainedSubcells(const shards::CellTopology & cellTopo,unsigned dim,
163  const std::vector<unsigned> & nodes,
164  std::set<std::pair<unsigned,unsigned> > & subCells)
165  {
166  unsigned subCellCount = cellTopo.getSubcellCount(dim);
167  for(unsigned subCellOrd=0;subCellOrd<subCellCount;++subCellOrd) {
168  // get all nodes in sub cell
169  std::vector<unsigned> subCellNodes;
170  getSubcellNodes(cellTopo,dim,subCellOrd,subCellNodes);
171 
172  // if subCellNodes \subset nodes => add (dim,subCellOrd) to subCells
173  bool isSubset = std::includes( nodes.begin(), nodes.end(),
174  subCellNodes.begin(), subCellNodes.end());
175  if(isSubset)
176  subCells.insert(std::make_pair(dim,subCellOrd));
177 
178  }
179 
180  // stop recursion base case
181  if (dim==0) return;
182 
183  // find subcells in next sub dimension
184  findContainedSubcells(cellTopo,dim-1,nodes,subCells);
185  }
186 
187  void
189  buildSubcellClosure(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
190  std::set<std::pair<unsigned,unsigned> > & closure)
191  {
192  switch(dim) {
193  case 0:
194  closure.insert(std::make_pair(0,subCell));
195  break;
196  case 1:
197  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,0)));
198  closure.insert(std::make_pair(0,cellTopo.getNodeMap(dim,subCell,1)));
199  closure.insert(std::make_pair(1,subCell));
200  break;
201  case 2:
202  {
203  unsigned cnt = (shards::CellTopology(cellTopo.getCellTopologyData(dim,subCell))).getSubcellCount(dim-1);
204  for(unsigned i=0;i<cnt;i++) {
205  int edge = mapCellFaceEdge(cellTopo.getCellTopologyData(),subCell,i);
206  buildSubcellClosure(cellTopo,dim-1,edge,closure);
207  }
208  closure.insert(std::make_pair(2,subCell));
209  }
210  break;
211  default:
212  // beyond a two dimension surface this thing crashes!
213  TEUCHOS_ASSERT(false);
214  };
215  }
216 
217  bool
220  {
221  // we no longer use CoordsInterface
222  return true;
223  }
224 
230  void
232  getInterpolatoryCoordinates(Kokkos::DynRankView<double,PHX::Device> & coords) const
233  {
234  // this may not be efficient if coords is allocated every time this function is called
235  coords = Kokkos::DynRankView<double,PHX::Device>("coords",intrepidBasis_->getCardinality(),getDimension());
236  intrepidBasis_->getDofCoords(coords);
237  }
238 
244  void
246  getInterpolatoryCoordinates(const Kokkos::DynRankView<double,PHX::Device> & cellNodes,
247  Kokkos::DynRankView<double,PHX::Device> & coords,
248  Teuchos::RCP<const shards::CellTopology> meshCellTopology) const
249  {
250  TEUCHOS_ASSERT(cellNodes.rank()==3);
251 
252  int numCells = cellNodes.extent(0);
253 
254  // grab the local coordinates
255  Kokkos::DynRankView<double,PHX::Device> localCoords;
256  getInterpolatoryCoordinates(localCoords);
257 
258  // resize the coordinates field container
259  coords = Kokkos::DynRankView<double,PHX::Device>("coords",numCells,localCoords.extent(0),getDimension());
260 
261  if(numCells>0) {
262  Intrepid2::CellTools<PHX::Device> cellTools;
263  // For backwards compatability, allow the FE basis to supply the mesh cell topology (via the FEM base cell topo)
264  // This will occur if no meshCellTopology is provided (defaults to Teuchos::null)
265  // If provided, use the mesh topology directly
266  if (meshCellTopology==Teuchos::null) {
267  cellTools.mapToPhysicalFrame(coords,localCoords,cellNodes,intrepidBasis_->getBaseCellTopology());
268  } else {
269  cellTools.mapToPhysicalFrame(coords,localCoords,cellNodes,*meshCellTopology);
270  }
271  }
272  }
273 
276  { return intrepidBasis_; }
277 
278 }
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)
virtual shards::CellTopology getCellTopology() const
bool supportsInterpolatoryCoordinates() const
Does this field pattern support interpolatory coordinates?
Intrepid2FieldPattern(const Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > &intrepidBasis)
static void getSubcellNodes(const shards::CellTopology &cellTopo, unsigned dim, unsigned subCell, std::vector< unsigned > &nodes)
virtual void getSubcellClosureIndices(int dim, int cellIndex, std::vector< int > &indices) const
Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > getIntrepidBasis() const
Returns the underlying Intrepid2::Basis object.
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)
void getInterpolatoryCoordinates(Kokkos::DynRankView< double, PHX::Device > &coords) const
#define TEUCHOS_ASSERT(assertion_test)
virtual int getSubcellCount(int dim) const
Teuchos::RCP< Intrepid2::Basis< PHX::Device, double, double > > intrepidBasis_