Intrepid2
Intrepid2_ProjectionStructDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
14#ifndef __INTREPID2_PROJECTIONSTRUCTDEF_HPP__
15#define __INTREPID2_PROJECTIONSTRUCTDEF_HPP__
16
17#include "Intrepid2_Utils.hpp"
18#include "Shards_CellTopology.hpp"
19#include "Shards_BasicTopologies.hpp"
20
21#include "Intrepid2_HGRAD_QUAD_Cn_FEM.hpp"
22#include "Intrepid2_HGRAD_TRI_Cn_FEM.hpp"
23#include "Intrepid2_HCURL_HEX_In_FEM.hpp"
24#include "Intrepid2_HCURL_TET_In_FEM.hpp"
25
26#include "Intrepid2_DefaultCubatureFactory.hpp"
27
28
29namespace Intrepid2 {
30
31template<typename DeviceType, typename ValueType>
32template<typename BasisPtrType>
33void ProjectionStruct<DeviceType,ValueType>::createL2ProjectionStruct(const BasisPtrType cellBasis,
34 const ordinal_type targetCubDegree) {
35 const auto cellTopo = cellBasis->getBaseCellTopology();
36 std::string name = cellBasis->getName();
37 ordinal_type dim = cellTopo.getDimension();
38 numBasisEvalPoints = 0;
39 numBasisDerivEvalPoints = 0;
40 numTargetEvalPoints = 0;
41 numTargetDerivEvalPoints = 0;
42 const ordinal_type edgeDim = 1;
43 const ordinal_type faceDim = 2;
44
45 ordinal_type basisCubDegree = cellBasis->getDegree();
46 ordinal_type edgeBasisCubDegree, faceBasisCubDegree;
47
48 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
49 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
50 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
51
52 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
53
54 INTREPID2_TEST_FOR_ABORT( (numVertices > maxSubCellsCount) || (numFaces > maxSubCellsCount) || (numEdges > maxSubCellsCount),
55 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with this cub cells count");
56
57
58 numBasisEvalPoints += numVertices;
59 numTargetEvalPoints += numVertices;
60
61 basisPointsRange = range_tag("basisPointsRange", 4,maxSubCellsCount);
62 targetPointsRange = range_tag("targetPointsRange", 4,maxSubCellsCount);
63 basisDerivPointsRange = range_tag("basisDerivPointsRange", 4,maxSubCellsCount);
64 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
65 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
66
67 maxNumBasisEvalPoints = numVertices; maxNumTargetEvalPoints = numVertices;
68
69 // The set of the eval points on the reference vertex contains only point (0.0).
70 // Not very useful, updating these for completeness
71 for(ordinal_type iv=0; iv<numVertices; ++iv) {
72 basisPointsRange(0,iv) = range_type(iv, iv+1);
73 basisCubPoints[0][iv] = host_view_type("basisCubPoints",1,1);
74 targetPointsRange(0,iv) = range_type(iv, iv+1);
75 targetCubPoints[0][iv] = host_view_type("targetCubPoints",1,1);
76 }
77
78 if (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HCURL) {
79 edgeBasisCubDegree = basisCubDegree - 1;
80 faceBasisCubDegree = basisCubDegree;
81 }
82 else if (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HDIV) {
83 edgeBasisCubDegree = basisCubDegree - 1;
84 faceBasisCubDegree = basisCubDegree -1;
85 }
86 else {
87 edgeBasisCubDegree = basisCubDegree;
88 faceBasisCubDegree = basisCubDegree;
89 }
90
91 DefaultCubatureFactory cub_factory;
92 for(ordinal_type ie=0; ie<numEdges; ++ie) {
93 ordinal_type cub_degree = 2*edgeBasisCubDegree;
94 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
95 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
96 basisPointsRange(edgeDim,ie) = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
97 numBasisEvalPoints += edgeBasisCub->getNumPoints();
98 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, edgeBasisCub->getNumPoints());
99 basisCubPoints[edgeDim][ie] = host_view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
100 basisCubWeights[edgeDim][ie] = host_view_type("basisCubWeights",edgeBasisCub->getNumPoints());
101 edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
102
103 cub_degree = edgeBasisCubDegree + targetCubDegree;
104 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
105 targetPointsRange(edgeDim,ie) = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
106 numTargetEvalPoints += edgeTargetCub->getNumPoints();
107 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, edgeTargetCub->getNumPoints());
108 targetCubPoints[edgeDim][ie] = host_view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
109 targetCubWeights[edgeDim][ie] = host_view_type("targetCubWeights",edgeTargetCub->getNumPoints());
110 edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
111 }
112
113 for(ordinal_type iface=0; iface<numFaces; ++iface) {
114 ordinal_type cub_degree = 2*faceBasisCubDegree;
115 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
116 auto faceBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
117 basisPointsRange(faceDim,iface) = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
118 numBasisEvalPoints += faceBasisCub->getNumPoints();
119 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, faceBasisCub->getNumPoints());
120 basisCubPoints[faceDim][iface] = host_view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
121 basisCubWeights[faceDim][iface] = host_view_type("basisCubWeights",faceBasisCub->getNumPoints());
122 faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
123
124 cub_degree = faceBasisCubDegree + targetCubDegree;
125 auto faceTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
126 targetPointsRange(faceDim,iface) = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
127 numTargetEvalPoints += faceTargetCub->getNumPoints();
128 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, faceTargetCub->getNumPoints());
129 targetCubPoints[faceDim][iface] = host_view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
130 targetCubWeights[faceDim][iface] = host_view_type("targetCubWeights",faceTargetCub->getNumPoints());
131 faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
132 }
133 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
134 if(hasCellDofs) {
135 ordinal_type cub_degree = 2*basisCubDegree;
136 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
137 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
138 numBasisEvalPoints += elemBasisCub->getNumPoints();
139 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
140 basisCubPoints[dim][0] = host_view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
141 basisCubWeights[dim][0] = host_view_type("basisCubWeights",elemBasisCub->getNumPoints());
142 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
143
144 cub_degree = basisCubDegree + targetCubDegree;
145 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
146 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
147 numTargetEvalPoints += elemTargetCub->getNumPoints();
148 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
149 targetCubPoints[dim][0] = host_view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
150 targetCubWeights[dim][0] = host_view_type("targetCubWeights",elemTargetCub->getNumPoints());
151 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
152
153
154 }
155
156 allBasisEPoints = view_type("allBasisPoints", numBasisEvalPoints, dim);
157 allTargetEPoints = view_type("allTargetPoints", numTargetEvalPoints, dim);
158 allBasisDerivEPoints = view_type("allBasisDerivPoints", numBasisDerivEvalPoints, dim);
159 allTargetDerivEPoints = view_type("allTargetDerivPoints", numTargetDerivEvalPoints, dim);
160
161 if(numVertices>0) {
162 for(ordinal_type iv=0; iv<numVertices; ++iv) {
163 CellTools<DeviceType>::getReferenceVertex(Kokkos::subview(allBasisEPoints, iv, Kokkos::ALL()), cellTopo, iv);
164 CellTools<DeviceType>::getReferenceVertex(Kokkos::subview(allTargetEPoints, iv, Kokkos::ALL()), cellTopo, iv);
165 }
166 }
167
168 if(numEdges>0) {
169 auto subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
170 for(ordinal_type ie=0; ie<numEdges; ++ie) {
171 auto edgeBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[edgeDim][ie]);
172 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisEPoints,basisPointsRange(edgeDim,ie),Kokkos::ALL()), edgeBasisEPoints, subcellParamEdge, ie);
173
174 auto edgeTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[edgeDim][ie]);
175 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetEPoints,targetPointsRange(edgeDim,ie),Kokkos::ALL()), edgeTargetEPoints, subcellParamEdge, ie);
176 }
177 }
178
179 if(numFaces>0) {
180 auto subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
181 for(ordinal_type iface=0; iface<numFaces; ++iface) {
182 auto faceBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[faceDim][iface]);
183 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisEPoints,basisPointsRange(faceDim,iface),Kokkos::ALL()), faceBasisEPoints, subcellParamFace, iface);
184
185 auto faceTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[faceDim][iface]);
186 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetEPoints,targetPointsRange(faceDim,iface),Kokkos::ALL()), faceTargetEPoints, subcellParamFace, iface);
187 }
188 }
189
190 if(hasCellDofs) {
191 auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[dim][0]);
192 Kokkos::deep_copy(Kokkos::subview(allBasisEPoints,basisPointsRange(dim,0), Kokkos::ALL()), cellBasisEPoints);
193
194 auto cellTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[dim][0]);
195 Kokkos::deep_copy(Kokkos::subview(allTargetEPoints,targetPointsRange(dim,0),Kokkos::ALL()), cellTargetEPoints);
196 }
197}
198
199template<typename DeviceType, typename ValueType>
200template<typename BasisPtrType>
201void ProjectionStruct<DeviceType,ValueType>::createHGradProjectionStruct(const BasisPtrType cellBasis,
202 const ordinal_type targetCubDegree,
203 const ordinal_type targetGradCubDegree) {
204 const auto cellTopo = cellBasis->getBaseCellTopology();
205 std::string name = cellBasis->getName();
206 ordinal_type dim = cellTopo.getDimension();
207 numBasisEvalPoints = 0;
208 numBasisDerivEvalPoints = 0;
209 numTargetEvalPoints = 0;
210 numTargetDerivEvalPoints = 0;
211 const ordinal_type edgeDim = 1;
212 const ordinal_type faceDim = 2;
213
214 ordinal_type basisCubDegree = cellBasis->getDegree();
215 ordinal_type edgeBasisCubDegree = basisCubDegree;
216 ordinal_type faceBasisCubDegree = basisCubDegree;
217
218 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
219 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount(): 0;
220 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
221
222 INTREPID2_TEST_FOR_ABORT( (numFaces > maxSubCellsCount) || (numEdges > maxSubCellsCount),
223 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with this cub cells count");
224
225
226 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
227
228 maxNumBasisEvalPoints = numVertices; maxNumTargetEvalPoints = numVertices;
229 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
230
231 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
232 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
233 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
234 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
235 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
236
237 numBasisEvalPoints += numVertices;
238 numTargetEvalPoints += numVertices;
239
240 // The set of the eval points on the reference vertex contains only point (0.0).
241 // Not very useful, updating these for completeness
242 for(ordinal_type iv=0; iv<numVertices; ++iv) {
243 basisPointsRange(0,iv) = range_type(iv, iv+1);
244 basisCubPoints[0][iv] = host_view_type("basisCubPoints",1,1);
245 targetPointsRange(0,iv) = range_type(iv, iv+1);
246 targetCubPoints[0][iv] = host_view_type("targetCubPoints",1,1);
247 }
248
249 DefaultCubatureFactory cub_factory;
250 for(ordinal_type ie=0; ie<numEdges; ++ie) {
251 ordinal_type cub_degree = 2*edgeBasisCubDegree;
252 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
253 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
254 basisDerivPointsRange(edgeDim,ie) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+edgeBasisCub->getNumPoints());
255 numBasisDerivEvalPoints += edgeBasisCub->getNumPoints();
256 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, edgeBasisCub->getNumPoints());
257 basisDerivCubPoints[edgeDim][ie] = host_view_type("basisDerivCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
258 basisDerivCubWeights[edgeDim][ie] = host_view_type("basisDerivCubWeights",edgeBasisCub->getNumPoints());
259 edgeBasisCub->getCubature(basisDerivCubPoints[edgeDim][ie], basisDerivCubWeights[edgeDim][ie]);
260
261 cub_degree = edgeBasisCubDegree + targetGradCubDegree;
262 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
263 targetDerivPointsRange(edgeDim,ie) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+edgeTargetCub->getNumPoints());
264 numTargetDerivEvalPoints += edgeTargetCub->getNumPoints();
265 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, edgeTargetCub->getNumPoints());
266 targetDerivCubPoints[edgeDim][ie] = host_view_type("targetDerivCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
267 targetDerivCubWeights[edgeDim][ie] = host_view_type("targetDerivCubWeights",edgeTargetCub->getNumPoints());
268 edgeTargetCub->getCubature(targetDerivCubPoints[edgeDim][ie], targetDerivCubWeights[edgeDim][ie]);
269 }
270
271 for(ordinal_type iface=0; iface<numFaces; ++iface) {
272 ordinal_type cub_degree = 2*faceBasisCubDegree;
273 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
274 auto faceBasisGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
275 basisDerivPointsRange(faceDim,iface) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisGradCub->getNumPoints());
276 numBasisDerivEvalPoints += faceBasisGradCub->getNumPoints();
277 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, faceBasisGradCub->getNumPoints());
278 basisDerivCubPoints[faceDim][iface] = host_view_type("basisDerivCubPoints",faceBasisGradCub->getNumPoints(),faceDim);
279 basisDerivCubWeights[faceDim][iface] = host_view_type("basisDerivCubWeights",faceBasisGradCub->getNumPoints());
280 faceBasisGradCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
281
282 cub_degree = faceBasisCubDegree + targetGradCubDegree;
283 auto faceTargetDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
284 targetDerivPointsRange(faceDim,iface) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
285 numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
286 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, faceTargetDerivCub->getNumPoints());
287 targetDerivCubPoints[faceDim][iface] = host_view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
288 targetDerivCubWeights[faceDim][iface] = host_view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
289 faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
290 }
291 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
292 if(hasCellDofs) {
293 ordinal_type cub_degree = 2*basisCubDegree;
294 auto elemBasisGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
295 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisGradCub->getNumPoints());
296 numBasisDerivEvalPoints += elemBasisGradCub->getNumPoints();
297 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisGradCub->getNumPoints());
298 basisDerivCubPoints[dim][0] = host_view_type("basisDerivCubPoints",elemBasisGradCub->getNumPoints(),dim);
299 basisDerivCubWeights[dim][0] = host_view_type("basisDerivCubWeights",elemBasisGradCub->getNumPoints());
300 elemBasisGradCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
301
302 cub_degree = basisCubDegree + targetGradCubDegree;
303 auto elemTargetGradCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
304 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetGradCub->getNumPoints());
305 numTargetDerivEvalPoints += elemTargetGradCub->getNumPoints();
306 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetGradCub->getNumPoints());
307 targetDerivCubPoints[dim][0] = host_view_type("targetDerivCubPoints",elemTargetGradCub->getNumPoints(),dim);
308 targetDerivCubWeights[dim][0] = host_view_type("targetDerivCubWeights",elemTargetGradCub->getNumPoints());
309 elemTargetGradCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
310 }
311
312 allBasisEPoints = view_type("allBasisPoints", numBasisEvalPoints, dim);
313 allTargetEPoints = view_type("allTargetPoints", numTargetEvalPoints, dim);
314 allBasisDerivEPoints = view_type("allBasisDerivPoints", numBasisDerivEvalPoints, dim);
315 allTargetDerivEPoints = view_type("allTargetDerivPoints", numTargetDerivEvalPoints, dim);
316
317 if(numVertices>0) {
318 for(ordinal_type iv=0; iv<numVertices; ++iv) {
319 CellTools<DeviceType>::getReferenceVertex(Kokkos::subview(allBasisEPoints, iv, Kokkos::ALL()), cellTopo, iv);
320 CellTools<DeviceType>::getReferenceVertex(Kokkos::subview(allTargetEPoints, iv, Kokkos::ALL()), cellTopo, iv);
321 }
322 }
323
324 if(numEdges>0) {
325 auto subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
326 for(ordinal_type ie=0; ie<numEdges; ++ie) {
327 auto edgeBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[edgeDim][ie]);
328 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(edgeDim,ie),Kokkos::ALL()), edgeBasisDerivEPoints, subcellParamEdge, ie);
329
330 auto edgeTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[edgeDim][ie]);
331 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(edgeDim,ie),Kokkos::ALL()), edgeTargetDerivEPoints, subcellParamEdge, ie);
332 }
333 }
334
335 if(numFaces>0) {
336 auto subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
337 for(ordinal_type iface=0; iface<numFaces; ++iface) {
338 auto faceBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[faceDim][iface]);
339 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(faceDim,iface),Kokkos::ALL()), faceBasisDerivEPoints, subcellParamFace, iface);
340
341 auto faceTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[faceDim][iface]);
342 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(faceDim,iface),Kokkos::ALL()), faceTargetDerivEPoints, subcellParamFace, iface);
343 }
344 }
345
346 if(hasCellDofs) {
347 auto cellBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[dim][0]);
348 Kokkos::deep_copy(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(dim,0), Kokkos::ALL()), cellBasisDerivEPoints);
349
350 auto cellTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[dim][0]);
351 Kokkos::deep_copy(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(dim,0),Kokkos::ALL()), cellTargetDerivEPoints);
352 }
353}
354
355template<typename DeviceType, typename ValueType>
356template<typename BasisPtrType>
357void ProjectionStruct<DeviceType,ValueType>::createHCurlProjectionStruct(const BasisPtrType cellBasis,
358 const ordinal_type targetCubDegree,
359 const ordinal_type targetCurlCubDegre) {
360 const auto cellTopo = cellBasis->getBaseCellTopology();
361 std::string name = cellBasis->getName();
362 ordinal_type dim = cellTopo.getDimension();
363 numBasisEvalPoints = 0;
364 numBasisDerivEvalPoints = 0;
365 numTargetEvalPoints = 0;
366 numTargetDerivEvalPoints = 0;
367 const ordinal_type edgeDim = 1;
368 const ordinal_type faceDim = 2;
369
370 ordinal_type basisCubDegree = cellBasis->getDegree();
371 ordinal_type edgeBasisCubDegree = basisCubDegree - 1;
372 ordinal_type faceBasisCubDegree = basisCubDegree;
373 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
374 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
375 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
376
377 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints = 0;
378 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
379
380 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
381 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
382 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
383 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
384 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
385
386 DefaultCubatureFactory cub_factory;
387 for(ordinal_type ie=0; ie<numEdges; ++ie) {
388 ordinal_type cub_degree = 2*edgeBasisCubDegree;
389 subCellTopologyKey(edgeDim,ie) = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
390 auto edgeBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(edgeDim,ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
391 basisPointsRange(edgeDim,ie) = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
392 numBasisEvalPoints += edgeBasisCub->getNumPoints();
393 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, edgeBasisCub->getNumPoints());
394 basisCubPoints[edgeDim][ie] = host_view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
395 basisCubWeights[edgeDim][ie] = host_view_type("basisCubWeights",edgeBasisCub->getNumPoints());
396 edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
397
398 cub_degree = edgeBasisCubDegree + targetCubDegree;
399 auto edgeTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(edgeDim,ie), cub_degree, EPolyType::POLYTYPE_MAX, true);
400 targetPointsRange(edgeDim,ie) = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
401 numTargetEvalPoints += edgeTargetCub->getNumPoints();
402 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, edgeTargetCub->getNumPoints());
403 targetCubPoints[edgeDim][ie] = host_view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
404 targetCubWeights[edgeDim][ie] = host_view_type("targetCubWeights",edgeTargetCub->getNumPoints());
405 edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
406 }
407
408 for(ordinal_type iface=0; iface<numFaces; ++iface) {
409 ordinal_type cub_degree = 2*faceBasisCubDegree;
410 subCellTopologyKey(faceDim,iface) = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
411 auto faceBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
412 basisPointsRange(faceDim,iface) = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
413 numBasisEvalPoints += faceBasisCub->getNumPoints();
414 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, faceBasisCub->getNumPoints());
415 basisCubPoints[faceDim][iface] = host_view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
416 basisCubWeights[faceDim][iface] = host_view_type("basisCubWeights",faceBasisCub->getNumPoints());
417 faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
418
419 auto faceBasisDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
420 basisDerivPointsRange(faceDim,iface) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisCub->getNumPoints());
421 numBasisDerivEvalPoints += faceBasisCub->getNumPoints();
422 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, faceBasisCub->getNumPoints());
423 basisDerivCubPoints[faceDim][iface] = host_view_type("basisDerivCubPoints",faceBasisCub->getNumPoints(),faceDim);
424 basisDerivCubWeights[faceDim][iface] = host_view_type("basisDerivCubWeights",faceBasisCub->getNumPoints());
425 faceBasisCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
426
427 cub_degree = faceBasisCubDegree + targetCubDegree;
428 auto faceTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
429 targetPointsRange(faceDim,iface) = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
430 numTargetEvalPoints += faceTargetCub->getNumPoints();
431 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, faceTargetCub->getNumPoints());
432 targetCubPoints[faceDim][iface] = host_view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
433 targetCubWeights[faceDim][iface] = host_view_type("targetCubWeights",faceTargetCub->getNumPoints());
434 faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
435
436 cub_degree = faceBasisCubDegree + targetCurlCubDegre;
437 auto faceTargetDerivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(faceDim,iface), cub_degree, EPolyType::POLYTYPE_MAX, true);
438 targetDerivPointsRange(faceDim,iface) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
439 numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
440 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, faceTargetDerivCub->getNumPoints());
441 targetDerivCubPoints[faceDim][iface] = host_view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
442 targetDerivCubWeights[faceDim][iface] = host_view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
443 faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
444 }
445
446 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
447 if(hasCellDofs) {
448 ordinal_type cub_degree = 2*basisCubDegree;
449 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
450 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
451 numBasisEvalPoints += elemBasisCub->getNumPoints();
452 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
453 basisCubPoints[dim][0] = host_view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
454 basisCubWeights[dim][0] = host_view_type("basisCubWeights",elemBasisCub->getNumPoints());
455 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
456
457 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisCub->getNumPoints());
458 numBasisDerivEvalPoints += elemBasisCub->getNumPoints();
459 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisCub->getNumPoints());
460 basisDerivCubPoints[dim][0] = host_view_type("basisDerivCubPoints",elemBasisCub->getNumPoints(),dim);
461 basisDerivCubWeights[dim][0] = host_view_type("basisDerivCubWeights",elemBasisCub->getNumPoints());
462 elemBasisCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
463
464 cub_degree = basisCubDegree + targetCubDegree;
465 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
466 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
467 numTargetEvalPoints += elemTargetCub->getNumPoints();
468 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
469 targetCubPoints[dim][0] = host_view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
470 targetCubWeights[dim][0] = host_view_type("targetCubWeights",elemTargetCub->getNumPoints());
471 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
472
473 cub_degree = basisCubDegree + targetCurlCubDegre;
474 auto elemTargetCurlCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
475 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetCurlCub->getNumPoints());
476 numTargetDerivEvalPoints += elemTargetCurlCub->getNumPoints();
477 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetCurlCub->getNumPoints());
478 targetDerivCubPoints[dim][0] = host_view_type("targetDerivCubPoints",elemTargetCurlCub->getNumPoints(),dim);
479 targetDerivCubWeights[dim][0] = host_view_type("targetDerivCubWeights",elemTargetCurlCub->getNumPoints());
480 elemTargetCurlCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
481 }
482
483 allBasisEPoints = view_type("allBasisPoints", numBasisEvalPoints, dim);
484 allTargetEPoints = view_type("allTargetPoints", numTargetEvalPoints, dim);
485 allBasisDerivEPoints = view_type("allBasisDerivPoints", numBasisDerivEvalPoints, dim);
486 allTargetDerivEPoints = view_type("allTargetDerivPoints", numTargetDerivEvalPoints, dim);
487
488 auto subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
489 for(ordinal_type ie=0; ie<numEdges; ++ie) {
490 auto edgeBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[edgeDim][ie]);
491 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisEPoints,basisPointsRange(edgeDim,ie),Kokkos::ALL()), edgeBasisEPoints, subcellParamEdge, ie);
492
493 auto edgeTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[edgeDim][ie]);
494 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetEPoints,targetPointsRange(edgeDim,ie),Kokkos::ALL()), edgeTargetEPoints, subcellParamEdge, ie);
495 }
496
497 if(numFaces>0) {
498 auto subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
499 for(ordinal_type iface=0; iface<numFaces; ++iface) {
500 auto faceBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[faceDim][iface]);
501 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisEPoints,basisPointsRange(faceDim,iface),Kokkos::ALL()), faceBasisEPoints, subcellParamFace, iface);
502
503 auto faceBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[faceDim][iface]);
504 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(faceDim,iface),Kokkos::ALL()), faceBasisDerivEPoints, subcellParamFace, iface);
505
506 auto faceTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[faceDim][iface]);
507 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetEPoints,targetPointsRange(faceDim,iface),Kokkos::ALL()), faceTargetEPoints, subcellParamFace, iface);
508
509 auto faceTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[faceDim][iface]);
510 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(faceDim,iface),Kokkos::ALL()), faceTargetDerivEPoints, subcellParamFace, iface);
511 }
512 }
513
514 if(hasCellDofs) {
515 auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[dim][0]);
516 Kokkos::deep_copy(Kokkos::subview(allBasisEPoints,basisPointsRange(dim,0), Kokkos::ALL()), cellBasisEPoints);
517
518 auto cellBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[dim][0]);
519 Kokkos::deep_copy(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(dim,0), Kokkos::ALL()), cellBasisDerivEPoints);
520
521 auto cellTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[dim][0]);
522 Kokkos::deep_copy(Kokkos::subview(allTargetEPoints,targetPointsRange(dim,0),Kokkos::ALL()), cellTargetEPoints);
523
524 auto cellTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[dim][0]);
525 Kokkos::deep_copy(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(dim,0),Kokkos::ALL()), cellTargetDerivEPoints);
526 }
527}
528
529template<typename DeviceType, typename ValueType>
530template<typename BasisPtrType>
531void ProjectionStruct<DeviceType,ValueType>::createHDivProjectionStruct(const BasisPtrType cellBasis,
532 const ordinal_type targetCubDegree,
533 const ordinal_type targetDivCubDegre) {
534 const auto cellTopo = cellBasis->getBaseCellTopology();
535 std::string name = cellBasis->getName();
536 ordinal_type dim = cellTopo.getDimension();
537 numBasisEvalPoints = 0;
538 numBasisDerivEvalPoints = 0;
539 numTargetEvalPoints = 0;
540 numTargetDerivEvalPoints = 0;
541 const ordinal_type sideDim = dim - 1;
542
543 ordinal_type basisCubDegree = cellBasis->getDegree();
544 ordinal_type sideBasisCubDegree = basisCubDegree - 1;
545 ordinal_type numSides = cellTopo.getSideCount()*ordinal_type(cellBasis->getDofCount(sideDim, 0) > 0);
546 ordinal_type hasCellDofs = (cellBasis->getDofCount(dim, 0) > 0);
547
548 INTREPID2_TEST_FOR_ABORT( numSides > maxSubCellsCount,
549 ">>> ERROR (Intrepid2::ProjectionStruct:createHDivProjectionStruct, Projections do not support a cell topology with so many sides");
550
551 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints = 0;
552 maxNumBasisDerivEvalPoints = 0; maxNumTargetDerivEvalPoints = 0;
553
554 basisPointsRange = range_tag("basisPointsRange", numberSubCellDims,maxSubCellsCount);
555 targetPointsRange = range_tag("targetPointsRange", numberSubCellDims,maxSubCellsCount);
556 basisDerivPointsRange = range_tag("basisDerivPointsRange", numberSubCellDims,maxSubCellsCount);
557 targetDerivPointsRange = range_tag("targetDerivPointsRange", numberSubCellDims,maxSubCellsCount);
558 subCellTopologyKey = key_tag("subCellTopologyKey",numberSubCellDims,maxSubCellsCount);
559
560 Basis<HostDeviceType,ValueType,ValueType> *hcurlBasis = NULL;
561 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
562 hcurlBasis = new Basis_HCURL_HEX_In_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
563 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
564 hcurlBasis = new Basis_HCURL_TET_In_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
565 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Wedge<6> >()->key)
566 hcurlBasis = new typename DerivedNodalBasisFamily<HostDeviceType,ValueType,ValueType>::HCURL_WEDGE(cellBasis->getDegree());
567// TODO: uncomment the next two lines once H(curl) pyramid implemented
568// else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Pyramid<5> >()->key)
569// hcurlBasis = new typename DerivedNodalBasisFamily<DeviceType,ValueType,ValueType>::HCURL_PYR(cellBasis->getDegree());
570 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
571 hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
572 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
573 hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM<HostDeviceType,ValueType,ValueType>(cellBasis->getDegree());
574 else {
575 std::stringstream ss;
576 ss << ">>> ERROR (Intrepid2::ProjectionTools::createHDivProjectionStruct): "
577 << "Method not implemented for basis " << name;
578 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
579 }
580
581 bool haveHCurlConstraint = (hcurlBasis != NULL) && (hcurlBasis->getDofCount(dim,0) > 0);
582 if(hcurlBasis != NULL) delete hcurlBasis;
583
584 DefaultCubatureFactory cub_factory;
585
586 for(ordinal_type is=0; is<numSides; ++is) {
587 ordinal_type cub_degree = 2*sideBasisCubDegree;
588 subCellTopologyKey(sideDim,is) = cellBasis->getBaseCellTopology().getKey(sideDim, is);
589 auto sideBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(sideDim,is), cub_degree, EPolyType::POLYTYPE_MAX, true);
590 basisPointsRange(sideDim,is) = range_type(numBasisEvalPoints, numBasisEvalPoints+sideBasisCub->getNumPoints());
591 numBasisEvalPoints += sideBasisCub->getNumPoints();
592 basisCubPoints[sideDim][is] = host_view_type("basisCubPoints",sideBasisCub->getNumPoints(),sideDim);
593 basisCubWeights[sideDim][is] = host_view_type("basisCubWeights",sideBasisCub->getNumPoints());
594 sideBasisCub->getCubature(basisCubPoints[sideDim][is], basisCubWeights[sideDim][is]);
595 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, sideBasisCub->getNumPoints());
596
597 cub_degree = sideBasisCubDegree + targetCubDegree;
598 auto sideTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(sideDim,is), cub_degree, EPolyType::POLYTYPE_MAX, true);
599 targetPointsRange(sideDim,is) = range_type(numTargetEvalPoints, numTargetEvalPoints+sideTargetCub->getNumPoints());
600 numTargetEvalPoints += sideTargetCub->getNumPoints();
601 targetCubPoints[sideDim][is] = host_view_type("targetCubPoints",sideTargetCub->getNumPoints(),sideDim);
602 targetCubWeights[sideDim][is] = host_view_type("targetCubWeights",sideTargetCub->getNumPoints());
603 sideTargetCub->getCubature(targetCubPoints[sideDim][is], targetCubWeights[sideDim][is]);
604 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, sideTargetCub->getNumPoints());
605 }
606
607 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
608 if(hasCellDofs) {
609 ordinal_type cub_degree = 2*basisCubDegree - 1;
610 auto elemBasisDivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
611 basisDerivPointsRange(dim,0) = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisDivCub->getNumPoints());
612 numBasisDerivEvalPoints += elemBasisDivCub->getNumPoints();
613 basisDerivCubPoints[dim][0] = host_view_type("basisDerivCubPoints",elemBasisDivCub->getNumPoints(),dim);
614 basisDerivCubWeights[dim][0] = host_view_type("basisDerivCubWeights",elemBasisDivCub->getNumPoints());
615 elemBasisDivCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
616 maxNumBasisDerivEvalPoints = std::max(maxNumBasisDerivEvalPoints, elemBasisDivCub->getNumPoints());
617
618 cub_degree = basisCubDegree - 1 + targetDivCubDegre;
619 auto elemTargetDivCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
620 targetDerivPointsRange(dim,0) = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetDivCub->getNumPoints());
621 numTargetDerivEvalPoints += elemTargetDivCub->getNumPoints();
622 targetDerivCubPoints[dim][0] = host_view_type("targetDerivCubPoints",elemTargetDivCub->getNumPoints(),dim);
623 targetDerivCubWeights[dim][0] = host_view_type("targetDerivCubWeights",elemTargetDivCub->getNumPoints());
624 elemTargetDivCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
625 maxNumTargetDerivEvalPoints = std::max(maxNumTargetDerivEvalPoints, elemTargetDivCub->getNumPoints());
626
627 if(haveHCurlConstraint)
628 {
629 cub_degree = 2*basisCubDegree;
630 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
631 basisPointsRange(dim,0) = range_type(numBasisEvalPoints, numBasisEvalPoints + elemBasisCub->getNumPoints());
632 numBasisEvalPoints += elemBasisCub->getNumPoints();
633 basisCubPoints[dim][0] = host_view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
634 basisCubWeights[dim][0] = host_view_type("basisCubWeights",elemBasisCub->getNumPoints());
635 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
636 maxNumBasisEvalPoints = std::max(maxNumBasisEvalPoints, elemBasisCub->getNumPoints());
637
638 cub_degree = basisCubDegree + targetCubDegree;
639 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(subCellTopologyKey(dim,0), cub_degree, EPolyType::POLYTYPE_MAX, true);
640 targetPointsRange(dim,0) = range_type(numTargetEvalPoints, numTargetEvalPoints + elemTargetCub->getNumPoints());
641 numTargetEvalPoints += elemTargetCub->getNumPoints();
642 targetCubPoints[dim][0] = host_view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
643 targetCubWeights[dim][0] = host_view_type("targetCubWeights",elemTargetCub->getNumPoints());
644 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
645 maxNumTargetEvalPoints = std::max(maxNumTargetEvalPoints, elemTargetCub->getNumPoints());
646 }
647 }
648
649 allBasisEPoints = view_type("allBasisPoints", numBasisEvalPoints, dim);
650 allTargetEPoints = view_type("allTargetPoints", numTargetEvalPoints, dim);
651 allBasisDerivEPoints = view_type("allBasisDerivPoints", numBasisDerivEvalPoints, dim);
652 allTargetDerivEPoints = view_type("allTargetDerivPoints", numTargetDerivEvalPoints, dim);
653
654
655 auto subcellParamSide = RefSubcellParametrization<DeviceType>::get(sideDim, cellTopo.getKey());
656 for(ordinal_type is=0; is<numSides; ++is) {
657 auto sideBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[sideDim][is]);
658 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allBasisEPoints,basisPointsRange(sideDim,is),Kokkos::ALL()), sideBasisEPoints, subcellParamSide, is);
659
660 auto sideTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[sideDim][is]);
661 CellTools<DeviceType>::mapToReferenceSubcell(Kokkos::subview(allTargetEPoints,targetPointsRange(sideDim,is),Kokkos::ALL()), sideTargetEPoints, subcellParamSide, is);
662 }
663
664 if(hasCellDofs) {
665 if(haveHCurlConstraint) {
666 auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[dim][0]);
667 Kokkos::deep_copy(Kokkos::subview(allBasisEPoints,basisPointsRange(dim,0), Kokkos::ALL()), cellBasisEPoints);
668
669 auto cellTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[dim][0]);
670 Kokkos::deep_copy(Kokkos::subview(allTargetEPoints,targetPointsRange(dim,0),Kokkos::ALL()), cellTargetEPoints);
671 }
672
673 auto cellBasisDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisDerivCubPoints[dim][0]);
674 Kokkos::deep_copy(Kokkos::subview(allBasisDerivEPoints,basisDerivPointsRange(dim,0), Kokkos::ALL()), cellBasisDerivEPoints);
675
676 auto cellTargetDerivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetDerivCubPoints[dim][0]);
677 Kokkos::deep_copy(Kokkos::subview(allTargetDerivEPoints,targetDerivPointsRange(dim,0),Kokkos::ALL()), cellTargetDerivEPoints);
678 }
679}
680
681template<typename DeviceType, typename ValueType>
682template<typename BasisPtrType>
683void ProjectionStruct<DeviceType,ValueType>::createHVolProjectionStruct(const BasisPtrType cellBasis,
684 const ordinal_type targetCubDegree) {
685 const auto cellTopo = cellBasis->getBaseCellTopology();
686 ordinal_type dim = cellTopo.getDimension();
687 numBasisEvalPoints = 0;
688 numBasisDerivEvalPoints = 0;
689 numTargetEvalPoints = 0;
690 numTargetDerivEvalPoints = 0;
691
692 basisPointsRange = range_tag("basisPointsRange", 4,maxSubCellsCount);
693 targetPointsRange = range_tag("targetPointsRange", 4,maxSubCellsCount);
694 basisDerivPointsRange = range_tag("basisDerivPointsRange", 4,maxSubCellsCount);
695 targetDerivPointsRange = range_tag("targetDerivPointsRange", 4,maxSubCellsCount);
696 subCellTopologyKey = key_tag("subCellTopologyKey",4,maxSubCellsCount);
697
698 ordinal_type basisCubDegree = cellBasis->getDegree();
699 DefaultCubatureFactory cub_factory;
700 subCellTopologyKey(dim,0) = cellBasis->getBaseCellTopology().getBaseKey();
701
702 maxNumBasisEvalPoints = 0; maxNumTargetEvalPoints =0;
703
704 ordinal_type cub_degree = 2*basisCubDegree;
705 auto elemBasisCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree, EPolyType::POLYTYPE_MAX, true);
706 basisPointsRange(dim,0) = range_type(0, elemBasisCub->getNumPoints());
707 numBasisEvalPoints += elemBasisCub->getNumPoints();
708 maxNumBasisEvalPoints = elemBasisCub->getNumPoints();
709 basisCubPoints[dim][0] = host_view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
710 basisCubWeights[dim][0] = host_view_type("basisCubWeights",elemBasisCub->getNumPoints());
711 elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
712
713 cub_degree = basisCubDegree + targetCubDegree;
714 auto elemTargetCub = cub_factory.create<HostDeviceType, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree, EPolyType::POLYTYPE_MAX, true);
715 targetPointsRange(dim,0) = range_type(0, elemTargetCub->getNumPoints());
716 numTargetEvalPoints += elemTargetCub->getNumPoints();
717 maxNumTargetEvalPoints = elemTargetCub->getNumPoints();
718 targetCubPoints[dim][0] = host_view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
719 targetCubWeights[dim][0] = host_view_type("targetCubWeights",elemTargetCub->getNumPoints());
720 elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
721
722 allBasisEPoints = view_type("allBasisPoints", numBasisEvalPoints, dim);
723 auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),basisCubPoints[dim][0]);
724 Kokkos::deep_copy(Kokkos::subview(allBasisEPoints,basisPointsRange(dim,0), Kokkos::ALL()), cellBasisEPoints);
725
726 allTargetEPoints = view_type("allTargetPoints", numTargetEvalPoints, dim);
727 auto cellTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),targetCubPoints[dim][0]);
728 Kokkos::deep_copy(Kokkos::subview(allTargetEPoints,targetPointsRange(dim,0),Kokkos::ALL()), cellTargetEPoints);
729
730 allBasisDerivEPoints = view_type("allBasisDerivPoints", numBasisDerivEvalPoints, dim);
731 allTargetDerivEPoints = view_type("allTargetDerivPoints", numTargetDerivEvalPoints, dim);
732}
733
734} // Intrepid2 namespace
735#endif
736
737
738
739
740
Intrepid2_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
Intrepid2_HCURL_HEX_In_FEM.hpp
Header file for the Intrepid2::Basis_HCURL_HEX_In_FEM class.
Intrepid2_HCURL_TET_In_FEM.hpp
Header file for the Intrepid2::Basis_HCURL_TET_In_FEM class.
Intrepid2_HGRAD_QUAD_Cn_FEM.hpp
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
Intrepid2_HGRAD_TRI_Cn_FEM.hpp
Header file for the Intrepid2::Basis_HGRAD_TRI_Cn_FEM class.
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
Intrepid2::Basis_Derived_HCURL_WEDGE
Definition Intrepid2_DerivedBasis_HCURL_WEDGE.hpp:356
Intrepid2::Basis_HCURL_HEX_In_FEM
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Definition Intrepid2_HCURL_HEX_In_FEM.hpp:149
Intrepid2::Basis_HCURL_TET_In_FEM
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
Definition Intrepid2_HCURL_TET_In_FEM.hpp:175
Intrepid2::Basis_HGRAD_QUAD_Cn_FEM
Implementation of the default H(grad)-compatible FEM basis of degree n on Quadrilateral cell Implemen...
Definition Intrepid2_HGRAD_QUAD_Cn_FEM.hpp:135
Intrepid2::Basis_HGRAD_TRI_Cn_FEM
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
Definition Intrepid2_HGRAD_TRI_Cn_FEM.hpp:156
Intrepid2::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition Intrepid2_Basis.hpp:89
Intrepid2::Basis::getDofCount
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
Definition Intrepid2_Basis.hpp:866
Intrepid2::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(refSubcellViewType refSubcellPoints, const paramPointViewType paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Intrepid2::CellTools::getReferenceVertex
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.
Definition Intrepid2_CellToolsDefNodeInfo.hpp:68
Intrepid2::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition Intrepid2_DefaultCubatureFactory.hpp:45
Intrepid2::DefaultCubatureFactory::create
static Teuchos::RCP< Cubature< DeviceType, pointValueType, weightValueType > > create(unsigned topologyKey, const std::vector< ordinal_type > &degree, const EPolyType polytype=POLYTYPE_MAX, const bool symmetric=false)
Factory method.
Intrepid2::ProjectionStruct::createHVolProjectionStruct
void createHVolProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for HVOL (local-L2) projection.
Definition Intrepid2_ProjectionStructDef.hpp:683
Intrepid2::ProjectionStruct::createHCurlProjectionStruct
void createHCurlProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetCurlCubDegre)
Initialize the ProjectionStruct for HCURL projections.
Definition Intrepid2_ProjectionStructDef.hpp:357
Intrepid2::ProjectionStruct::createHGradProjectionStruct
void createHGradProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetGradCubDegre)
Initialize the ProjectionStruct for HGRAD projections.
Definition Intrepid2_ProjectionStructDef.hpp:201
Intrepid2::ProjectionStruct::createHDivProjectionStruct
void createHDivProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree, const ordinal_type targetDivCubDegre)
Initialize the ProjectionStruct for HDIV projections.
Definition Intrepid2_ProjectionStructDef.hpp:531
Intrepid2::ProjectionStruct::createL2ProjectionStruct
void createL2ProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for L2 projections.
Definition Intrepid2_ProjectionStructDef.hpp:33
Intrepid2::RefSubcellParametrization::get
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
Definition Intrepid2_CellDataDef.hpp:69
Generated by doxygen 1.9.8