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
56namespace 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.getSubcellCount(1);
70 const ordinal_type numFaces = cellTopo.getSubcellCount(2);
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 } else if (basis->getFunctionSpace() == FUNCTION_SPACE_HVOL) {
97 init_HVOL(matData, basis);
98 }
99 return matData;
100 }
101
102 //
103 // HGRAD elements
104 //
105 template<typename DT>
106 template<typename BasisHostType>
107 void
110 BasisHostType const *cellBasis) {
111
112 const auto cellTopo = cellBasis->getBaseCellTopology();
113 const ordinal_type numEdges = cellTopo.getSubcellCount(1);
114 const ordinal_type numFaces = cellTopo.getSubcellCount(2);
115 Intrepid2::BasisPtr<typename BasisHostType::DeviceType,
116 typename BasisHostType::OutputValueType,
117 typename BasisHostType::PointValueType>
118 basisPtr;
119 BasisHostType const *subcellBasis;
120
121 { //edges
122 subcellBasis = cellBasis; // if (dim==1)
123 const ordinal_type numOrt = 2;
124 for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
125 if(cellBasis->getDofCount(1, edgeId) < 2) continue;
126 if(cellTopo.getDimension()!=1) {
127 basisPtr = cellBasis->getSubCellRefBasis(1,edgeId);
128 subcellBasis = basisPtr.get();
129 }
130
131 for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
132 auto mat = Kokkos::subview(matData,
133 edgeId, edgeOrt,
134 Kokkos::ALL(), Kokkos::ALL());
136 (mat,
137 *subcellBasis, *cellBasis,
138 edgeId, edgeOrt);
139 }
140 }
141 }
142 { //faces
143 subcellBasis = cellBasis; // if(dim==2)
144 for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
145 // this works for triangles (numOrt=6) and quadratures (numOrt=8)
146 const ordinal_type numOrt = 2*cellTopo.getSideCount(2,faceId);
147 if(cellBasis->getDofCount(2, faceId) < 1) continue;
148 if(cellTopo.getDimension()!=2) {
149 basisPtr = cellBasis->getSubCellRefBasis(2,faceId);
150 subcellBasis = basisPtr.get();
151 }
152 for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
153 auto mat = Kokkos::subview(matData,
154 numEdges+faceId, faceOrt,
155 Kokkos::ALL(), Kokkos::ALL());
157 (mat,
158 *subcellBasis, *cellBasis,
159 faceId, faceOrt);
160 }
161 }
162 }
163 }
164
165 //
166 // HCURL elements
167 //
168 template<typename DT>
169 template<typename BasisHostType>
170 void
173 BasisHostType const *cellBasis) {
174 const auto cellTopo = cellBasis->getBaseCellTopology();
175 const ordinal_type numEdges = cellTopo.getSubcellCount(1);
176 const ordinal_type numFaces = cellTopo.getSubcellCount(2);
177 Intrepid2::BasisPtr<typename BasisHostType::DeviceType,
178 typename BasisHostType::OutputValueType,
179 typename BasisHostType::PointValueType>
180 basisPtr;
181 BasisHostType const* subcellBasis;
182
183 { // edges
184 subcellBasis = cellBasis; // if (dim==1)
185 const ordinal_type numOrt = 2;
186 for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
187 if(cellBasis->getDofCount(1, edgeId) < 1) continue;
188 if(cellTopo.getDimension()!=1) {
189 basisPtr = cellBasis->getSubCellRefBasis(1,edgeId);
190 subcellBasis = basisPtr.get();
191 }
192 for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
193 auto mat = Kokkos::subview(matData,
194 edgeId, edgeOrt,
195 Kokkos::ALL(), Kokkos::ALL());
197 *subcellBasis, *cellBasis,
198 edgeId, edgeOrt);
199 }
200 }
201 }
202 { //faces
203 subcellBasis = cellBasis; // if (dim==2)
204 for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
205 // this works for triangles (numOrt=6) and quadratures (numOrt=8)
206 const ordinal_type numOrt = 2*cellTopo.getSideCount(2,faceId);
207 if(cellBasis->getDofCount(2, faceId) < 1) continue;
208 if(cellTopo.getDimension()!=2) {
209 basisPtr = cellBasis->getSubCellRefBasis(2,faceId);
210 subcellBasis = basisPtr.get();
211 }
212 for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
213 auto mat = Kokkos::subview(matData,
214 numEdges+faceId, faceOrt,
215 Kokkos::ALL(), Kokkos::ALL());
217 (mat,
218 *subcellBasis, *cellBasis,
219 faceId, faceOrt);
220 }
221 }
222 }
223 }
224
225 //
226 // HDIV elements
227 //
228 template<typename DT>
229 template<typename BasisHostType>
230 void
233 BasisHostType const *cellBasis) {
234 const auto cellTopo = cellBasis->getBaseCellTopology();
235 const ordinal_type numSides = cellTopo.getSideCount();
236 const ordinal_type sideDim = cellTopo.getDimension()-1;
237 Intrepid2::BasisPtr<typename BasisHostType::DeviceType,
238 typename BasisHostType::OutputValueType,
239 typename BasisHostType::PointValueType>
240 subcellBasisPtr;
241
242 {
243 for (ordinal_type sideId=0;sideId<numSides;++sideId) {
244 if(cellBasis->getDofCount(sideDim, sideId) < 1) continue;
245 const ordinal_type numOrt = (sideDim == 1) ? 2 : 2*cellTopo.getSideCount(sideDim,sideId);
246 subcellBasisPtr = cellBasis->getSubCellRefBasis(sideDim,sideId);
247 for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
248 auto mat = Kokkos::subview(matData,
249 sideId, faceOrt,
250 Kokkos::ALL(), Kokkos::ALL());
252 *subcellBasisPtr, *cellBasis,
253 sideId, faceOrt);
254 }
255 }
256 }
257 }
258
259 //
260 // HVOL elements (used for 2D and 1D side cells)
261 //
262 template<typename DT>
263 template<typename BasisHostType>
264 void
267 BasisHostType const *cellBasis) {
268
269 const auto cellTopo = cellBasis->getBaseCellTopology();
270 const ordinal_type numEdges = (cellTopo.getDimension()==1);
271 const ordinal_type numFaces = (cellTopo.getDimension()==2);
272
273 {
274 const ordinal_type numOrt = 2;
275 for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
276 if(cellBasis->getDofCount(1, edgeId) < 1) continue;
277 for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
278 auto mat = Kokkos::subview(matData,
279 edgeId, edgeOrt,
280 Kokkos::ALL(), Kokkos::ALL());
282 (mat, *cellBasis, edgeOrt);
283 }
284 }
285 }
286 {
287 for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
288 // this works for triangles (numOrt=6) and quadratures (numOrt=8)
289 const ordinal_type numOrt = 2*cellTopo.getSideCount(2,faceId);
290 if(cellBasis->getDofCount(2, faceId) < 1) continue;
291 for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
292 auto mat = Kokkos::subview(matData,
293 numEdges+faceId, faceOrt,
294 Kokkos::ALL(), Kokkos::ALL());
296 (mat, *cellBasis, faceOrt);
297 }
298 }
299 }
300 }
301
302 template<typename DT>
303 template<typename BasisType>
306 Kokkos::push_finalize_hook( [=] {
307 ortCoeffData=std::map<std::pair<std::string,ordinal_type>, typename OrientationTools<DT>::CoeffMatrixDataViewType>();
308 });
309
310 const std::pair<std::string,ordinal_type> key(basis->getName(), basis->getDegree());
311 const auto found = ortCoeffData.find(key);
312
314 if (found == ortCoeffData.end()) {
315 {
316 auto basis_host = basis->getHostBasis();
317 matData = createCoeffMatrixInternal(basis_host.getRawPtr());
318 ortCoeffData.insert(std::make_pair(key, matData));
319 }
320 } else {
321 matData = found->second;
322 }
323
324 return matData;
325 }
326
327 template<typename DT>
329 ortCoeffData.clear();
330 }
331}
332
333#endif
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
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 void getCoeffMatrix_HVOL(OutputViewType &output, const cellBasisHostType &cellBasis, const ordinal_type cellOrt)
Compute orientation matrix for HVOL basis for a given (2D or 1D) cell 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 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 CoeffMatrixDataViewType createCoeffMatrix(const BasisType *basis)
Create coefficient matrix.
static void init_HDIV(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HDIV basis.
Kokkos::View< double ****, DeviceType > CoeffMatrixDataViewType
subcell ordinal, orientation, matrix m x n
static void init_HVOL(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HVOL basis.
static void clearCoeffMatrix()
Clear coefficient matrix.
static void init_HCURL(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HCURL basis.
static void init_HGRAD(CoeffMatrixDataViewType matData, BasisHostType const *cellBasis)
Compute orientation matrix for HGRAD basis.