Intrepid2
Intrepid2_ProjectionToolsDefL2.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 
49 #ifndef __INTREPID2_PROJECTIONTOOLSDEFL2_HPP__
50 #define __INTREPID2_PROJECTIONTOOLSDEFL2_HPP__
51 
53 #include "Intrepid2_ArrayTools.hpp"
55 
56 
57 namespace Intrepid2 {
58 namespace Experimental {
59 
60 
61 template<typename ViewType1, typename ViewType2, typename ViewType3,
62 typename ViewType4, typename ViewType5>
64  ViewType1 basisCoeffs_;
65  const ViewType2 tagToOrdinal_;
66  const ViewType3 targetEPointsRange_;
67  const ViewType4 targetAtTargetEPoints_;
68  const ViewType5 basisAtTargetEPoints_;
69  ordinal_type numVertices_;
70 
71 
72  ComputeBasisCoeffsOnVertices_L2(ViewType1 basisCoeffs, ViewType2 tagToOrdinal, ViewType3 targetEPointsRange,
73  ViewType4 targetAtTargetEPoints, ViewType5 basisAtTargetEPoints, ordinal_type numVertices) :
74  basisCoeffs_(basisCoeffs), tagToOrdinal_(tagToOrdinal), targetEPointsRange_(targetEPointsRange),
75  targetAtTargetEPoints_(targetAtTargetEPoints), basisAtTargetEPoints_(basisAtTargetEPoints), numVertices_(numVertices) {}
76 
77  void
78  KOKKOS_INLINE_FUNCTION
79  operator()(const ordinal_type ic) const {
80  for(ordinal_type iv=0; iv<numVertices_; ++iv) {
81  ordinal_type idof = tagToOrdinal_(0, iv, 0);
82  ordinal_type pt = targetEPointsRange_(0,iv).first;
83  //the value of the basis at the vertex might be different than 1; HGrad basis, so the function is scalar
84  basisCoeffs_(ic,idof) = targetAtTargetEPoints_(ic,pt)/basisAtTargetEPoints_(ic,idof,pt,0);
85  }
86  }
87 };
88 
89 
90 template<typename ViewType1, typename ViewType2, typename ViewType3,
91 typename ViewType4, typename ViewType5, typename ViewType6>
93  const ViewType1 basisCoeffs_;
94  const ViewType2 negPartialProj_;
95  const ViewType2 basisDofDofAtBasisEPoints_;
96  const ViewType2 basisAtBasisEPoints_;
97  const ViewType3 basisEWeights_;
98  const ViewType2 wBasisDofAtBasisEPoints_;
99  const ViewType3 targetEWeights_;
100  const ViewType2 basisAtTargetEPoints_;
101  const ViewType2 wBasisDofAtTargetEPoints_;
102  const ViewType4 computedDofs_;
103  const ViewType5 tagToOrdinal_;
104  const ViewType6 targetAtTargetEPoints_;
105  const ViewType2 targetTanAtTargetEPoints_;
106  const ViewType2 refEdgesVec_;
107  ordinal_type fieldDim_;
108  ordinal_type edgeCardinality_;
109  ordinal_type offsetBasis_;
110  ordinal_type offsetTarget_;
111  ordinal_type numVertexDofs_;
112  ordinal_type edgeDim_;
113  ordinal_type iedge_;
114 
115  ComputeBasisCoeffsOnEdges_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 basisDofDofAtBasisEPoints,
116  const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisDofAtBasisEPoints, const ViewType3 targetEWeights,
117  const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal,
118  const ViewType6 targetAtTargetEPoints, const ViewType2 targetTanAtTargetEPoints, const ViewType2 refEdgesVec,
119  ordinal_type fieldDim, ordinal_type edgeCardinality, ordinal_type offsetBasis,
120  ordinal_type offsetTarget, ordinal_type numVertexDofs, ordinal_type edgeDim, ordinal_type iedge) :
121  basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), basisDofDofAtBasisEPoints_(basisDofDofAtBasisEPoints),
122  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights),
123  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
124  computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), targetAtTargetEPoints_(targetAtTargetEPoints),
125  targetTanAtTargetEPoints_(targetTanAtTargetEPoints), refEdgesVec_(refEdgesVec),
126  fieldDim_(fieldDim), edgeCardinality_(edgeCardinality), offsetBasis_(offsetBasis),
127  offsetTarget_(offsetTarget), numVertexDofs_(numVertexDofs), edgeDim_(edgeDim), iedge_(iedge)
128  {}
129 
130  void
131  KOKKOS_INLINE_FUNCTION
132  operator()(const ordinal_type ic) const {
133  for(ordinal_type j=0; j <edgeCardinality_; ++j) {
134  ordinal_type jdof =tagToOrdinal_(edgeDim_, iedge_, j);
135  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
136  for(ordinal_type d=0; d <fieldDim_; ++d)
137  basisDofDofAtBasisEPoints_(ic,j,iq) += basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d)*refEdgesVec_(iedge_,d);
138  wBasisDofAtBasisEPoints_(ic,j,iq) = basisDofDofAtBasisEPoints_(ic,j,iq)*basisEWeights_(iq);
139  }
140  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
141  for(ordinal_type d=0; d <fieldDim_; ++d)
142  wBasisDofAtTargetEPoints_(ic,j,iq) += basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d)*refEdgesVec_(iedge_,d)*targetEWeights_(iq);
143  }
144  }
145 
146  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq)
147  for(ordinal_type d=0; d <fieldDim_; ++d)
148  targetTanAtTargetEPoints_(ic,iq) += targetAtTargetEPoints_(ic,offsetTarget_+iq,d)*refEdgesVec_(iedge_,d);
149 
150  for(ordinal_type j=0; j <numVertexDofs_; ++j) {
151  ordinal_type jdof = computedDofs_(j);
152  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
153  for(ordinal_type d=0; d <fieldDim_; ++d)
154  negPartialProj_(ic,iq) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d)*refEdgesVec_(iedge_,d);
155  }
156  }
157 };
158 
159 template<typename ViewType1, typename ViewType2, typename ViewType3,
160 typename ViewType4, typename ViewType5, typename ViewType6, typename ViewType7, typename ViewType8>
162  const ViewType1 basisCoeffs_;
163  const ViewType2 negPartialProj_;
164  const ViewType2 faceBasisDofAtBasisEPoints_;
165  const ViewType2 basisAtBasisEPoints_;
166  const ViewType3 basisEWeights_;
167  const ViewType2 wBasisDofAtBasisEPoints_;
168  const ViewType3 targetEWeights_;
169  const ViewType2 basisAtTargetEPoints_;
170  const ViewType2 wBasisDofAtTargetEPoints_;
171  const ViewType4 computedDofs_;
172  const ViewType5 tagToOrdinal_;
173  const ViewType6 orts_;
174  const ViewType7 targetAtTargetEPoints_;
175  const ViewType2 targetDofAtTargetEPoints_;
176  const ViewType2 faceCoeff_;
177  const ViewType8 faceParametrization_;
178  ordinal_type fieldDim_;
179  ordinal_type faceCardinality_;
180  ordinal_type offsetBasis_;
181  ordinal_type offsetTarget_;
182  ordinal_type numVertexEdgeDofs_;
183  ordinal_type numFaces_;
184  ordinal_type faceDim_;
185  ordinal_type faceDofDim_;
186  ordinal_type dim_;
187  ordinal_type iface_;
188  unsigned topoKey_;
189  bool isHCurlBasis_, isHDivBasis_;
190 
191  ComputeBasisCoeffsOnFaces_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 faceBasisDofAtBasisEPoints,
192  const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisDofAtBasisEPoints, const ViewType3 targetEWeights,
193  const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal,
194  const ViewType6 orts, const ViewType7 targetAtTargetEPoints, const ViewType2 targetDofAtTargetEPoints, const ViewType2 faceCoeff,
195  const ViewType8 faceParametrization, ordinal_type fieldDim, ordinal_type faceCardinality, ordinal_type offsetBasis,
196  ordinal_type offsetTarget, ordinal_type numVertexEdgeDofs, ordinal_type numFaces, ordinal_type faceDim, ordinal_type faceDofDim,
197  ordinal_type dim, ordinal_type iface, unsigned topoKey, bool isHCurlBasis, bool isHDivBasis) :
198  basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), faceBasisDofAtBasisEPoints_(faceBasisDofAtBasisEPoints),
199  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights),
200  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
201  computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), orts_(orts), targetAtTargetEPoints_(targetAtTargetEPoints),
202  targetDofAtTargetEPoints_(targetDofAtTargetEPoints), faceCoeff_(faceCoeff),
203  faceParametrization_(faceParametrization),
204  fieldDim_(fieldDim), faceCardinality_(faceCardinality), offsetBasis_(offsetBasis),
205  offsetTarget_(offsetTarget), numVertexEdgeDofs_(numVertexEdgeDofs), numFaces_(numFaces),
206  faceDim_(faceDim), faceDofDim_(faceDofDim), dim_(dim), iface_(iface), topoKey_(topoKey),
207  isHCurlBasis_(isHCurlBasis), isHDivBasis_(isHDivBasis)
208  {}
209 
210  void
211  KOKKOS_INLINE_FUNCTION
212  operator()(const ordinal_type ic) const {
213 
214  ordinal_type fOrt[6];
215  orts_(ic).getFaceOrientation(fOrt, numFaces_);
216  ordinal_type ort = fOrt[iface_];
217  //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection
218 
219  typename ViewType3::value_type data[3*3];
220  auto tangentsAndNormal = ViewType3(data, dim_, dim_);
221 
222  if(isHCurlBasis_ || isHDivBasis_)
223  Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,topoKey_, iface_, ort);
224 
225  if(isHCurlBasis_) {
226  for(ordinal_type d=0; d <dim_; ++d)
227  for(ordinal_type itan=0; itan <faceDim_; ++itan) {
228  faceCoeff_(ic,d,itan) = tangentsAndNormal(itan,d);
229  }
230  } else if (isHDivBasis_) {
231  for(ordinal_type d=0; d <dim_; ++d)
232  faceCoeff_(ic,d,0) = tangentsAndNormal(dim_-1,d);
233  } else
234  faceCoeff_(ic,0,0) = 1;
235  for(ordinal_type j=0; j <faceCardinality_; ++j) {
236  ordinal_type jdof = tagToOrdinal_(faceDim_, iface_, j);
237  for(ordinal_type itan=0; itan <faceDofDim_; ++itan) {
238  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
239  for(ordinal_type d=0; d <fieldDim_; ++d)
240  faceBasisDofAtBasisEPoints_(ic,j,iq,itan) += faceCoeff_(ic,d, itan)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
241  wBasisDofAtBasisEPoints_(ic,j,iq,itan) = faceBasisDofAtBasisEPoints_(ic,j,iq,itan) * basisEWeights_(iq);
242  }
243  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
244  typename ViewType2::value_type sum=0;
245  for(ordinal_type d=0; d <fieldDim_; ++d)
246  sum += faceCoeff_(ic, d, itan)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
247  wBasisDofAtTargetEPoints_(ic,j,iq,itan) = sum * targetEWeights_(iq);
248  }
249  }
250  }
251 
252  for(ordinal_type d=0; d <fieldDim_; ++d)
253  for(ordinal_type itan=0; itan <faceDofDim_; ++itan) {
254  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq)
255  targetDofAtTargetEPoints_(ic,iq,itan) += faceCoeff_(ic, d, itan)*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
256  }
257 
258  for(ordinal_type j=0; j <numVertexEdgeDofs_; ++j) {
259  ordinal_type jdof = computedDofs_(j);
260  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
261  for(ordinal_type d=0; d <fieldDim_; ++d)
262  for(ordinal_type itan=0; itan <faceDofDim_; ++itan)
263  negPartialProj_(ic,iq,itan) -= basisCoeffs_(ic,jdof)*faceCoeff_(ic, d, itan)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
264  }
265  }
266 };
267 
268 
269 template<typename ViewType1, typename ViewType2, typename ViewType3,
270 typename ViewType4, typename ViewType5>
272  const ViewType1 basisCoeffs_;
273  const ViewType2 negPartialProj_;
274  const ViewType2 internalBasisAtBasisEPoints_;
275  const ViewType2 basisAtBasisEPoints_;
276  const ViewType3 basisEWeights_;
277  const ViewType2 wBasisAtBasisEPoints_;
278  const ViewType3 targetEWeights_;
279  const ViewType2 basisAtTargetEPoints_;
280  const ViewType2 wBasisDofAtTargetEPoints_;
281  const ViewType4 computedDofs_;
282  const ViewType5 elemDof_;
283  ordinal_type fieldDim_;
284  ordinal_type numElemDofs_;
285  ordinal_type offsetBasis_;
286  ordinal_type offsetTarget_;
287  ordinal_type numVertexEdgeFaceDofs_;
288 
289  ComputeBasisCoeffsOnCells_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 internalBasisAtBasisEPoints,
290  const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisAtBasisEPoints, const ViewType3 targetEWeights,
291  const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 elemDof,
292  ordinal_type fieldDim, ordinal_type numElemDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexEdgeFaceDofs) :
293  basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), internalBasisAtBasisEPoints_(internalBasisAtBasisEPoints),
294  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
295  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
296  computedDofs_(computedDofs), elemDof_(elemDof), fieldDim_(fieldDim), numElemDofs_(numElemDofs), offsetBasis_(offsetBasis),
297  offsetTarget_(offsetTarget), numVertexEdgeFaceDofs_(numVertexEdgeFaceDofs) {}
298 
299  void
300  KOKKOS_INLINE_FUNCTION
301  operator()(const ordinal_type ic) const {
302 
303  for(ordinal_type j=0; j <numElemDofs_; ++j) {
304  ordinal_type idof = elemDof_(j);
305  for(ordinal_type d=0; d <fieldDim_; ++d) {
306  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
307  internalBasisAtBasisEPoints_(ic,j,iq,d) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
308  wBasisAtBasisEPoints_(ic,j,iq,d) = internalBasisAtBasisEPoints_(ic,j,iq,d) * basisEWeights_(iq);
309  }
310  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
311  wBasisDofAtTargetEPoints_(ic,j,iq,d) = basisAtTargetEPoints_(ic,idof,offsetTarget_+iq,d)* targetEWeights_(iq);
312  }
313  }
314  }
315  for(ordinal_type j=0; j < numVertexEdgeFaceDofs_; ++j) {
316  ordinal_type jdof = computedDofs_(j);
317  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
318  for(ordinal_type d=0; d <fieldDim_; ++d) {
319  negPartialProj_(ic,iq,d) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
320  }
321  }
322  }
323 };
324 
325 
326 template<typename DeviceType>
327 template<typename BasisType,
328 typename ortValueType, class ...ortProperties>
329 void
330 ProjectionTools<DeviceType>::getL2EvaluationPoints(typename BasisType::ScalarViewType ePoints,
331  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
332  const BasisType* cellBasis,
334  const EvalPointsType ePointType) {
335  typedef typename BasisType::scalarType scalarType;
336  typedef Kokkos::DynRankView<scalarType,ortProperties...> ScalarViewType;
337  const auto cellTopo = cellBasis->getBaseCellTopology();
338  //const auto cellTopoKey = cellBasis->getBaseCellTopology().getKey();
339  ordinal_type dim = cellTopo.getDimension();
340  ordinal_type numCells = ePoints.extent(0);
341  const ordinal_type edgeDim = 1;
342  const ordinal_type faceDim = 2;
343 
344  ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
345  ordinal_type numEdges = (cellBasis->getDofCount(edgeDim, 0) > 0) ? cellTopo.getEdgeCount() : 0;
346  ordinal_type numFaces = (cellBasis->getDofCount(faceDim, 0) > 0) ? cellTopo.getFaceCount() : 0;
347  ordinal_type numVols = (cellBasis->getDofCount(dim, 0) > 0);
348 
349  auto ePointsRange = projStruct->getPointsRange(ePointType);
350 
351  typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamEdge, subcellParamFace;
352  if(numEdges>0)
353  subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
354  if(numFaces>0)
355  subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
356 
357  auto refTopologyKey = projStruct->getTopologyKey();
358 
359  ScalarViewType workView("workView", numCells, projStruct->getMaxNumEvalPoints(ePointType), dim-1);
360 
361  if(numVertices>0) {
362  for(ordinal_type iv=0; iv<numVertices; ++iv) {
363  auto vertexEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(0,iv,ePointType));
364  RealSpaceTools<DeviceType>::clone(Kokkos::subview(ePoints, Kokkos::ALL(),
365  ePointsRange(0, iv), Kokkos::ALL()), vertexEPoints);
366  }
367  }
368 
369  for(ordinal_type ie=0; ie<numEdges; ++ie) {
370  auto edgePointsRange = ePointsRange(edgeDim, ie);
371  auto edgeEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(edgeDim,ie,ePointType));
372 
373  const auto topoKey = refTopologyKey(edgeDim,ie);
374  Kokkos::parallel_for
375  ("Evaluate Points",
376  Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
377  KOKKOS_LAMBDA (const size_t ic) {
378 
379  ordinal_type eOrt[12];
380  orts(ic).getEdgeOrientation(eOrt, numEdges);
381  ordinal_type ort = eOrt[ie];
382 
383  Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(ePoints,ic,edgePointsRange,Kokkos::ALL()),
384  edgeEPoints, subcellParamEdge, topoKey, ie, ort);
385  });
386  }
387 
388  for(ordinal_type iface=0; iface<numFaces; ++iface) {
389  auto facePointsRange = ePointsRange(faceDim, iface);
390  auto faceEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(faceDim,iface,ePointType));
391 
392  const auto topoKey = refTopologyKey(faceDim,iface);
393  Kokkos::parallel_for
394  ("Evaluate Points",
395  Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
396  KOKKOS_LAMBDA (const size_t ic) {
397  ordinal_type fOrt[6];
398  orts(ic).getFaceOrientation(fOrt, numFaces);
399  ordinal_type ort = fOrt[iface];
400 
401  Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(ePoints, ic, facePointsRange, Kokkos::ALL()),
402  faceEPoints, subcellParamFace, topoKey, iface, ort);
403  });
404  }
405 
406 
407  if(numVols > 0) {
408  auto pointsRange = ePointsRange(dim, 0);
409  auto cellEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,ePointType));
410  RealSpaceTools<DeviceType>::clone(Kokkos::subview(ePoints, Kokkos::ALL(), pointsRange, Kokkos::ALL()), cellEPoints);
411  }
412 }
413 
414 template<typename DeviceType>
415 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
416 typename funValsValueType, class ...funValsProperties,
417 typename BasisType,
418 typename ortValueType,class ...ortProperties>
419 void
420 ProjectionTools<DeviceType>::getL2BasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
421  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtTargetEPoints,
422  const typename BasisType::ScalarViewType targetEPoints,
423  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
424  const BasisType* cellBasis,
426 
427  typedef typename BasisType::scalarType scalarType;
428  typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
429  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
430  const auto cellTopo = cellBasis->getBaseCellTopology();
431  ordinal_type dim = cellTopo.getDimension();
432  ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
433  ordinal_type basisCardinality = cellBasis->getCardinality();
434  ordinal_type numCells = targetAtTargetEPoints.extent(0);
435  const ordinal_type edgeDim = 1;
436  const ordinal_type faceDim = 2;
437  const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
438 
439  ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
440  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
441  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
442 
443  ScalarViewType refEdgesVec("refEdgesVec", numEdges, dim);
444  ScalarViewType refFacesTangents("refFaceTangents", numFaces, dim, 2);
445  ScalarViewType refFacesNormal("refFaceNormal", numFaces, dim);
446 
447  ordinal_type numVertexDofs = numVertices;
448 
449  ordinal_type numEdgeDofs(0);
450  for(ordinal_type ie=0; ie<numEdges; ++ie)
451  numEdgeDofs += cellBasis->getDofCount(edgeDim,ie);
452 
453  ordinal_type numFaceDofs(0);
454  for(ordinal_type iface=0; iface<numFaces; ++iface)
455  numFaceDofs += cellBasis->getDofCount(faceDim,iface);
456 
457  Kokkos::View<ordinal_type*, DeviceType> computedDofs("computedDofs", numVertexDofs+numEdgeDofs+numFaceDofs);
458 
459  auto basisEPointsRange = projStruct->getBasisPointsRange();
460 
461  auto refTopologyKey = projStruct->getTopologyKey();
462 
463  ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
464  ScalarViewType basisEPoints("basisEPoints",numCells,numTotalBasisEPoints, dim);
465  getL2EvaluationPoints(basisEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
466 
467  auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal());
468 
469  ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim);
470  ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
471  {
472  if(fieldDim == 1) {
473  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints);
474  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints);
475  for(ordinal_type ic=0; ic<numCells; ++ic) {
476  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
477  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
478  }
479  OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(),
480  Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis);
481  OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(),
482  Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis);
483  }
484  else {
485  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints,fieldDim);
486  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints,fieldDim);
487  for(ordinal_type ic=0; ic<numCells; ++ic) {
488  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
489  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
490  }
491  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
492  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
493  }
494  }
495 
496  {
497  auto hostComputedDofs = Kokkos::create_mirror_view(computedDofs);
498  ordinal_type computedDofsCount = 0;
499  for(ordinal_type iv=0; iv<numVertices; ++iv)
500  hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(0, iv, 0);
501 
502  for(ordinal_type ie=0; ie<numEdges; ++ie) {
503  ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
504  for(ordinal_type i=0; i<edgeCardinality; ++i)
505  hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(edgeDim, ie, i);
506  }
507 
508  for(ordinal_type iface=0; iface<numFaces; ++iface) {
509  ordinal_type faceCardinality = cellBasis->getDofCount(faceDim,iface);
510  for(ordinal_type i=0; i<faceCardinality; ++i)
511  hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(faceDim, iface, i);
512  }
513  Kokkos::deep_copy(computedDofs,hostComputedDofs);
514  }
515 
516  bool isHGradBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HGRAD);
517  bool isHCurlBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HCURL);
518  bool isHDivBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HDIV);
519  ordinal_type faceDofDim = isHCurlBasis ? 2 : 1;
520  ScalarViewType edgeCoeff("edgeCoeff", fieldDim);
521 
522 
523  const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
524 
525  if(isHGradBasis) {
526 
527  auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(MemSpaceType(), projStruct->getTargetPointsRange());
528  typedef ComputeBasisCoeffsOnVertices_L2<decltype(basisCoeffs), decltype(tagToOrdinal), decltype(targetEPointsRange),
529  decltype(targetAtTargetEPoints), decltype(basisAtTargetEPoints)> functorType;
530  Kokkos::parallel_for(policy, functorType(basisCoeffs, tagToOrdinal, targetEPointsRange,
531  targetAtTargetEPoints, basisAtTargetEPoints, numVertices));
532  }
533 
534  auto targetEPointsRange = projStruct->getTargetPointsRange();
535  for(ordinal_type ie=0; ie<numEdges; ++ie) {
536  auto edgeVec = Kokkos::subview(refEdgesVec, ie, Kokkos::ALL());
537  //auto edgeVecHost = Kokkos::create_mirror_view(edgeVec);
538 
539  /*if(isHCurlBasis) {
540  CellTools<DeviceType>::getReferenceEdgeTangent(edgeVecHost,ie, cellTopo);
541  } else if(isHDivBasis) {
542  CellTools<DeviceType>::getReferenceSideNormal(edgeVecHost, ie, cellTopo);
543  } else {
544  edgeVecHost(0) = 1;
545  }
546  Kokkos::deep_copy(edgeVec,edgeVecHost);*/
547  if(isHCurlBasis) {
548  CellTools<DeviceType>::getReferenceEdgeTangent(edgeVec, ie, cellTopo);
549  } else if(isHDivBasis) {
550  CellTools<DeviceType>::getReferenceSideNormal(edgeVec, ie, cellTopo);
551  } else {
552  deep_copy(edgeVec, 1.0);
553  }
554 
555  ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
556  ordinal_type numBasisEPoints = range_size(basisEPointsRange(edgeDim, ie));
557  ordinal_type numTargetEPoints = range_size(targetEPointsRange(edgeDim, ie));
558 
559  ScalarViewType basisDofAtBasisEPoints("BasisDofAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
560  ScalarViewType tragetDofAtTargetEPoints("TargetDofAtTargetEPoints",numCells, numTargetEPoints);
561  ScalarViewType weightedBasisAtBasisEPoints("weightedTanBasisAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
562  ScalarViewType weightedBasisAtTargetEPoints("weightedTanBasisAtTargetEPoints",numCells,edgeCardinality, numTargetEPoints);
563  ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints);
564 
565  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(edgeDim,ie));
566  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(edgeDim,ie));
567 
568  //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection
569  ordinal_type offsetBasis = basisEPointsRange(edgeDim, ie).first;
570  ordinal_type offsetTarget = targetEPointsRange(edgeDim, ie).first;
571 
572 
573  typedef ComputeBasisCoeffsOnEdges_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
574  decltype(computedDofs), decltype(tagToOrdinal), decltype(targetAtTargetEPoints)> functorTypeEdge;
575 
576  Kokkos::parallel_for(policy, functorTypeEdge(basisCoeffs, negPartialProj, basisDofAtBasisEPoints,
577  basisAtBasisEPoints, basisEWeights, weightedBasisAtBasisEPoints, targetEWeights,
578  basisAtTargetEPoints, weightedBasisAtTargetEPoints, computedDofs, tagToOrdinal,
579  targetAtTargetEPoints,tragetDofAtTargetEPoints, refEdgesVec, fieldDim,
580  edgeCardinality, offsetBasis, offsetTarget, numVertexDofs, edgeDim, ie));
581 
582 
583  ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality, edgeCardinality),
584  edgeRhsMat_("rhsMat_", numCells, edgeCardinality);
585 
586  FunctionSpaceTools<DeviceType >::integrate(edgeMassMat_, basisDofAtBasisEPoints, weightedBasisAtBasisEPoints);
587  FunctionSpaceTools<DeviceType >::integrate(edgeRhsMat_, tragetDofAtTargetEPoints, weightedBasisAtTargetEPoints);
588  FunctionSpaceTools<DeviceType >::integrate(edgeRhsMat_, negPartialProj, weightedBasisAtBasisEPoints,true);
589 
590 
591  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
592  ScalarViewType t_("t",numCells, edgeCardinality);
593  WorkArrayViewType w_("w",numCells, edgeCardinality);
594 
595  auto edgeDof = Kokkos::subview(tagToOrdinal, edgeDim, ie, Kokkos::ALL());
596 
597  ElemSystem edgeSystem("edgeSystem", false);
598  edgeSystem.solve(basisCoeffs, edgeMassMat_, edgeRhsMat_, t_, w_, edgeDof, edgeCardinality);
599  }
600 
601  typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamFace;
602  if(numFaces>0)
603  subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellBasis->getBaseCellTopology().getKey());
604 
605  ScalarViewType faceCoeff("faceCoeff", numCells, fieldDim, faceDofDim);
606  for(ordinal_type iface=0; iface<numFaces; ++iface) {
607  const auto topoKey = refTopologyKey(faceDim,iface);
608  ordinal_type faceCardinality = cellBasis->getDofCount(faceDim,iface);
609 
610  ordinal_type numTargetEPoints = range_size(targetEPointsRange(faceDim, iface));
611  ordinal_type numBasisEPoints = range_size(basisEPointsRange(faceDim, iface));
612 
613  ScalarViewType faceBasisDofAtBasisEPoints("normaBasisAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim);
614  ScalarViewType wBasisDofAtBasisEPoints("weightedNormalBasisAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim);
615 
616  ScalarViewType faceBasisAtTargetEPoints("normalBasisAtTargetEPoints",numCells,faceCardinality, numTargetEPoints,faceDofDim);
617  ScalarViewType wBasisDofAtTargetEPoints("weightedNormalBasisAtTargetEPoints",numCells,faceCardinality, numTargetEPoints,faceDofDim);
618 
619  ScalarViewType targetDofAtTargetEPoints("targetDofAtTargetEPoints",numCells, numTargetEPoints,faceDofDim);
620  ScalarViewType negPartialProj("mNormalComputedProjection", numCells,numBasisEPoints,faceDofDim);
621 
622  ordinal_type offsetBasis = basisEPointsRange(faceDim, iface).first;
623  ordinal_type offsetTarget = targetEPointsRange(faceDim, iface).first;
624  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(faceDim,iface));
625  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(faceDim,iface));
626 
627 
628  typedef ComputeBasisCoeffsOnFaces_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
629  decltype(computedDofs), decltype(tagToOrdinal), decltype(orts), decltype(targetAtTargetEPoints), decltype(subcellParamFace)> functorTypeFace;
630 
631  Kokkos::parallel_for(policy, functorTypeFace(basisCoeffs, negPartialProj,faceBasisDofAtBasisEPoints,
632  basisAtBasisEPoints, basisEWeights, wBasisDofAtBasisEPoints, targetEWeights,
633  basisAtTargetEPoints, wBasisDofAtTargetEPoints, computedDofs, tagToOrdinal,
634  orts, targetAtTargetEPoints,targetDofAtTargetEPoints, faceCoeff,
635  subcellParamFace, fieldDim, faceCardinality, offsetBasis,
636  offsetTarget, numVertexDofs+numEdgeDofs, numFaces, faceDim,faceDofDim,
637  dim, iface, topoKey, isHCurlBasis, isHDivBasis));
638 
639  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
640  ScalarViewType faceMassMat_("faceMassMat_", numCells, faceCardinality, faceCardinality),
641  faceRhsMat_("rhsMat_", numCells, faceCardinality);
642 
643  FunctionSpaceTools<DeviceType >::integrate(faceMassMat_, faceBasisDofAtBasisEPoints, wBasisDofAtBasisEPoints);
644  FunctionSpaceTools<DeviceType >::integrate(faceRhsMat_, targetDofAtTargetEPoints, wBasisDofAtTargetEPoints);
645  FunctionSpaceTools<DeviceType >::integrate(faceRhsMat_, negPartialProj, wBasisDofAtBasisEPoints,true);
646 
647  ScalarViewType t_("t",numCells, faceCardinality);
648  WorkArrayViewType w_("w",numCells,faceCardinality);
649 
650  auto faceDof = Kokkos::subview(tagToOrdinal, faceDim, iface, Kokkos::ALL());
651 
652  ElemSystem faceSystem("faceSystem", false);
653  faceSystem.solve(basisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDof, faceCardinality);
654  }
655 
656  ordinal_type numElemDofs = cellBasis->getDofCount(dim,0);
657 
658 
659  if(numElemDofs>0) {
660 
661  auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
662 
663  range_type cellPointsRange = targetEPointsRange(dim, 0);
664 
665  ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
666  ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
667 
668  ScalarViewType internalBasisAtBasisEPoints("internalBasisAtBasisEPoints",numCells,numElemDofs, numBasisEPoints, fieldDim);
669  ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints, fieldDim);
670 
671  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
672  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
673  ordinal_type offsetBasis = basisEPointsRange(dim, 0).first;
674  ordinal_type offsetTarget = targetEPointsRange(dim, 0).first;
675 
676 
677  ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",numCells,numElemDofs, numBasisEPoints,fieldDim);
678  ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTargetEPoints,fieldDim);
679 
681  Kokkos::parallel_for(policy, functorType( basisCoeffs, negPartialProj, internalBasisAtBasisEPoints,
682  basisAtBasisEPoints, basisEWeights, wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints,
683  computedDofs, cellDofs, fieldDim, numElemDofs, offsetBasis, offsetTarget, numVertexDofs+numEdgeDofs+numFaceDofs));
684 
685  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
686  ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs),
687  cellRhsMat_("rhsMat_", numCells, numElemDofs);
688 
689  FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, internalBasisAtBasisEPoints, wBasisAtBasisEPoints);
690  if(fieldDim==1)
691  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()),
692  Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
693  else
694  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), wBasisDofAtTargetEPoints);
695  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, negPartialProj, wBasisAtBasisEPoints, true);
696 
697  ScalarViewType t_("t",numCells, numElemDofs);
698  WorkArrayViewType w_("w",numCells,numElemDofs);
699  ElemSystem cellSystem("cellSystem", false);
700  cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
701  }
702 }
703 
704 template<typename ViewType1, typename ViewType2>
706  const ViewType1 basisAtBasisEPoints_;
707  const ViewType2 basisEWeights_;
708  const ViewType1 wBasisAtBasisEPoints_;
709  const ViewType2 targetEWeights_;
710  const ViewType1 basisAtTargetEPoints_;
711  const ViewType1 wBasisDofAtTargetEPoints_;
712  ordinal_type fieldDim_;
713  ordinal_type numElemDofs_;
714 
715  MultiplyBasisByWeights(const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisAtBasisEPoints, const ViewType2 targetEWeights,
716  const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEPoints,
717  ordinal_type fieldDim, ordinal_type numElemDofs) :
718  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
719  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
720  fieldDim_(fieldDim), numElemDofs_(numElemDofs) {}
721 
722  void
723  KOKKOS_INLINE_FUNCTION
724  operator()(const ordinal_type ic) const {
725 
726  for(ordinal_type j=0; j <numElemDofs_; ++j) {
727  for(ordinal_type d=0; d <fieldDim_; ++d) {
728  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
729  wBasisAtBasisEPoints_(ic,j,iq,d) = basisAtBasisEPoints_(ic,j,iq,d) * basisEWeights_(iq);
730  }
731  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
732  wBasisDofAtTargetEPoints_(ic,j,iq,d) = basisAtTargetEPoints_(ic,j,iq,d)* targetEWeights_(iq);
733  }
734  }
735  }
736  }
737 };
738 
739 template<typename DeviceType>
740 template<typename BasisType>
741 void
742 ProjectionTools<DeviceType>::getL2DGEvaluationPoints(typename BasisType::ScalarViewType ePoints,
743  const BasisType* cellBasis,
745  const EvalPointsType ePointType) {
746 
747  ordinal_type dim = cellBasis->getBaseCellTopology().getDimension();
748  auto cellEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,ePointType));
749  RealSpaceTools<DeviceType>::clone(ePoints, cellEPoints);
750 }
751 
752 template<typename DeviceType>
753 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
754 typename funValsValueType, class ...funValsProperties,
755 typename BasisType,
756 typename ortValueType,class ...ortProperties>
757 void
758 ProjectionTools<DeviceType>::getL2DGBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
759  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtTargetEPoints,
760  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
761  const BasisType* cellBasis,
763 
764  typedef typename BasisType::scalarType scalarType;
765  typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
766  const auto cellTopo = cellBasis->getBaseCellTopology();
767 
768  ordinal_type dim = cellTopo.getDimension();
769 
770  auto basisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
771  projStruct->getEvalPoints(dim,0,EvalPointsType::BASIS));
772  auto targetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
773  projStruct->getEvalPoints(dim,0,EvalPointsType::TARGET));
774 
775 
776  ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
777  ordinal_type basisCardinality = cellBasis->getCardinality();
778  ordinal_type numCells = targetAtTargetEPoints.extent(0);
779  const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
780 
781  ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
782  ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim);
783  ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
784  {
785  if(fieldDim == 1) {
786  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints);
787  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints);
788  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
789  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
790 
791  RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtTargetEPoints, Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL()));
792  RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtBasisEPoints, Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL()));
793  OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(),
794  Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis);
795  OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(),
796  Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis);
797  }
798  else {
799  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints,fieldDim);
800  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints,fieldDim);
801  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
802  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
803 
804  RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtTargetEPoints, Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
805  RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtBasisEPoints, Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
806  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
807  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
808  }
809  }
810 
811  const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
812  ordinal_type numElemDofs = cellBasis->getCardinality();
813 
814  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
815  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
816 
817  ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",numCells,numElemDofs, numTotalBasisEPoints,fieldDim);
818  ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTotalTargetEPoints,fieldDim);
819 
821  Kokkos::parallel_for( "Multiply basis by weights", policy, functorType(basisAtBasisEPoints, basisEWeights,
822  wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints, fieldDim, numElemDofs));// )){
823 
824  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
825  ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs),
826  cellRhsMat_("rhsMat_", numCells, numElemDofs);
827 
828  FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, basisAtBasisEPoints, wBasisAtBasisEPoints);
829  if(fieldDim==1)
830  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints,
831  Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
832  else
833  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints, wBasisDofAtTargetEPoints);
834 
835  Kokkos::DynRankView<ordinal_type,Kokkos::HostSpace> hCellDofs("cellDoFs", numElemDofs);
836  for(ordinal_type i=0; i<numElemDofs; ++i) hCellDofs(i) = i;
837  auto cellDofs = Kokkos::create_mirror_view_and_copy(MemSpaceType(),hCellDofs);
838 
839  ScalarViewType t_("t",numCells, numElemDofs);
840  WorkArrayViewType w_("w",numCells,numElemDofs);
841  ElemSystem cellSystem("cellSystem", false);
842  cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
843 }
844 
845 template<typename DeviceType>
846 template<typename basisViewType, typename targetViewType, typename BasisType>
847 void
849  const targetViewType targetAtTargetEPoints,
850  const BasisType* cellBasis,
852 
853  typedef typename BasisType::scalarType scalarType;
854  typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
855  const auto cellTopo = cellBasis->getBaseCellTopology();
856 
857  ordinal_type dim = cellTopo.getDimension();
858 
859  auto basisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
860  projStruct->getEvalPoints(dim,0,EvalPointsType::BASIS));
861  auto targetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
862  projStruct->getEvalPoints(dim,0,EvalPointsType::TARGET));
863 
864  ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
865  ordinal_type basisCardinality = cellBasis->getCardinality();
866  ordinal_type numCells = targetAtTargetEPoints.extent(0);
867  const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
868 
869  ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
870  ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",1,basisCardinality, numTotalBasisEPoints, fieldDim);
871  ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
872  {
873  if(fieldDim == 1) {
874  cellBasis->getValues(Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),0), targetEPoints);
875  cellBasis->getValues(Kokkos::subview(basisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),0), basisEPoints);
876 
877  RealSpaceTools<DeviceType>::clone(basisAtTargetEPoints, Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
878  }
879  else {
880  cellBasis->getValues(Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
881  cellBasis->getValues(Kokkos::subview(basisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
882 
883  RealSpaceTools<DeviceType>::clone(basisAtTargetEPoints, Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
884  }
885  }
886 
887  const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
888  ordinal_type numElemDofs = cellBasis->getCardinality();
889 
890  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
891  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
892 
893  ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",1,numElemDofs, numTotalBasisEPoints,fieldDim);
894  ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTotalTargetEPoints,fieldDim);
895  Kokkos::DynRankView<ordinal_type, DeviceType> cellDofs("cellDoFs", numElemDofs);
896 
897  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,numElemDofs),
898  KOKKOS_LAMBDA (const int &j) {
899  for(ordinal_type d=0; d <fieldDim; ++d) {
900  for(ordinal_type iq=0; iq < numTotalBasisEPoints; ++iq)
901  wBasisAtBasisEPoints(0,j,iq,d) = basisAtBasisEPoints(0,j,iq,d) * basisEWeights(iq);
902  for(ordinal_type iq=0; iq <numTotalTargetEPoints; ++iq) {
903  wBasisDofAtTargetEPoints(0,j,iq,d) = basisAtTargetEPoints(0,j,iq,d)* targetEWeights(iq);
904  }
905  }
906  cellDofs(j) = j;
907  });
908  RealSpaceTools<DeviceType>::clone(wBasisDofAtTargetEPoints, Kokkos::subview(wBasisDofAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
909 
910  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
911  ScalarViewType cellMassMat_("cellMassMat_", 1, numElemDofs, numElemDofs),
912  cellRhsMat_("rhsMat_", numCells, numElemDofs);
913 
914  FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, basisAtBasisEPoints, wBasisAtBasisEPoints);
915  if(fieldDim==1)
916  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints,
917  Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
918  else
919  FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints, wBasisDofAtTargetEPoints);
920 
921  ScalarViewType t_("t",1, numElemDofs);
922  WorkArrayViewType w_("w",numCells,numElemDofs);
923  ElemSystem cellSystem("cellSystem", true);
924  cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
925 }
926 
927 }
928 }
929 
930 #endif
931 
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
static void getL2BasisCoeffs(Kokkos::DynRankView< basisCoeffsValueType, basisCoeffsProperties... > basisCoeffs, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetAtEvalPoints, const typename BasisType::ScalarViewType evaluationPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct)
Computes the basis coefficients of the L2 projection of the target function.
static void modifyBasisByOrientation(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input, const OrientationViewType orts, const BasisType *basis)
Modify basis due to orientation.
static void getL2DGBasisCoeffs(Kokkos::DynRankView< basisCoeffsValueType, basisCoeffsProperties... > basisCoeffs, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetAtEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct)
Computes evaluation points for local L2 projection for broken HGRAD HCURL HDIV and HVOL spaces...
static KOKKOS_INLINE_FUNCTION void mapSubcellCoordsToRefCell(coordsViewType cellCoords, const subcellCoordsViewType subCellCoords, const ParamViewType subcellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Maps points defined on the subCell manifold into the parent Cell.
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...
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties... > refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
const range_tag getBasisPointsRange() const
Returns the range tag of the basis evaluation points subcells.
static void getL2DGEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=EvalPointsType::TARGET)
Computes evaluation points for local L2 projection for broken HGRAD HCURL HDIV and HVOL spaces...
view_type getBasisEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis evaluation weights on a subcell.
ordinal_type getMaxNumEvalPoints(const EvalPointsType type) const
Returns the maximum number of evaluation points across all the subcells.
static void getReferenceSideNormal(Kokkos::DynRankView< refSideNormalValueType, refSideNormalProperties... > refSideNormal, const ordinal_type sideOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
Header file for the Intrepid2::FunctionSpaceTools class.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
Class to solve a square system A x = b on each cell A is expected to be saddle a point (KKT) matrix o...
ordinal_type getNumBasisEvalPoints()
Returns number of basis evaluation points.
Header file for Intrepid2::ArrayTools class providing utilities for array operations.
static void integrate(Kokkos::DynRankView< outputValueValueType, outputValueProperties... > outputValues, const Kokkos::DynRankView< leftValueValueType, leftValueProperties... > leftValues, const Kokkos::DynRankView< rightValueValueType, rightValueProperties... > rightValues, const bool sumInto=false)
Contracts leftValues and rightValues arrays on the point and possibly space dimensions and stores the...
const range_tag getTargetPointsRange() const
Returns the range of the target function evaluation points on subcells.
void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau, ViewType3 w, const ViewType4 elemDof, ordinal_type n, ordinal_type m=0)
Solve the system and returns the basis coefficients solve the system either using Kokkos Kernel QR or...
An helper class to compute the evaluation points and weights needed for performing projections...
static KOKKOS_INLINE_FUNCTION void getRefSideTangentsAndNormal(TanNormViewType tangentsAndNormal, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) tangents and normal of the side subCell The normals are only defined for side...
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
const key_tag getTopologyKey() const
Returns the key tag for subcells.
static void getL2EvaluationPoints(typename BasisType::ScalarViewType evaluationPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=EvalPointsType::TARGET)
Computes evaluation points for L2 projection.
view_type getEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the basis/target evaluation points on a subcell.
view_type getTargetEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function evaluation weights on a subcell.