Intrepid2
Intrepid2_OrientationToolsDefMatrixData.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
43 
48 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MATRIX_DATA_HPP__
49 #define __INTREPID2_ORIENTATIONTOOLS_DEF_MATRIX_DATA_HPP__
50 
51 // disable clang warnings
52 #if defined (__clang__) && !defined (__INTEL_COMPILER)
53 #pragma clang system_header
54 #endif
55 
56 namespace Intrepid2 {
57 
58  template<typename DT>
59  template<typename BasisHostType>
61  OrientationTools<DT>::createCoeffMatrixInternal(const BasisHostType* basis) {
62  const std::string name(basis->getName());
63  CoeffMatrixDataViewType matData;
64 
65  //
66  // High order HGRAD Elements
67  //
68  const auto cellTopo = basis->getBaseCellTopology();
69  const ordinal_type numEdges = cellTopo.getEdgeCount();
70  const ordinal_type numFaces = cellTopo.getFaceCount();
71  ordinal_type matDim = 0, matDim1 = 0, matDim2 = 0, numOrts = 0, numSubCells;
72  for(ordinal_type i=0; i<numEdges; ++i) {
73  matDim1 = std::max(matDim1, basis->getDofCount(1,i));
74  numOrts = std::max(numOrts,2);
75  }
76  for(ordinal_type i=0; i<numFaces; ++i) {
77  matDim2 = std::max(matDim2, basis->getDofCount(2,i));
78  numOrts = std::max(numOrts,2*ordinal_type(cellTopo.getSideCount(2,i)));
79  }
80  matDim = std::max(matDim1,matDim2);
81  numSubCells = (matDim1>0)*numEdges + (matDim2>0)*numFaces;
82 
83 
84  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::"+name,
85  numSubCells,
86  numOrts,
87  matDim,
88  matDim);
89 
90  if(basis->getFunctionSpace() == FUNCTION_SPACE_HGRAD) {
91  init_HGRAD(matData, basis);
92  } else if (basis->getFunctionSpace() == FUNCTION_SPACE_HCURL) {
93  init_HCURL(matData, basis);
94  } else if (basis->getFunctionSpace() == FUNCTION_SPACE_HDIV) {
95  init_HDIV(matData, basis);
96  }
97  return matData;
98  }
99 
100  //
101  // HGRAD elements
102  //
103  template<typename DT>
104  template<typename BasisHostType>
105  void
108  BasisHostType const *cellBasis) {
109 
110  const auto cellTopo = cellBasis->getBaseCellTopology();
111  const ordinal_type numEdges = cellTopo.getEdgeCount();
112  const ordinal_type numFaces = cellTopo.getFaceCount();
113 
114  {
115  const ordinal_type numOrt = 2;
116  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
117  if(cellBasis->getDofCount(1, edgeId) < 2) continue;
118  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
119  auto mat = Kokkos::subview(matData,
120  edgeId, edgeOrt,
121  Kokkos::ALL(), Kokkos::ALL());
123  (mat,
127  *cellBasis->getSubCellRefBasis(1,edgeId), *cellBasis,
128  edgeId, edgeOrt);
129  }
130  }
131  }
132  {
133 
134  for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
135  // this works for triangles (numOrt=6) and quadratures (numOrt=8)
136  const ordinal_type numOrt = 2*cellTopo.getSideCount(2,faceId);
137  if(cellBasis->getDofCount(2, faceId) < 1) continue;
138  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
139  auto mat = Kokkos::subview(matData,
140  numEdges+faceId, faceOrt,
141  Kokkos::ALL(), Kokkos::ALL());
146  (mat,
147  *cellBasis->getSubCellRefBasis(2,faceId), *cellBasis,
148  faceId, faceOrt);
149  }
150  }
151  }
152  }
153 
154  //
155  // HCURL elements
156  //
157  template<typename DT>
158  template<typename BasisHostType>
159  void
162  BasisHostType const *cellBasis) {
163  const auto cellTopo = cellBasis->getBaseCellTopology();
164  const ordinal_type numEdges = cellTopo.getEdgeCount();
165  const ordinal_type numFaces = cellTopo.getFaceCount();
166 
167  {
168  const ordinal_type numOrt = 2;
169  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
170  if(cellBasis->getDofCount(1, edgeId) < 1) continue;
171  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
172  auto mat = Kokkos::subview(matData,
173  edgeId, edgeOrt,
174  Kokkos::ALL(), Kokkos::ALL());
176  *cellBasis->getSubCellRefBasis(1,edgeId), *cellBasis,
177  edgeId, edgeOrt);
178  }
179  }
180  }
181  {
182  for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
183  // this works for triangles (numOrt=6) and quadratures (numOrt=8)
184  const ordinal_type numOrt = 2*cellTopo.getSideCount(2,faceId);
185  if(cellBasis->getDofCount(2, faceId) < 1) continue;
186  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
187  auto mat = Kokkos::subview(matData,
188  numEdges+faceId, faceOrt,
189  Kokkos::ALL(), Kokkos::ALL());
191  (mat,
192  *cellBasis->getSubCellRefBasis(2,faceId), *cellBasis,
193  faceId, faceOrt);
194  }
195  }
196  }
197  }
198 
199  //
200  // HDIV elements
201  //
202  template<typename DT>
203  template<typename BasisHostType>
204  void
207  BasisHostType const *cellBasis) {
208  const auto cellTopo = cellBasis->getBaseCellTopology();
209  const ordinal_type numSides = cellTopo.getSideCount();
210  const ordinal_type sideDim = cellTopo.getDimension()-1;
211 
212  {
213  for (ordinal_type sideId=0;sideId<numSides;++sideId) {
214  if(cellBasis->getDofCount(sideDim, sideId) < 1) continue;
215  const ordinal_type numOrt = (sideDim == 1) ? 2 : 2*cellTopo.getSideCount(sideDim,sideId);
216  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
217  auto mat = Kokkos::subview(matData,
218  sideId, faceOrt,
219  Kokkos::ALL(), Kokkos::ALL());
221  *cellBasis->getSubCellRefBasis(sideDim,sideId), *cellBasis,
222  sideId, faceOrt);
223  }
224  }
225  }
226  }
227 
228  template<typename DT>
229  template<typename BasisType>
231  OrientationTools<DT>::createCoeffMatrix(const BasisType* basis) {
232 #ifdef HAVE_INTREPID2_DEBUG
233  INTREPID2_TEST_FOR_EXCEPTION( !basis->requireOrientation(), std::invalid_argument,
234  ">>> ERROR (OrientationTools::createCoeffMatrix): basis does not require orientations." );
235 #endif
236  Kokkos::push_finalize_hook( [=] {
237  ortCoeffData=std::map<std::pair<std::string,ordinal_type>, typename OrientationTools<DT>::CoeffMatrixDataViewType>();
238  });
239 
240  const std::pair<std::string,ordinal_type> key(basis->getName(), basis->getDegree());
241  const auto found = ortCoeffData.find(key);
242 
243  CoeffMatrixDataViewType matData;
244  if (found == ortCoeffData.end()) {
245  {
246  auto basis_host = basis->getHostBasis();
247  matData = createCoeffMatrixInternal(basis_host.getRawPtr());
248  ortCoeffData.insert(std::make_pair(key, matData));
249  }
250  } else {
251  matData = found->second;
252  }
253 
254  return matData;
255  }
256 
257  template<typename DT>
259  ortCoeffData.clear();
260  }
261 }
262 
263 #endif
Kokkos::View< double ****, DeviceType > CoeffMatrixDataViewType
subcell ordinal, orientation, matrix m x n
static void getCoeffMatrix_HDIV(OutputViewType &output, const subcellBasisHostType &subcellBasis, const cellBasisHostType &cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute orientation matrix for HDIV basis for a given subcell and its reference basis.
static void getCoeffMatrix_HCURL(OutputViewType &output, const subcellBasisHostType &subcellBasis, const cellBasisHostType &cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute orientation matrix for HCURL basis for a given subcell and its reference basis.
static void clearCoeffMatrix()
Clear coefficient matrix.
static void init_HCURL(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HCURL basis.
static void getCoeffMatrix_HGRAD(OutputViewType &output, const subcellBasisHostType &subcellBasis, const cellBasisHostType &cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute orientation matrix for HGRAD basis for a given subcell and its reference basis.
static CoeffMatrixDataViewType createCoeffMatrix(const BasisType *basis)
Create coefficient matrix.
static void init_HGRAD(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HGRAD basis.
static void init_HDIV(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HDIV basis.