Intrepid2
Intrepid2_ProjectionStructDef.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
47#ifndef __INTREPID2_PROJECTIONSTRUCTDEF_HPP__
48#define __INTREPID2_PROJECTIONSTRUCTDEF_HPP__
49
50#include "Intrepid2_Utils.hpp"
51#include "Shards_CellTopology.hpp"
52#include "Shards_BasicTopologies.hpp"
53
58
60
61
62namespace Intrepid2 {
63
64namespace Experimental {
65template<typename DeviceType, typename ValueType>
66template<typename BasisPtrType>
68 const ordinal_type targetCubDegree) {
69 const auto cellTopo = cellBasis->getBaseCellTopology();
70 std::string name = cellBasis->getName();
71 ordinal_type dim = cellTopo.getDimension();
72 numBasisEvalPoints = 0;
73 numBasisDerivEvalPoints = 0;
74 numTargetEvalPoints = 0;
75 numTargetDerivEvalPoints = 0;
76 const ordinal_type edgeDim = 1;
77 const ordinal_type faceDim = 2;
78
79 ordinal_type basisCubDegree = cellBasis->getDegree();
80 ordinal_type edgeBasisCubDegree, faceBasisCubDegree;
81
82 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
83 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
84 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
85
86 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
87
88 INTREPID2_TEST_FOR_ABORT( (numVertices > maxSubCellsCount) || (numFaces > maxSubCellsCount) || (numEdges > maxSubCellsCount),
89 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with this cub cells count");
90
91
92 numBasisEvalPoints += numVertices;
93 numTargetEvalPoints += numVertices;
94 view_type coord("vertex_coord", dim);
95
96 basisPointsRange = range_tag("basisPointsRange", 4,maxSubCellsCount);
97 targetPointsRange = range_tag("targetPointsRange", 4,maxSubCellsCount);
98 basisDerivPointsRange = range_tag("basisDerivPointsRange", 4,maxSubCellsCount);
99 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
100 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
101
102 maxNumBasisEvalPoints = numVertices; maxNumTargetEvalPoints = numVertices;
103 for(ordinal_type iv=0; iv<numVertices; ++iv) {
104 basisPointsRange(0,iv) = range_type(iv, iv+1);
105 basisCubPoints[0][iv] = view_type("basisCubPoints",1,dim);
106 targetPointsRange(0,iv) = range_type(iv, iv+1);
107 targetCubPoints[0][iv] = view_type("targetCubPoints",1,dim);
109 for(ordinal_type d=0; d<dim; d++) {
110 basisCubPoints[0][iv](0,d) = coord(d);
111 targetCubPoints[0][iv](0,d) = coord(d);
112 }
113 }
114
115 if (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HCURL) {
116 edgeBasisCubDegree = basisCubDegree - 1;
117 faceBasisCubDegree = basisCubDegree;
118 }
119 else if (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HDIV) {
120 edgeBasisCubDegree = basisCubDegree - 1;
121 faceBasisCubDegree = basisCubDegree -1;
122 }
123 else {
124 edgeBasisCubDegree = basisCubDegree;
125 faceBasisCubDegree = basisCubDegree;
126 }
127
128 DefaultCubatureFactory cub_factory;
129 for(ordinal_type ie=0; ie<numEdges; ++ie) {
130 ordinal_type cub_degree = 2*edgeBasisCubDegree;
131 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
132 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
133 basisPointsRange(edgeDim,ie) = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
134 numBasisEvalPoints += edgeBasisCub->getNumPoints();
135 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, edgeBasisCub->getNumPoints());
136 basisCubPoints[edgeDim][ie] = view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
137 basisCubWeights[edgeDim][ie] = view_type("basisCubWeights",edgeBasisCub->getNumPoints());
138 edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
139
140 cub_degree = edgeBasisCubDegree + targetCubDegree;
141 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
142 targetPointsRange(edgeDim,ie) = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
143 numTargetEvalPoints += edgeTargetCub->getNumPoints();
144 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, edgeTargetCub->getNumPoints());
145 targetCubPoints[edgeDim][ie] = view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
146 targetCubWeights[edgeDim][ie] = view_type("targetCubWeights",edgeTargetCub->getNumPoints());
147 edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
148 }
149
150 for(ordinal_type iface=0; iface<numFaces; ++iface) {
151 ordinal_type cub_degree = 2*faceBasisCubDegree;
152 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
153 auto faceBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
154 basisPointsRange(faceDim,iface) = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
155 numBasisEvalPoints += faceBasisCub->getNumPoints();
156 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, faceBasisCub->getNumPoints());
157 basisCubPoints[faceDim][iface] = view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
158 basisCubWeights[faceDim][iface] = view_type("basisCubWeights",faceBasisCub->getNumPoints());
159 faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
160
161 cub_degree = faceBasisCubDegree + targetCubDegree;
162 auto faceTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
163 targetPointsRange(faceDim,iface) = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
164 numTargetEvalPoints += faceTargetCub->getNumPoints();
165 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, faceTargetCub->getNumPoints());
166 targetCubPoints[faceDim][iface] = view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
167 targetCubWeights[faceDim][iface] = view_type("targetCubWeights",faceTargetCub->getNumPoints());
168 faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
169 }
170 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
171 if(hasCellDofs) {
172 ordinal_type cub_degree = 2*basisCubDegree;
173 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
174 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
175 numBasisEvalPoints += elemBasisCub->getNumPoints();
176 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
177 basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
178 basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
179 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
180
181 cub_degree = basisCubDegree + targetCubDegree;
182 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
183 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
184 numTargetEvalPoints += elemTargetCub->getNumPoints();
185 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
186 targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
187 targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
188 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
189 }
190}
191
192template<typename DeviceType, typename ValueType>
193template<typename BasisPtrType>
195 const ordinal_type targetCubDegree,
196 const ordinal_type targetGradCubDegree) {
197 const auto cellTopo = cellBasis->getBaseCellTopology();
198 std::string name = cellBasis->getName();
199 ordinal_type dim = cellTopo.getDimension();
200 numBasisEvalPoints = 0;
201 numBasisDerivEvalPoints = 0;
202 numTargetEvalPoints = 0;
203 numTargetDerivEvalPoints = 0;
204 const ordinal_type edgeDim = 1;
205 const ordinal_type faceDim = 2;
206
207 ordinal_type basisCubDegree = cellBasis->getDegree();
208 ordinal_type edgeBasisCubDegree = basisCubDegree;
209 ordinal_type faceBasisCubDegree = basisCubDegree;
210
211 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
212 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount(): 0;
213 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
214
215 INTREPID2_TEST_FOR_ABORT( (numFaces > maxSubCellsCount) || (numEdges > maxSubCellsCount),
216 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with this cub cells count");
217
218
219 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
220
221 maxNumBasisEvalPoints = numVertices; maxNumTargetEvalPoints = numVertices;
222 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
223
224 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
225 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
226 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
227 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
228 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
229
230 numBasisEvalPoints += numVertices;
231 numTargetEvalPoints += numVertices;
232 view_type coord("vertex_coord", dim);
233 for(ordinal_type iv=0; iv<numVertices; ++iv) {
234 basisPointsRange(0,iv) = range_type(iv, iv+1);
235 basisCubPoints[0][iv] = view_type("basisCubPoints",1,dim);
236 targetPointsRange(0,iv) = range_type(iv, iv+1);
237 targetCubPoints[0][iv] = view_type("targetCubPoints",1,dim);
239 for(ordinal_type d=0; d<dim; d++) {
240 basisCubPoints[0][iv](0,d) = coord(d);
241 targetCubPoints[0][iv](0,d) = coord(d);
242 }
243 }
244
245 DefaultCubatureFactory cub_factory;
246 for(ordinal_type ie=0; ie<numEdges; ++ie) {
247 ordinal_type cub_degree = 2*edgeBasisCubDegree;
248 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
249 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
250 basisDerivPointsRange(edgeDim,ie) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+edgeBasisCub->getNumPoints());
251 numBasisDerivEvalPoints += edgeBasisCub->getNumPoints();
252 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, edgeBasisCub->getNumPoints());
253 basisDerivCubPoints[edgeDim][ie] = view_type("basisDerivCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
254 basisDerivCubWeights[edgeDim][ie] = view_type("basisDerivCubWeights",edgeBasisCub->getNumPoints());
255 edgeBasisCub->getCubature(basisDerivCubPoints[edgeDim][ie], basisDerivCubWeights[edgeDim][ie]);
256
257 cub_degree = edgeBasisCubDegree + targetGradCubDegree;
258 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
259 targetDerivPointsRange(edgeDim,ie) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+edgeTargetCub->getNumPoints());
260 numTargetDerivEvalPoints += edgeTargetCub->getNumPoints();
261 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, edgeTargetCub->getNumPoints());
262 targetDerivCubPoints[edgeDim][ie] = view_type("targetDerivCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
263 targetDerivCubWeights[edgeDim][ie] = view_type("targetDerivCubWeights",edgeTargetCub->getNumPoints());
264 edgeTargetCub->getCubature(targetDerivCubPoints[edgeDim][ie], targetDerivCubWeights[edgeDim][ie]);
265 }
266
267 for(ordinal_type iface=0; iface<numFaces; ++iface) {
268 ordinal_type cub_degree = 2*faceBasisCubDegree;
269 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
270 auto faceBasisGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
271 basisDerivPointsRange(faceDim,iface) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisGradCub->getNumPoints());
272 numBasisDerivEvalPoints += faceBasisGradCub->getNumPoints();
273 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, faceBasisGradCub->getNumPoints());
274 basisDerivCubPoints[faceDim][iface] = view_type("basisDerivCubPoints",faceBasisGradCub->getNumPoints(),faceDim);
275 basisDerivCubWeights[faceDim][iface] = view_type("basisDerivCubWeights",faceBasisGradCub->getNumPoints());
276 faceBasisGradCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
277
278 cub_degree = faceBasisCubDegree + targetGradCubDegree;
279 auto faceTargetDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
280 targetDerivPointsRange(faceDim,iface) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
281 numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
282 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, faceTargetDerivCub->getNumPoints());
283 targetDerivCubPoints[faceDim][iface] = view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
284 targetDerivCubWeights[faceDim][iface] = view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
285 faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
286 }
287 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
288 if(hasCellDofs) {
289 ordinal_type cub_degree = 2*basisCubDegree;
290 auto elemBasisGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
291 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisGradCub->getNumPoints());
292 numBasisDerivEvalPoints += elemBasisGradCub->getNumPoints();
293 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisGradCub->getNumPoints());
294 basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisGradCub->getNumPoints(),dim);
295 basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisGradCub->getNumPoints());
296 elemBasisGradCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
297
298 cub_degree = basisCubDegree + targetGradCubDegree;
299 auto elemTargetGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
300 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetGradCub->getNumPoints());
301 numTargetDerivEvalPoints += elemTargetGradCub->getNumPoints();
302 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetGradCub->getNumPoints());
303 targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetGradCub->getNumPoints(),dim);
304 targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetGradCub->getNumPoints());
305 elemTargetGradCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
306 }
307}
308
309template<typename DeviceType, typename ValueType>
310template<typename BasisPtrType>
312 const ordinal_type targetCubDegree,
313 const ordinal_type targetCurlCubDegre) {
314 const auto cellTopo = cellBasis->getBaseCellTopology();
315 std::string name = cellBasis->getName();
316 ordinal_type dim = cellTopo.getDimension();
317 numBasisEvalPoints = 0;
318 numBasisDerivEvalPoints = 0;
319 numTargetEvalPoints = 0;
320 numTargetDerivEvalPoints = 0;
321 const ordinal_type edgeDim = 1;
322 const ordinal_type faceDim = 2;
323
324 ordinal_type basisCubDegree = cellBasis->getDegree();
325 ordinal_type edgeBasisCubDegree = basisCubDegree - 1;
326 ordinal_type faceBasisCubDegree = basisCubDegree;
327 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
328 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
329 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
330
331 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints = 0;
332 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
333
334 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
335 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
336 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
337 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
338 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
339
340 DefaultCubatureFactory cub_factory;
341 for(ordinal_type ie=0; ie<numEdges; ++ie) {
342 ordinal_type cub_degree = 2*edgeBasisCubDegree;
343 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
344 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(edgeDim,ie), cub_degree);
345 basisPointsRange(edgeDim,ie) = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
346 numBasisEvalPoints += edgeBasisCub->getNumPoints();
347 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, edgeBasisCub->getNumPoints());
348 basisCubPoints[edgeDim][ie] = view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
349 basisCubWeights[edgeDim][ie] = view_type("basisCubWeights",edgeBasisCub->getNumPoints());
350 edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
351
352 cub_degree = edgeBasisCubDegree + targetCubDegree;
353 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(edgeDim,ie), cub_degree);
354 targetPointsRange(edgeDim,ie) = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
355 numTargetEvalPoints += edgeTargetCub->getNumPoints();
356 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, edgeTargetCub->getNumPoints());
357 targetCubPoints[edgeDim][ie] = view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
358 targetCubWeights[edgeDim][ie] = view_type("targetCubWeights",edgeTargetCub->getNumPoints());
359 edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
360 }
361
362 for(ordinal_type iface=0; iface<numFaces; ++iface) {
363 ordinal_type cub_degree = 2*faceBasisCubDegree;
364 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
365 auto faceBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
366 basisPointsRange(faceDim,iface) = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
367 numBasisEvalPoints += faceBasisCub->getNumPoints();
368 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, faceBasisCub->getNumPoints());
369 basisCubPoints[faceDim][iface] = view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
370 basisCubWeights[faceDim][iface] = view_type("basisCubWeights",faceBasisCub->getNumPoints());
371 faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
372
373 auto faceBasisDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
374 basisDerivPointsRange(faceDim,iface) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisCub->getNumPoints());
375 numBasisDerivEvalPoints += faceBasisCub->getNumPoints();
376 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, faceBasisCub->getNumPoints());
377 basisDerivCubPoints[faceDim][iface] = view_type("basisDerivCubPoints",faceBasisCub->getNumPoints(),faceDim);
378 basisDerivCubWeights[faceDim][iface] = view_type("basisDerivCubWeights",faceBasisCub->getNumPoints());
379 faceBasisCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
380
381 cub_degree = faceBasisCubDegree + targetCubDegree;
382 auto faceTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
383 targetPointsRange(faceDim,iface) = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
384 numTargetEvalPoints += faceTargetCub->getNumPoints();
385 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, faceTargetCub->getNumPoints());
386 targetCubPoints[faceDim][iface] = view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
387 targetCubWeights[faceDim][iface] = view_type("targetCubWeights",faceTargetCub->getNumPoints());
388 faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
389
390 cub_degree = faceBasisCubDegree + targetCurlCubDegre;
391 auto faceTargetDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree);
392 targetDerivPointsRange(faceDim,iface) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
393 numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
394 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, faceTargetDerivCub->getNumPoints());
395 targetDerivCubPoints[faceDim][iface] = view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
396 targetDerivCubWeights[faceDim][iface] = view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
397 faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
398 }
399
400 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
401 if(hasCellDofs) {
402 ordinal_type cub_degree = 2*basisCubDegree;
403 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
404 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
405 numBasisEvalPoints += elemBasisCub->getNumPoints();
406 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
407 basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
408 basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
409 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
410
411 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisCub->getNumPoints());
412 numBasisDerivEvalPoints += elemBasisCub->getNumPoints();
413 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisCub->getNumPoints());
414 basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisCub->getNumPoints(),dim);
415 basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisCub->getNumPoints());
416 elemBasisCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
417
418 cub_degree = basisCubDegree + targetCubDegree;
419 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
420 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
421 numTargetEvalPoints += elemTargetCub->getNumPoints();
422 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
423 targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
424 targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
425 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
426
427 cub_degree = basisCubDegree + targetCurlCubDegre;
428 auto elemTargetCurlCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
429 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetCurlCub->getNumPoints());
430 numTargetDerivEvalPoints += elemTargetCurlCub->getNumPoints();
431 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetCurlCub->getNumPoints());
432 targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetCurlCub->getNumPoints(),dim);
433 targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetCurlCub->getNumPoints());
434 elemTargetCurlCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
435 }
436}
437
438template<typename DeviceType, typename ValueType>
439template<typename BasisPtrType>
441 const ordinal_type targetCubDegree,
442 const ordinal_type targetDivCubDegre) {
443 const auto cellTopo = cellBasis->getBaseCellTopology();
444 std::string name = cellBasis->getName();
445 ordinal_type dim = cellTopo.getDimension();
446 numBasisEvalPoints = 0;
447 numBasisDerivEvalPoints = 0;
448 numTargetEvalPoints = 0;
449 numTargetDerivEvalPoints = 0;
450 const ordinal_type sideDim = dim - 1;
451
452 ordinal_type basisCubDegree = cellBasis->getDegree();
453 ordinal_type sideBasisCubDegree = basisCubDegree - 1;
454 ordinal_type numSides = cellTopo.getSideCount()*ordinal_type(cellBasis->getDofCount(sideDim, 0) > 0);
455 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
456
457 INTREPID2_TEST_FOR_ABORT( numSides > maxSubCellsCount,
458 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with so many sides");
459
460 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints = 0;
461 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
462
463 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
464 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
465 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
466 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
467 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
468
470 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
471 hcurlBasis = new Basis_HCURL_HEX_In_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
472 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
473 hcurlBasis = new Basis_HCURL_TET_In_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
474 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
475 hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
476 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
477 hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
478 else {
479 std::stringstream ss;
480 ss << ">>> ERROR (Intrepid2::ProjectionTools::createHDivProjectionStruct): "
481 << "Method not implemented for basis " << name;
482 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
483 }
484
485 bool haveHCurlConstraint = (hcurlBasis != NULL) && (hcurlBasis->getDofCount(dim,0) > 0);
486 if(hcurlBasis != NULL) delete hcurlBasis;
487
488 DefaultCubatureFactory cub_factory;
489
490 for(ordinal_type is=0; is<numSides; ++is) {
491 ordinal_type cub_degree = 2*sideBasisCubDegree;
492 subCellTopologyKey(sideDim,is) = cellBasis->getBaseCellTopology().getKey(sideDim, is);
493 auto sideBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(sideDim,is), cub_degree);
494 basisPointsRange(sideDim,is) = range_type(numBasisEvalPoints, numBasisEvalPoints+sideBasisCub->getNumPoints());
495 numBasisEvalPoints += sideBasisCub->getNumPoints();
496 basisCubPoints[sideDim][is] = view_type("basisCubPoints",sideBasisCub->getNumPoints(),sideDim);
497 basisCubWeights[sideDim][is] = view_type("basisCubWeights",sideBasisCub->getNumPoints());
498 sideBasisCub->getCubature(basisCubPoints[sideDim][is], basisCubWeights[sideDim][is]);
499 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, sideBasisCub->getNumPoints());
500
501 cub_degree = sideBasisCubDegree + targetCubDegree;
502 auto sideTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(sideDim,is), cub_degree);
503 targetPointsRange(sideDim,is) = range_type(numTargetEvalPoints, numTargetEvalPoints+sideTargetCub->getNumPoints());
504 numTargetEvalPoints += sideTargetCub->getNumPoints();
505 targetCubPoints[sideDim][is] = view_type("targetCubPoints",sideTargetCub->getNumPoints(),sideDim);
506 targetCubWeights[sideDim][is] = view_type("targetCubWeights",sideTargetCub->getNumPoints());
507 sideTargetCub->getCubature(targetCubPoints[sideDim][is], targetCubWeights[sideDim][is]);
508 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, sideTargetCub->getNumPoints());
509 }
510
511 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
512 if(hasCellDofs) {
513 ordinal_type cub_degree = 2*basisCubDegree - 1;
514 auto elemBasisDivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
515 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisDivCub->getNumPoints());
516 numBasisDerivEvalPoints += elemBasisDivCub->getNumPoints();
517 basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisDivCub->getNumPoints(),dim);
518 basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisDivCub->getNumPoints());
519 elemBasisDivCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
520 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisDivCub->getNumPoints());
521
522 cub_degree = basisCubDegree - 1 + targetDivCubDegre;
523 auto elemTargetDivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
524 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetDivCub->getNumPoints());
525 numTargetDerivEvalPoints += elemTargetDivCub->getNumPoints();
526 targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetDivCub->getNumPoints(),dim);
527 targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetDivCub->getNumPoints());
528 elemTargetDivCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
529 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetDivCub->getNumPoints());
530
531 if(haveHCurlConstraint)
532 {
533 cub_degree = 2*basisCubDegree;
534 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
535 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints + elemBasisCub->getNumPoints());
536 numBasisEvalPoints += elemBasisCub->getNumPoints();
537 basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
538 basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
539 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
540 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
541
542 cub_degree = basisCubDegree + targetCubDegree;
543 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree);
544 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints + elemTargetCub->getNumPoints());
545 numTargetEvalPoints += elemTargetCub->getNumPoints();
546 targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
547 targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
548 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
549 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
550 }
551 }
552}
553
554template<typename DeviceType, typename ValueType>
555template<typename BasisPtrType>
557 const ordinal_type targetCubDegree) {
558 const auto cellTopo = cellBasis->getBaseCellTopology();
559 ordinal_type dim = cellTopo.getDimension();
560 numBasisEvalPoints = 0;
561 numBasisDerivEvalPoints = 0;
562 numTargetEvalPoints = 0;
563 numTargetDerivEvalPoints = 0;
564
565 basisPointsRange = range_tag("basisPointsRange", 4,maxSubCellsCount);
566 targetPointsRange = range_tag("targetPointsRange", 4,maxSubCellsCount);
567 basisDerivPointsRange = range_tag("basisDerivPointsRange", 4,maxSubCellsCount);
568 targetDerivPointsRange = range_tag("targetDerivPointsRange", 4,maxSubCellsCount);
569 subCellTopologyKey = key_tag("subCellTopologyKey",4,maxSubCellsCount);
570
571 ordinal_type basisCubDegree = cellBasis->getDegree();
572 DefaultCubatureFactory cub_factory;
573 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
574
575 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints =0;
576
577 ordinal_type cub_degree = 2*basisCubDegree;
578 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree);
579 basisPointsRange(dim,0) = range_type(0, elemBasisCub->getNumPoints());
580 numBasisEvalPoints += elemBasisCub->getNumPoints();
581 maxNumBasisEvalPoints = elemBasisCub->getNumPoints();
582 basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
583 basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
584 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
585
586 cub_degree = basisCubDegree + targetCubDegree;
587 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree);
588 targetPointsRange(dim,0) = range_type(0, elemTargetCub->getNumPoints());
589 numTargetEvalPoints += elemTargetCub->getNumPoints();
590 maxNumTargetEvalPoints = elemTargetCub->getNumPoints();
591 targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
592 targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
593 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
594}
595
596}
597}
598#endif
599
600
601
602
603
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
Header file for the Intrepid2::Basis_HCURL_HEX_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_TET_In_FEM class.
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
Header file for the Intrepid2::Basis_HGRAD_TRI_Cn_FEM class.
Header function for Intrepid2::Util class and other utility functions.
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
Implementation of the default H(grad)-compatible FEM basis of degree n on Quadrilateral cell Implemen...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
static void getReferenceVertex(Kokkos::DynRankView< cellVertexValueType, cellVertexProperties... > cellVertex, const shards::CellTopology cell, const ordinal_type vertexOrd)
Retrieves the Cartesian coordinates of a reference cell vertex.
A factory class that generates specific instances of cubatures.
static Teuchos::RCP< Cubature< DeviceType, pointValueType, weightValueType > > create(unsigned topologyKey, const std::vector< ordinal_type > &degree, const EPolyType polytype=POLYTYPE_MAX)
Factory method.
void createHCurlProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetCurlCubDegre)
Initialize the ProjectionStruct for HCURL projections.
void createHDivProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetDivCubDegre)
Initialize the ProjectionStruct for HDIV projections.
void createL2ProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for L2 projections.
void createHVolProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for HVOL (local-L2) projection.
void createHGradProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetGradCubDegre)
Initialize the ProjectionStruct for HGRAD projections.