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 
29 namespace Intrepid2 {
30 
31 template<typename DeviceType, typename ValueType>
32 template<typename BasisPtrType>
33 void 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 
199 template<typename DeviceType, typename ValueType>
200 template<typename BasisPtrType>
201 void 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 
355 template<typename DeviceType, typename ValueType>
356 template<typename BasisPtrType>
357 void 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 
529 template<typename DeviceType, typename ValueType>
530 template<typename BasisPtrType>
531 void 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 
681 template<typename DeviceType, typename ValueType>
682 template<typename BasisPtrType>
683 void 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::DefaultCubatureFactory
A factory class that generates specific instances of cubatures.
Definition: Intrepid2_DefaultCubatureFactory.hpp:45
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::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:174
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::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
Intrepid2::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid2_Basis.hpp:38
Intrepid2::Basis_Derived_HCURL_WEDGE
Definition: Intrepid2_DerivedBasis_HCURL_WEDGE.hpp:354
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::Basis_HCURL_HEX_In_FEM
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Definition: Intrepid2_HCURL_HEX_In_FEM.hpp:148
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
Intrepid2_HGRAD_TRI_Cn_FEM.hpp
Header file for the Intrepid2::Basis_HGRAD_TRI_Cn_FEM class.
Intrepid2_HCURL_TET_In_FEM.hpp
Header file for the Intrepid2::Basis_HCURL_TET_In_FEM class.
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_HGRAD_QUAD_Cn_FEM.hpp
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
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::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::ProjectionStruct::createL2ProjectionStruct
void createL2ProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for L2 projections.
Definition: Intrepid2_ProjectionStructDef.hpp:33
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:155
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::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::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:134
Generated by   doxygen 1.8.5