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 namespace Experimental {
65 template<typename SpT, typename ValueType>
66 template<typename BasisPtrType>
67 void ProjectionStruct<SpT,ValueType>::createL2ProjectionStruct(const BasisPtrType cellBasis,
68  const ordinal_type targetCubDegree) {
69  const auto cellTopo = cellBasis->getBaseCellTopology();
70  std::string name = cellBasis->getName();
71  ordinal_type dim = cellTopo.getDimension();
72  numBasisEvalPoints = 0;
73  numBasisDerivEvalPoints = 0;
74  numTargetEvalPoints = 0;
75  numTargetDerivEvalPoints = 0;
76  const ordinal_type edgeDim = 1;
77  const ordinal_type faceDim = 2;
78 
79  ordinal_type basisCubDegree = cellBasis->getDegree();
80  ordinal_type edgeBasisCubDegree, faceBasisCubDegree;
81 
82  ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
83  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
84  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
85 
86  numBasisEvalPoints += numVertices;
87  numTargetEvalPoints += numVertices;
88  view_type dofCoords("dofCoords", cellBasis->getCardinality(), dim);
89  for(ordinal_type iv=0; iv<numVertices; ++iv) {
90  basisPointsRange[0][iv] = range_type(iv, iv+1);
91  basisCubPoints[0][iv] = view_type("basisCubPoints",1,dim);
92  targetPointsRange[0][iv] = range_type(iv, iv+1);
93  targetCubPoints[0][iv] = view_type("targetCubPoints",1,dim);
94  ordinal_type idof = cellBasis->getDofOrdinal(0, iv, 0);
95  for(ordinal_type d=0; d<dim; d++) {
96  basisCubPoints[0][iv](0,d) = dofCoords(idof,d);
97  targetCubPoints[0][iv](0,d) = dofCoords(idof,d);
98  }
99  }
100 
101  if ((name == "Intrepid2_HCURL_QUAD_In_FEM") || (name == "Intrepid2_HCURL_HEX_In_FEM") ||
102  (name == "Intrepid2_HCURL_TRI_In_FEM") || (name == "Intrepid2_HCURL_TET_In_FEM")) {
103  edgeBasisCubDegree = basisCubDegree - 1;
104  faceBasisCubDegree = basisCubDegree;
105  }
106  else if ((name == "Intrepid2_HDIV_QUAD_In_FEM") || (name == "Intrepid2_HDIV_HEX_In_FEM") ||
107  (name == "Intrepid2_HDIV_TRI_In_FEM") || (name == "Intrepid2_HDIV_TET_In_FEM")) {
108  edgeBasisCubDegree = basisCubDegree - 1;
109  faceBasisCubDegree = basisCubDegree -1;
110  }
111  else {
112  edgeBasisCubDegree = basisCubDegree;
113  faceBasisCubDegree = basisCubDegree;
114  }
115 
116  DefaultCubatureFactory cub_factory;
117  for(ordinal_type ie=0; ie<numEdges; ++ie) {
118  ordinal_type cub_degree = 2*edgeBasisCubDegree;
119  subCellTopologyKey[edgeDim][ie] = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
120  auto edgeBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
121  basisPointsRange[edgeDim][ie] = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
122  numBasisEvalPoints += edgeBasisCub->getNumPoints();
123  basisCubPoints[edgeDim][ie] = view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
124  basisCubWeights[edgeDim][ie] = view_type("basisCubWeights",edgeBasisCub->getNumPoints());
125  edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
126 
127  cub_degree = edgeBasisCubDegree + targetCubDegree;
128  auto edgeTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
129  targetPointsRange[edgeDim][ie] = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
130  numTargetEvalPoints += edgeTargetCub->getNumPoints();
131  targetCubPoints[edgeDim][ie] = view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
132  targetCubWeights[edgeDim][ie] = view_type("targetCubWeights",edgeTargetCub->getNumPoints());
133  edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
134  }
135 
136  for(ordinal_type iface=0; iface<numFaces; ++iface) {
137  ordinal_type cub_degree = 2*faceBasisCubDegree;
138  subCellTopologyKey[faceDim][iface] = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
139  auto faceBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
140  basisPointsRange[faceDim][iface] = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
141  numBasisEvalPoints += faceBasisCub->getNumPoints();
142  basisCubPoints[faceDim][iface] = view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
143  basisCubWeights[faceDim][iface] = view_type("basisCubWeights",faceBasisCub->getNumPoints());
144  faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
145 
146  cub_degree = faceBasisCubDegree + targetCubDegree;
147  auto faceTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
148  targetPointsRange[faceDim][iface] = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
149  numTargetEvalPoints += faceTargetCub->getNumPoints();
150  targetCubPoints[faceDim][iface] = view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
151  targetCubWeights[faceDim][iface] = view_type("targetCubWeights",faceTargetCub->getNumPoints());
152  faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
153  }
154  subCellTopologyKey[dim][0] = cellBasis->getBaseCellTopology().getBaseKey();
155  if(cellBasis->getDofCount(dim,0)>0) {
156  ordinal_type cub_degree = 2*basisCubDegree;
157  auto elemBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
158  basisPointsRange[dim][0] = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
159  numBasisEvalPoints += elemBasisCub->getNumPoints();
160  basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
161  basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
162  elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
163 
164  cub_degree = basisCubDegree + targetCubDegree;
165  auto elemTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
166  targetPointsRange[dim][0] = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
167  numTargetEvalPoints += elemTargetCub->getNumPoints();
168  targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
169  targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
170  elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
171  }
172 }
173 
174 template<typename SpT, typename ValueType>
175 template<typename BasisPtrType>
176 void ProjectionStruct<SpT,ValueType>::createHGradProjectionStruct(const BasisPtrType cellBasis,
177  const ordinal_type targetCubDegree,
178  const ordinal_type targetGradCubDegree) {
179  const auto cellTopo = cellBasis->getBaseCellTopology();
180  std::string name = cellBasis->getName();
181  ordinal_type dim = cellTopo.getDimension();
182  numBasisEvalPoints = 0;
183  numBasisDerivEvalPoints = 0;
184  numTargetEvalPoints = 0;
185  numTargetDerivEvalPoints = 0;
186  const ordinal_type edgeDim = 1;
187  const ordinal_type faceDim = 2;
188 
189  ordinal_type basisCubDegree = cellBasis->getDegree();
190  ordinal_type edgeBasisCubDegree = basisCubDegree;
191  ordinal_type faceBasisCubDegree = basisCubDegree;
192 
193  ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
194  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount(): 0;
195  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
196 
197  numBasisEvalPoints += numVertices;
198  numTargetEvalPoints += numVertices;
199  view_type dofCoords("dofCoords", cellBasis->getCardinality(), dim);
200  for(ordinal_type iv=0; iv<numVertices; ++iv) {
201  basisPointsRange[0][iv] = range_type(iv, iv+1);
202  basisCubPoints[0][iv] = view_type("basisCubPoints",1,dim);
203  targetPointsRange[0][iv] = range_type(iv, iv+1);
204  targetCubPoints[0][iv] = view_type("targetCubPoints",1,dim);
205  ordinal_type idof = cellBasis->getDofOrdinal(0, iv, 0);
206  for(ordinal_type d=0; d<dim; d++) {
207  basisCubPoints[0][iv](0,d) = dofCoords(idof,d);
208  targetCubPoints[0][iv](0,d) = dofCoords(idof,d);
209  }
210  }
211 
212  DefaultCubatureFactory cub_factory;
213  for(ordinal_type ie=0; ie<numEdges; ++ie) {
214  ordinal_type cub_degree = 2*edgeBasisCubDegree;
215  subCellTopologyKey[edgeDim][ie] = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
216  auto edgeBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
217  basisDerivPointsRange[edgeDim][ie] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+edgeBasisCub->getNumPoints());
218  numBasisDerivEvalPoints += edgeBasisCub->getNumPoints();
219  basisDerivCubPoints[edgeDim][ie] = view_type("basisDerivCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
220  basisDerivCubWeights[edgeDim][ie] = view_type("basisDerivCubWeights",edgeBasisCub->getNumPoints());
221  edgeBasisCub->getCubature(basisDerivCubPoints[edgeDim][ie], basisDerivCubWeights[edgeDim][ie]);
222 
223  cub_degree = edgeBasisCubDegree + targetGradCubDegree;
224  auto edgeTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellBasis->getBaseCellTopology().getKey(edgeDim, ie), cub_degree);
225  targetDerivPointsRange[edgeDim][ie] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+edgeTargetCub->getNumPoints());
226  numTargetDerivEvalPoints += edgeTargetCub->getNumPoints();
227  targetDerivCubPoints[edgeDim][ie] = view_type("targetDerivCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
228  targetDerivCubWeights[edgeDim][ie] = view_type("targetDerivCubWeights",edgeTargetCub->getNumPoints());
229  edgeTargetCub->getCubature(targetDerivCubPoints[edgeDim][ie], targetDerivCubWeights[edgeDim][ie]);
230  }
231 
232  for(ordinal_type iface=0; iface<numFaces; ++iface) {
233  ordinal_type cub_degree = 2*faceBasisCubDegree;
234  subCellTopologyKey[faceDim][iface] = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
235  auto faceBasisGradCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
236  basisDerivPointsRange[faceDim][iface] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisGradCub->getNumPoints());
237  numBasisDerivEvalPoints += faceBasisGradCub->getNumPoints();
238  basisDerivCubPoints[faceDim][iface] = view_type("basisDerivCubPoints",faceBasisGradCub->getNumPoints(),faceDim);
239  basisDerivCubWeights[faceDim][iface] = view_type("basisDerivCubWeights",faceBasisGradCub->getNumPoints());
240  faceBasisGradCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
241 
242  cub_degree = faceBasisCubDegree + targetGradCubDegree;
243  auto faceTargetDerivCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
244  targetDerivPointsRange[faceDim][iface] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
245  numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
246  targetDerivCubPoints[faceDim][iface] = view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
247  targetDerivCubWeights[faceDim][iface] = view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
248  faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
249  }
250  subCellTopologyKey[dim][0] = cellBasis->getBaseCellTopology().getBaseKey();
251  if(cellBasis->getDofCount(dim,0)>0) {
252  ordinal_type cub_degree = 2*basisCubDegree;
253  auto elemBasisGradCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
254  basisDerivPointsRange[dim][0] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisGradCub->getNumPoints());
255  numBasisDerivEvalPoints += elemBasisGradCub->getNumPoints();
256  basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisGradCub->getNumPoints(),dim);
257  basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisGradCub->getNumPoints());
258  elemBasisGradCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
259 
260  cub_degree = basisCubDegree + targetGradCubDegree;
261  auto elemTargetGradCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
262  targetDerivPointsRange[dim][0] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetGradCub->getNumPoints());
263  numTargetDerivEvalPoints += elemTargetGradCub->getNumPoints();
264  targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetGradCub->getNumPoints(),dim);
265  targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetGradCub->getNumPoints());
266  elemTargetGradCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
267  }
268 }
269 
270 template<typename SpT, typename ValueType>
271 template<typename BasisPtrType>
272 void ProjectionStruct<SpT,ValueType>::createHCurlProjectionStruct(const BasisPtrType cellBasis,
273  const ordinal_type targetCubDegree,
274  const ordinal_type targetCurlCubDegre) {
275  const auto cellTopo = cellBasis->getBaseCellTopology();
276  std::string name = cellBasis->getName();
277  ordinal_type dim = cellTopo.getDimension();
278  numBasisEvalPoints = 0;
279  numBasisDerivEvalPoints = 0;
280  numTargetEvalPoints = 0;
281  numTargetDerivEvalPoints = 0;
282  const ordinal_type edgeDim = 1;
283  const ordinal_type faceDim = 2;
284 
285  ordinal_type basisCubDegree = cellBasis->getDegree();
286  ordinal_type edgeBasisCubDegree = basisCubDegree - 1;
287  ordinal_type faceBasisCubDegree = basisCubDegree;
288  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
289  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
290 
291  DefaultCubatureFactory cub_factory;
292  for(ordinal_type ie=0; ie<numEdges; ++ie) {
293  ordinal_type cub_degree = 2*edgeBasisCubDegree;
294  subCellTopologyKey[edgeDim][ie] = cellBasis->getBaseCellTopology().getKey(edgeDim, ie);
295  auto edgeBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[edgeDim][ie], cub_degree);
296  basisPointsRange[edgeDim][ie] = range_type(numBasisEvalPoints, numBasisEvalPoints+edgeBasisCub->getNumPoints());
297  numBasisEvalPoints += edgeBasisCub->getNumPoints();
298  basisCubPoints[edgeDim][ie] = view_type("basisCubPoints",edgeBasisCub->getNumPoints(),edgeDim);
299  basisCubWeights[edgeDim][ie] = view_type("basisCubWeights",edgeBasisCub->getNumPoints());
300  edgeBasisCub->getCubature(basisCubPoints[edgeDim][ie], basisCubWeights[edgeDim][ie]);
301 
302  cub_degree = edgeBasisCubDegree + targetCubDegree;
303  auto edgeTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[edgeDim][ie], cub_degree);
304  targetPointsRange[edgeDim][ie] = range_type(numTargetEvalPoints, numTargetEvalPoints+edgeTargetCub->getNumPoints());
305  numTargetEvalPoints += edgeTargetCub->getNumPoints();
306  targetCubPoints[edgeDim][ie] = view_type("targetCubPoints",edgeTargetCub->getNumPoints(),edgeDim);
307  targetCubWeights[edgeDim][ie] = view_type("targetCubWeights",edgeTargetCub->getNumPoints());
308  edgeTargetCub->getCubature(targetCubPoints[edgeDim][ie], targetCubWeights[edgeDim][ie]);
309  }
310 
311  for(ordinal_type iface=0; iface<numFaces; ++iface) {
312  ordinal_type cub_degree = 2*faceBasisCubDegree;
313  subCellTopologyKey[faceDim][iface] = cellBasis->getBaseCellTopology().getKey(faceDim, iface);
314  auto faceBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
315  basisPointsRange[faceDim][iface] = range_type(numBasisEvalPoints, numBasisEvalPoints+faceBasisCub->getNumPoints());
316  numBasisEvalPoints += faceBasisCub->getNumPoints();
317  basisCubPoints[faceDim][iface] = view_type("basisCubPoints",faceBasisCub->getNumPoints(),faceDim);
318  basisCubWeights[faceDim][iface] = view_type("basisCubWeights",faceBasisCub->getNumPoints());
319  faceBasisCub->getCubature(basisCubPoints[faceDim][iface], basisCubWeights[faceDim][iface]);
320 
321  basisDerivPointsRange[faceDim][iface] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+faceBasisCub->getNumPoints());
322  numBasisDerivEvalPoints += faceBasisCub->getNumPoints();
323  basisDerivCubPoints[faceDim][iface] = view_type("basisDerivCubPoints",faceBasisCub->getNumPoints(),faceDim);
324  basisDerivCubWeights[faceDim][iface] = view_type("basisDerivCubWeights",faceBasisCub->getNumPoints());
325  faceBasisCub->getCubature(basisDerivCubPoints[faceDim][iface], basisDerivCubWeights[faceDim][iface]);
326 
327  cub_degree = faceBasisCubDegree + targetCubDegree;
328  auto faceTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
329  targetPointsRange[faceDim][iface] = range_type(numTargetEvalPoints, numTargetEvalPoints+faceTargetCub->getNumPoints());
330  numTargetEvalPoints += faceTargetCub->getNumPoints();
331  targetCubPoints[faceDim][iface] = view_type("targetCubPoints",faceTargetCub->getNumPoints(),faceDim);
332  targetCubWeights[faceDim][iface] = view_type("targetCubWeights",faceTargetCub->getNumPoints());
333  faceTargetCub->getCubature(targetCubPoints[faceDim][iface], targetCubWeights[faceDim][iface]);
334 
335  cub_degree = faceBasisCubDegree + targetCurlCubDegre;
336  auto faceTargetDerivCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[faceDim][iface], cub_degree);
337  targetDerivPointsRange[faceDim][iface] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+faceTargetDerivCub->getNumPoints());
338  numTargetDerivEvalPoints += faceTargetDerivCub->getNumPoints();
339  targetDerivCubPoints[faceDim][iface] = view_type("targetDerivCubPoints",faceTargetDerivCub->getNumPoints(),faceDim);
340  targetDerivCubWeights[faceDim][iface] = view_type("targetDerivCubWeights",faceTargetDerivCub->getNumPoints());
341  faceTargetDerivCub->getCubature(targetDerivCubPoints[faceDim][iface], targetDerivCubWeights[faceDim][iface]);
342  }
343 
344  subCellTopologyKey[dim][0] = cellBasis->getBaseCellTopology().getBaseKey();
345  if(cellBasis->getDofCount(dim,0)>0) {
346  ordinal_type cub_degree = 2*basisCubDegree;
347  auto elemBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
348  basisPointsRange[dim][0] = range_type(numBasisEvalPoints, numBasisEvalPoints+elemBasisCub->getNumPoints());
349  numBasisEvalPoints += elemBasisCub->getNumPoints();
350  basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
351  basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
352  elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
353 
354  basisDerivPointsRange[dim][0] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisCub->getNumPoints());
355  numBasisDerivEvalPoints += elemBasisCub->getNumPoints();
356  basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisCub->getNumPoints(),dim);
357  basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisCub->getNumPoints());
358  elemBasisCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
359 
360  cub_degree = basisCubDegree + targetCubDegree;
361  auto elemTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
362  targetPointsRange[dim][0] = range_type(numTargetEvalPoints, numTargetEvalPoints+elemTargetCub->getNumPoints());
363  numTargetEvalPoints += elemTargetCub->getNumPoints();
364  targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
365  targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
366  elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
367 
368  cub_degree = basisCubDegree + targetCurlCubDegre;
369  auto elemTargetCurlCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
370  targetDerivPointsRange[dim][0] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetCurlCub->getNumPoints());
371  numTargetDerivEvalPoints += elemTargetCurlCub->getNumPoints();
372  targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetCurlCub->getNumPoints(),dim);
373  targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetCurlCub->getNumPoints());
374  elemTargetCurlCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
375  }
376 
377 }
378 
379 template<typename SpT, typename ValueType>
380 template<typename BasisPtrType>
381 void ProjectionStruct<SpT,ValueType>::createHDivProjectionStruct(const BasisPtrType cellBasis,
382  const ordinal_type targetCubDegree,
383  const ordinal_type targetDivCubDegre) {
384  const auto cellTopo = cellBasis->getBaseCellTopology();
385  std::string name = cellBasis->getName();
386  ordinal_type dim = cellTopo.getDimension();
387  numBasisEvalPoints = 0;
388  numBasisDerivEvalPoints = 0;
389  numTargetEvalPoints = 0;
390  numTargetDerivEvalPoints = 0;
391  const ordinal_type sideDim = dim - 1;
392 
393  ordinal_type basisCubDegree = cellBasis->getDegree();
394  ordinal_type sideBasisCubDegree = basisCubDegree - 1;
395  ordinal_type numSides = cellTopo.getSideCount()*ordinal_type(cellBasis->getDofCount(sideDim, 0) > 0);
396 
397  Basis<host_space_type,ValueType,ValueType> *hcurlBasis = NULL;
398  if(name.find("HEX")!=std::string::npos)
399  hcurlBasis = new Basis_HCURL_HEX_In_FEM<host_space_type,ValueType,ValueType>(cellBasis->getDegree());
400  else if(name.find("TET")!=std::string::npos)
401  hcurlBasis = new Basis_HCURL_TET_In_FEM<host_space_type,ValueType,ValueType>(cellBasis->getDegree());
402  else if(name.find("QUAD")!=std::string::npos)
403  hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM<host_space_type,ValueType,ValueType>(cellBasis->getDegree());
404  else if(name.find("TRI")!=std::string::npos)
405  hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM<host_space_type,ValueType,ValueType>(cellBasis->getDegree());
406  else {
407  std::stringstream ss;
408  ss << ">>> ERROR (Intrepid2::ProjectionTools::createHDivProjectionStruct): "
409  << "Method not implemented for basis " << name;
410  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
411  }
412 
413  bool haveHCurlConstraint = (hcurlBasis != NULL) && (hcurlBasis->getDofCount(dim,0) > 0);
414  if(hcurlBasis != NULL) delete hcurlBasis;
415 
416  DefaultCubatureFactory cub_factory;
417 
418  for(ordinal_type is=0; is<numSides; ++is) {
419  ordinal_type cub_degree = 2*sideBasisCubDegree;
420  subCellTopologyKey[sideDim][is] = cellBasis->getBaseCellTopology().getKey(sideDim, is);
421  auto sideBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[sideDim][is], cub_degree);
422  basisPointsRange[sideDim][is] = range_type(numBasisEvalPoints, numBasisEvalPoints+sideBasisCub->getNumPoints());
423  numBasisEvalPoints += sideBasisCub->getNumPoints();
424  basisCubPoints[sideDim][is] = view_type("basisCubPoints",sideBasisCub->getNumPoints(),sideDim);
425  basisCubWeights[sideDim][is] = view_type("basisCubWeights",sideBasisCub->getNumPoints());
426  sideBasisCub->getCubature(basisCubPoints[sideDim][is], basisCubWeights[sideDim][is]);
427 
428  cub_degree = sideBasisCubDegree + targetCubDegree;
429  auto sideTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[sideDim][is], cub_degree);
430  targetPointsRange[sideDim][is] = range_type(numTargetEvalPoints, numTargetEvalPoints+sideTargetCub->getNumPoints());
431  numTargetEvalPoints += sideTargetCub->getNumPoints();
432  targetCubPoints[sideDim][is] = view_type("targetCubPoints",sideTargetCub->getNumPoints(),sideDim);
433  targetCubWeights[sideDim][is] = view_type("targetCubWeights",sideTargetCub->getNumPoints());
434  sideTargetCub->getCubature(targetCubPoints[sideDim][is], targetCubWeights[sideDim][is]);
435  }
436 
437  subCellTopologyKey[dim][0] = cellBasis->getBaseCellTopology().getBaseKey();
438  if(cellBasis->getDofCount(dim,0)>0) {
439  ordinal_type cub_degree = 2*basisCubDegree - 1;
440  auto elemBasisDivCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
441  basisDerivPointsRange[dim][0] = range_type(numBasisDerivEvalPoints, numBasisDerivEvalPoints+elemBasisDivCub->getNumPoints());
442  numBasisDerivEvalPoints += elemBasisDivCub->getNumPoints();
443  basisDerivCubPoints[dim][0] = view_type("basisDerivCubPoints",elemBasisDivCub->getNumPoints(),dim);
444  basisDerivCubWeights[dim][0] = view_type("basisDerivCubWeights",elemBasisDivCub->getNumPoints());
445  elemBasisDivCub->getCubature(basisDerivCubPoints[dim][0], basisDerivCubWeights[dim][0]);
446 
447  cub_degree = basisCubDegree - 1 + targetDivCubDegre;
448  auto elemTargetDivCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
449  targetDerivPointsRange[dim][0] = range_type(numTargetDerivEvalPoints, numTargetDerivEvalPoints+elemTargetDivCub->getNumPoints());
450  numTargetDerivEvalPoints += elemTargetDivCub->getNumPoints();
451  targetDerivCubPoints[dim][0] = view_type("targetDerivCubPoints",elemTargetDivCub->getNumPoints(),dim);
452  targetDerivCubWeights[dim][0] = view_type("targetDerivCubWeights",elemTargetDivCub->getNumPoints());
453  elemTargetDivCub->getCubature(targetDerivCubPoints[dim][0], targetDerivCubWeights[dim][0]);
454 
455  if(haveHCurlConstraint)
456  {
457  cub_degree = 2*basisCubDegree;
458  auto elemBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
459  basisPointsRange[dim][0] = range_type(numBasisEvalPoints, numBasisEvalPoints + elemBasisCub->getNumPoints());
460  numBasisEvalPoints += elemBasisCub->getNumPoints();
461  basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
462  basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
463  elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
464 
465  cub_degree = basisCubDegree + targetCubDegree;
466  auto elemTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(subCellTopologyKey[dim][0], cub_degree);
467  targetPointsRange[dim][0] = range_type(numTargetEvalPoints, numTargetEvalPoints + elemTargetCub->getNumPoints());
468  numTargetEvalPoints += elemTargetCub->getNumPoints();
469  targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
470  targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
471  elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
472  }
473  }
474 }
475 
476 template<typename SpT, typename ValueType>
477 template<typename BasisPtrType>
478 void ProjectionStruct<SpT,ValueType>::createHVolProjectionStruct(const BasisPtrType cellBasis,
479  const ordinal_type targetCubDegree) {
480  const auto cellTopo = cellBasis->getBaseCellTopology();
481  ordinal_type dim = cellTopo.getDimension();
482  numBasisEvalPoints = 0;
483  numBasisDerivEvalPoints = 0;
484  numTargetEvalPoints = 0;
485  numTargetDerivEvalPoints = 0;
486 
487  ordinal_type basisCubDegree = cellBasis->getDegree();
488  DefaultCubatureFactory cub_factory;
489  subCellTopologyKey[dim][0] = cellBasis->getBaseCellTopology().getBaseKey();
490  if(cellBasis->getDofCount(dim,0)>0) {
491 
492  ordinal_type cub_degree = 2*basisCubDegree;
493  auto elemBasisCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree);
494  basisPointsRange[dim][0] = range_type(0, elemBasisCub->getNumPoints());
495  numBasisEvalPoints += elemBasisCub->getNumPoints();
496  basisCubPoints[dim][0] = view_type("basisCubPoints",elemBasisCub->getNumPoints(),dim);
497  basisCubWeights[dim][0] = view_type("basisCubWeights",elemBasisCub->getNumPoints());
498  elemBasisCub->getCubature(basisCubPoints[dim][0], basisCubWeights[dim][0]);
499 
500  cub_degree = basisCubDegree + targetCubDegree;
501  auto elemTargetCub = cub_factory.create<host_space_type, ValueType, ValueType>(cellTopo.getBaseKey(), cub_degree);
502  targetPointsRange[dim][0] = range_type(0, elemTargetCub->getNumPoints());
503  numTargetEvalPoints += elemTargetCub->getNumPoints();
504  targetCubPoints[dim][0] = view_type("targetCubPoints",elemTargetCub->getNumPoints(),dim);
505  targetCubWeights[dim][0] = view_type("targetCubWeights",elemTargetCub->getNumPoints());
506  elemTargetCub->getCubature(targetCubPoints[dim][0], targetCubWeights[dim][0]);
507  }
508 }
509 
510 }
511 }
512 #endif
513 
514 
515 
516 
517 
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:76
Intrepid2::Experimental::ProjectionStruct::createL2ProjectionStruct
void createL2ProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for L2 projections.
Definition: Intrepid2_ProjectionStructDef.hpp:67
Intrepid2::Experimental::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:381
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:204
Intrepid2::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid2_Basis.hpp:93
Intrepid2::Basis_HCURL_HEX_In_FEM
Implementation of the default H(curl)-compatible FEM basis on Hexahedral cell.
Definition: Intrepid2_HCURL_HEX_In_FEM.hpp:181
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
Intrepid2::DefaultCubatureFactory::create
static Teuchos::RCP< Cubature< ExecSpaceType, pointValueType, weightValueType > > create(unsigned topologyKey, const std::vector< ordinal_type > &degree, const EPolyType polytype=POLYTYPE_MAX)
Factory method.
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_HGRAD_QUAD_Cn_FEM.hpp
Header file for the Intrepid2::Basis_HGRAD_QUAD_Cn_FEM class.
Intrepid2::Experimental::ProjectionStruct::createHVolProjectionStruct
void createHVolProjectionStruct(const BasisPtrType cellBasis, const ordinal_type targetCubDegree)
Initialize the ProjectionStruct for HVOL (local-L2) projection.
Definition: Intrepid2_ProjectionStructDef.hpp:478
Intrepid2::Experimental::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:272
Intrepid2::Experimental::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:176
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:181
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:166
Generated by   doxygen 1.8.5