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