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
55
56
57namespace Intrepid2 {
58namespace Experimental {
59
60
61template<typename ViewType1, typename ViewType2, typename ViewType3,
62typename 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
90template<typename ViewType1, typename ViewType2, typename ViewType3,
91typename 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
159template<typename ViewType1, typename ViewType2, typename ViewType3,
160typename 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
269template<typename ViewType1, typename ViewType2, typename ViewType3,
270typename 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
326template<typename DeviceType>
327template<typename BasisType,
328typename ortValueType, class ...ortProperties>
329void
330ProjectionTools<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
414template<typename DeviceType>
415template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
416typename funValsValueType, class ...funValsProperties,
417typename BasisType,
418typename ortValueType,class ...ortProperties>
419void
420ProjectionTools<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) {
549 } else if(isHDivBasis) {
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
680 typedef ComputeBasisCoeffsOnCells_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights), decltype(computedDofs), decltype(cellDofs)> functorType;
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
704template<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
739template<typename DeviceType>
740template<typename BasisType>
741void
742ProjectionTools<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
752template<typename DeviceType>
753template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
754typename funValsValueType, class ...funValsProperties,
755typename BasisType,
756typename ortValueType,class ...ortProperties>
757void
758ProjectionTools<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
820 typedef MultiplyBasisByWeights<decltype(basisAtBasisEPoints), decltype(basisEWeights)> functorType;
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
845template<typename DeviceType>
846template<typename basisViewType, typename targetViewType, typename BasisType>
847void
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 Intrepid2::ArrayTools class providing utilities for array operations.
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
Header file for the Intrepid2::FunctionSpaceTools class.
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.
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.
An helper class to compute the evaluation points and weights needed for performing projections.
ordinal_type getNumBasisEvalPoints()
Returns number of basis evaluation points.
const range_tag getBasisPointsRange() const
Returns the range tag of the basis evaluation points subcells.
const range_tag getTargetPointsRange() const
Returns the range of the target function evaluation points on subcells.
view_type getEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the basis/target evaluation points on a subcell.
const key_tag getTopologyKey() const
Returns the key tag for subcells.
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.
view_type getTargetEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function evaluation weights on a subcell.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
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 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.
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.
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 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...
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 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 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 clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
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...
Class to solve a square system A x = b on each cell A is expected to be saddle a point (KKT) matrix o...
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...