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 ViewType8 faceParametrization_;
177 ordinal_type fieldDim_;
178 ordinal_type faceCardinality_;
179 ordinal_type offsetBasis_;
180 ordinal_type offsetTarget_;
181 ordinal_type numVertexEdgeDofs_;
182 ordinal_type numFaces_;
183 ordinal_type faceDim_;
184 ordinal_type dim_;
185 ordinal_type iface_;
186 unsigned topoKey_;
187 bool isHCurlBasis_, isHDivBasis_;
188
189 ComputeBasisCoeffsOnFaces_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 faceBasisDofAtBasisEPoints,
190 const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisDofAtBasisEPoints, const ViewType3 targetEWeights,
191 const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 tagToOrdinal,
192 const ViewType6 orts, const ViewType7 targetAtTargetEPoints, const ViewType2 targetDofAtTargetEPoints,
193 const ViewType8 faceParametrization, ordinal_type fieldDim, ordinal_type faceCardinality, ordinal_type offsetBasis,
194 ordinal_type offsetTarget, ordinal_type numVertexEdgeDofs, ordinal_type numFaces, ordinal_type faceDim,
195 ordinal_type dim, ordinal_type iface, unsigned topoKey, bool isHCurlBasis, bool isHDivBasis) :
196 basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), faceBasisDofAtBasisEPoints_(faceBasisDofAtBasisEPoints),
197 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights),
198 basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
199 computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), orts_(orts), targetAtTargetEPoints_(targetAtTargetEPoints),
200 targetDofAtTargetEPoints_(targetDofAtTargetEPoints), faceParametrization_(faceParametrization),
201 fieldDim_(fieldDim), faceCardinality_(faceCardinality), offsetBasis_(offsetBasis),
202 offsetTarget_(offsetTarget), numVertexEdgeDofs_(numVertexEdgeDofs), numFaces_(numFaces),
203 faceDim_(faceDim), dim_(dim), iface_(iface), topoKey_(topoKey),
204 isHCurlBasis_(isHCurlBasis), isHDivBasis_(isHDivBasis)
205 {}
206
207 void
208 KOKKOS_INLINE_FUNCTION
209 operator()(const ordinal_type ic) const {
210
211 ordinal_type fOrt[6];
212 orts_(ic).getFaceOrientation(fOrt, numFaces_);
213 ordinal_type ort = fOrt[iface_];
214
215 if(isHCurlBasis_) {
216 typename ViewType3::value_type data[2*3];
217 auto tangents = ViewType3(data, 2, dim_);
218 Impl::OrientationTools::getRefSubcellTangents(tangents,faceParametrization_,topoKey_,iface_,ort);
219 typename ViewType3::value_type tmp[2] = {};
220 for(ordinal_type j=0; j <faceCardinality_; ++j) {
221 ordinal_type jdof = tagToOrdinal_(faceDim_, iface_, j);
222 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
223 tmp[0] = tmp[1] = 0;
224 for(ordinal_type d=0; d <fieldDim_; ++d) {
225 tmp[0] += tangents(0,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
226 tmp[1] += tangents(1,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
227 }
228 // u \times n = t_0 (u \dot t_1) - t1 (u \dot t_0)
229 for(ordinal_type d=0; d <fieldDim_; ++d) {
230 faceBasisDofAtBasisEPoints_(ic,j,iq,d) = tangents(0,d)*tmp[1]-tangents(1,d)*tmp[0];
231 wBasisDofAtBasisEPoints_(ic,j,iq,d) = faceBasisDofAtBasisEPoints_(ic,j,iq,d) * basisEWeights_(iq);
232 }
233 }
234 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
235 tmp[0] = tmp[1] = 0;
236 for(ordinal_type d=0; d <fieldDim_; ++d) {
237 tmp[0] += tangents(0,d)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
238 tmp[1] += tangents(1,d)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
239 }
240 // u \times n = t_0 (u \dot t_1) - t1 (u \dot t_0)
241 for(ordinal_type d=0; d <fieldDim_; ++d)
242 wBasisDofAtTargetEPoints_(ic,j,iq,d) = (tangents(0,d)*tmp[1]-tangents(1,d)*tmp[0]) * targetEWeights_(iq);
243 }
244 }
245
246 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
247 tmp[0] = tmp[1] = 0;
248 for(ordinal_type d=0; d <fieldDim_; ++d) {
249 tmp[0] += tangents(0,d)*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
250 tmp[1] += tangents(1,d)*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
251 }
252 // u \times n = t_0 (u \dot t_1) - t1 (u \dot t_0)
253 for(ordinal_type d=0; d <fieldDim_; ++d)
254 targetDofAtTargetEPoints_(ic,iq,d) = (tangents(0,d)*tmp[1]-tangents(1,d)*tmp[0]);
255 }
256
257 for(ordinal_type j=0; j <numVertexEdgeDofs_; ++j) {
258 ordinal_type jdof = computedDofs_(j);
259 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
260 tmp[0] = tmp[1] = 0;
261 for(ordinal_type d=0; d <fieldDim_; ++d) {
262 tmp[0] += tangents(0,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
263 tmp[1] += tangents(1,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
264 }
265
266 // u \times n = t_0 (u \dot t_1) - t1 (u \dot t_0)
267 for(ordinal_type d=0; d <fieldDim_; ++d)
268 negPartialProj_(ic,iq,d) -= (tangents(0,d)*tmp[1]-tangents(1,d)*tmp[0])*basisCoeffs_(ic,jdof);
269 }
270 }
271 } else {
272 typename ViewType3::value_type coeff[3];
273 if (isHDivBasis_) {
274 typename ViewType3::value_type data[3*3];
275 auto tangentsAndNormal = ViewType3(data, dim_, dim_);
276 Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,topoKey_, iface_, ort);
277 for(ordinal_type d=0; d <dim_; ++d)
278 coeff[d] = tangentsAndNormal(dim_-1,d);
279 }
280 else //isHGradBasis
281 coeff[0] = 1;
282
283 for(ordinal_type j=0; j <faceCardinality_; ++j) {
284 ordinal_type jdof = tagToOrdinal_(faceDim_, iface_, j);
285 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
286 faceBasisDofAtBasisEPoints_(ic,j,iq,0) = 0;
287 for(ordinal_type d=0; d <fieldDim_; ++d)
288 faceBasisDofAtBasisEPoints_(ic,j,iq,0) += coeff[d]*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
289 wBasisDofAtBasisEPoints_(ic,j,iq,0) = faceBasisDofAtBasisEPoints_(ic,j,iq,0) * basisEWeights_(iq);
290 }
291 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
292 typename ViewType2::value_type sum=0;
293 for(ordinal_type d=0; d <fieldDim_; ++d)
294 sum += coeff[d]*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
295 wBasisDofAtTargetEPoints_(ic,j,iq,0) = sum * targetEWeights_(iq);
296 }
297 }
298
299 for(ordinal_type d=0; d <fieldDim_; ++d)
300 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq)
301 targetDofAtTargetEPoints_(ic,iq,0) += coeff[d]*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
302
303 for(ordinal_type j=0; j <numVertexEdgeDofs_; ++j) {
304 ordinal_type jdof = computedDofs_(j);
305 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
306 for(ordinal_type d=0; d <fieldDim_; ++d)
307 negPartialProj_(ic,iq,0) -= basisCoeffs_(ic,jdof)*coeff[d]*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
308 }
309 }
310 }
311};
312
313
314template<typename ViewType1, typename ViewType2, typename ViewType3,
315typename ViewType4, typename ViewType5>
317 const ViewType1 basisCoeffs_;
318 const ViewType2 negPartialProj_;
319 const ViewType2 internalBasisAtBasisEPoints_;
320 const ViewType2 basisAtBasisEPoints_;
321 const ViewType3 basisEWeights_;
322 const ViewType2 wBasisAtBasisEPoints_;
323 const ViewType3 targetEWeights_;
324 const ViewType2 basisAtTargetEPoints_;
325 const ViewType2 wBasisDofAtTargetEPoints_;
326 const ViewType4 computedDofs_;
327 const ViewType5 elemDof_;
328 ordinal_type fieldDim_;
329 ordinal_type numElemDofs_;
330 ordinal_type offsetBasis_;
331 ordinal_type offsetTarget_;
332 ordinal_type numVertexEdgeFaceDofs_;
333
334 ComputeBasisCoeffsOnCells_L2(const ViewType1 basisCoeffs, ViewType2 negPartialProj, const ViewType2 internalBasisAtBasisEPoints,
335 const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisAtBasisEPoints, const ViewType3 targetEWeights,
336 const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisDofAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 elemDof,
337 ordinal_type fieldDim, ordinal_type numElemDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numVertexEdgeFaceDofs) :
338 basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), internalBasisAtBasisEPoints_(internalBasisAtBasisEPoints),
339 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
340 basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
341 computedDofs_(computedDofs), elemDof_(elemDof), fieldDim_(fieldDim), numElemDofs_(numElemDofs), offsetBasis_(offsetBasis),
342 offsetTarget_(offsetTarget), numVertexEdgeFaceDofs_(numVertexEdgeFaceDofs) {}
343
344 void
345 KOKKOS_INLINE_FUNCTION
346 operator()(const ordinal_type ic) const {
347
348 for(ordinal_type j=0; j <numElemDofs_; ++j) {
349 ordinal_type idof = elemDof_(j);
350 for(ordinal_type d=0; d <fieldDim_; ++d) {
351 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
352 internalBasisAtBasisEPoints_(ic,j,iq,d) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
353 wBasisAtBasisEPoints_(ic,j,iq,d) = internalBasisAtBasisEPoints_(ic,j,iq,d) * basisEWeights_(iq);
354 }
355 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
356 wBasisDofAtTargetEPoints_(ic,j,iq,d) = basisAtTargetEPoints_(ic,idof,offsetTarget_+iq,d)* targetEWeights_(iq);
357 }
358 }
359 }
360 for(ordinal_type j=0; j < numVertexEdgeFaceDofs_; ++j) {
361 ordinal_type jdof = computedDofs_(j);
362 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
363 for(ordinal_type d=0; d <fieldDim_; ++d) {
364 negPartialProj_(ic,iq,d) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
365 }
366 }
367 }
368};
369
370
371template<typename DeviceType>
372template<typename BasisType,
373typename ortValueType, class ...ortProperties>
374void
375ProjectionTools<DeviceType>::getL2EvaluationPoints(typename BasisType::ScalarViewType ePoints,
376 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
377 const BasisType* cellBasis,
379 const EvalPointsType ePointType) {
380 typedef typename BasisType::scalarType scalarType;
381 typedef Kokkos::DynRankView<scalarType,ortProperties...> ScalarViewType;
382 const auto cellTopo = cellBasis->getBaseCellTopology();
383 //const auto cellTopoKey = cellBasis->getBaseCellTopology().getKey();
384 ordinal_type dim = cellTopo.getDimension();
385 ordinal_type numCells = ePoints.extent(0);
386 const ordinal_type edgeDim = 1;
387 const ordinal_type faceDim = 2;
388
389 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
390 ordinal_type numEdges = (cellBasis->getDofCount(edgeDim, 0) > 0) ? cellTopo.getEdgeCount() : 0;
391 ordinal_type numFaces = (cellBasis->getDofCount(faceDim, 0) > 0) ? cellTopo.getFaceCount() : 0;
392 ordinal_type numVols = (cellBasis->getDofCount(dim, 0) > 0);
393
394 auto ePointsRange = projStruct->getPointsRange(ePointType);
395
396 typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamEdge, subcellParamFace;
397 if(numEdges>0)
398 subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
399 if(numFaces>0)
400 subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
401
402 auto refTopologyKey = projStruct->getTopologyKey();
403
404 ScalarViewType workView("workView", numCells, projStruct->getMaxNumEvalPoints(ePointType), dim-1);
405
406 if(numVertices>0) {
407 for(ordinal_type iv=0; iv<numVertices; ++iv) {
408 auto vertexEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(0,iv,ePointType));
409 RealSpaceTools<DeviceType>::clone(Kokkos::subview(ePoints, Kokkos::ALL(),
410 ePointsRange(0, iv), Kokkos::ALL()), vertexEPoints);
411 }
412 }
413
414 for(ordinal_type ie=0; ie<numEdges; ++ie) {
415 auto edgePointsRange = ePointsRange(edgeDim, ie);
416 auto edgeEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(edgeDim,ie,ePointType));
417
418 const auto topoKey = refTopologyKey(edgeDim,ie);
419 Kokkos::parallel_for
420 ("Evaluate Points",
421 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
422 KOKKOS_LAMBDA (const size_t ic) {
423
424 ordinal_type eOrt[12];
425 orts(ic).getEdgeOrientation(eOrt, numEdges);
426 ordinal_type ort = eOrt[ie];
427
428 Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(ePoints,ic,edgePointsRange,Kokkos::ALL()),
429 edgeEPoints, subcellParamEdge, topoKey, ie, ort);
430 });
431 }
432
433 for(ordinal_type iface=0; iface<numFaces; ++iface) {
434 auto facePointsRange = ePointsRange(faceDim, iface);
435 auto faceEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(faceDim,iface,ePointType));
436
437 const auto topoKey = refTopologyKey(faceDim,iface);
438 Kokkos::parallel_for
439 ("Evaluate Points",
440 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
441 KOKKOS_LAMBDA (const size_t ic) {
442 ordinal_type fOrt[6];
443 orts(ic).getFaceOrientation(fOrt, numFaces);
444 ordinal_type ort = fOrt[iface];
445
446 Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(ePoints, ic, facePointsRange, Kokkos::ALL()),
447 faceEPoints, subcellParamFace, topoKey, iface, ort);
448 });
449 }
450
451
452 if(numVols > 0) {
453 auto pointsRange = ePointsRange(dim, 0);
454 auto cellEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,ePointType));
455 if(dim == 3)
456 RealSpaceTools<DeviceType>::clone(Kokkos::subview(ePoints, Kokkos::ALL(), pointsRange, Kokkos::ALL()), cellEPoints);
457 else { //if the cell is a side also internal points need orientation
458 const auto topoKey = refTopologyKey(dim,0);
459 RealSpaceTools<DeviceType>::clone(Kokkos::subview(ePoints, Kokkos::ALL(), pointsRange, Kokkos::ALL()), cellEPoints);
460 Kokkos::parallel_for
461 ("Evaluate Points",
462 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
463 KOKKOS_LAMBDA (const size_t ic) {
464 ordinal_type ort = 0;
465 if(dim == 1)
466 orts(ic).getEdgeOrientation(&ort,1);
467 else if (dim == 2)
468 orts(ic).getFaceOrientation(&ort,1);
469 Impl::OrientationTools::mapToModifiedReference(Kokkos::subview(ePoints, ic, pointsRange, Kokkos::ALL()), cellEPoints,topoKey,ort);
470 });
471 }
472 }
473}
474
475template<typename DeviceType>
476template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
477typename funValsValueType, class ...funValsProperties,
478typename BasisType,
479typename ortValueType,class ...ortProperties>
480void
481ProjectionTools<DeviceType>::getL2BasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
482 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtTargetEPoints,
483 const typename BasisType::ScalarViewType targetEPoints,
484 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
485 const BasisType* cellBasis,
487
488 typedef typename BasisType::scalarType scalarType;
489 typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
490 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
491 const auto cellTopo = cellBasis->getBaseCellTopology();
492 ordinal_type dim = cellTopo.getDimension();
493 ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
494 ordinal_type basisCardinality = cellBasis->getCardinality();
495 ordinal_type numCells = targetAtTargetEPoints.extent(0);
496 const ordinal_type edgeDim = 1;
497 const ordinal_type faceDim = 2;
498 const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
499
500 ordinal_type numVertices = (cellBasis->getDofCount(0, 0) > 0) ? cellTopo.getVertexCount() : 0;
501 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
502 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
503
504 ScalarViewType refEdgesVec("refEdgesVec", numEdges, dim);
505 ScalarViewType refFacesTangents("refFaceTangents", numFaces, dim, 2);
506 ScalarViewType refFacesNormal("refFaceNormal", numFaces, dim);
507
508 ordinal_type numVertexDofs = numVertices;
509
510 ordinal_type numEdgeDofs(0);
511 for(ordinal_type ie=0; ie<numEdges; ++ie)
512 numEdgeDofs += cellBasis->getDofCount(edgeDim,ie);
513
514 ordinal_type numFaceDofs(0);
515 for(ordinal_type iface=0; iface<numFaces; ++iface)
516 numFaceDofs += cellBasis->getDofCount(faceDim,iface);
517
518 Kokkos::View<ordinal_type*, DeviceType> computedDofs("computedDofs", numVertexDofs+numEdgeDofs+numFaceDofs);
519
520 auto basisEPointsRange = projStruct->getBasisPointsRange();
521
522 auto refTopologyKey = projStruct->getTopologyKey();
523
524 ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
525 ScalarViewType basisEPoints("basisEPoints",numCells,numTotalBasisEPoints, dim);
526 getL2EvaluationPoints(basisEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
527
528 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal());
529
530 ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim);
531 ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
532 {
533 if(fieldDim == 1) {
534 ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints);
535 ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints);
536 for(ordinal_type ic=0; ic<numCells; ++ic) {
537 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
538 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
539 }
540 OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(),
541 Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis);
542 OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(),
543 Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis);
544 }
545 else {
546 ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints,fieldDim);
547 ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints,fieldDim);
548 for(ordinal_type ic=0; ic<numCells; ++ic) {
549 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
550 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
551 }
552 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
553 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
554 }
555 }
556
557 {
558 auto hostComputedDofs = Kokkos::create_mirror_view(computedDofs);
559 ordinal_type computedDofsCount = 0;
560 for(ordinal_type iv=0; iv<numVertices; ++iv)
561 hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(0, iv, 0);
562
563 for(ordinal_type ie=0; ie<numEdges; ++ie) {
564 ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
565 for(ordinal_type i=0; i<edgeCardinality; ++i)
566 hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(edgeDim, ie, i);
567 }
568
569 for(ordinal_type iface=0; iface<numFaces; ++iface) {
570 ordinal_type faceCardinality = cellBasis->getDofCount(faceDim,iface);
571 for(ordinal_type i=0; i<faceCardinality; ++i)
572 hostComputedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(faceDim, iface, i);
573 }
574 Kokkos::deep_copy(computedDofs,hostComputedDofs);
575 }
576
577 bool isHGradBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HGRAD);
578 bool isHCurlBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HCURL);
579 bool isHDivBasis = (cellBasis->getFunctionSpace() == FUNCTION_SPACE_HDIV);
580 ordinal_type faceDofDim = isHCurlBasis ? 3 : 1;
581 ScalarViewType edgeCoeff("edgeCoeff", fieldDim);
582
583
584 const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
585
586 if(isHGradBasis) {
587
588 auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(MemSpaceType(), projStruct->getTargetPointsRange());
589 typedef ComputeBasisCoeffsOnVertices_L2<decltype(basisCoeffs), decltype(tagToOrdinal), decltype(targetEPointsRange),
590 decltype(targetAtTargetEPoints), decltype(basisAtTargetEPoints)> functorType;
591 Kokkos::parallel_for(policy, functorType(basisCoeffs, tagToOrdinal, targetEPointsRange,
592 targetAtTargetEPoints, basisAtTargetEPoints, numVertices));
593 }
594
595 auto targetEPointsRange = projStruct->getTargetPointsRange();
596 for(ordinal_type ie=0; ie<numEdges; ++ie) {
597 auto edgeVec = Kokkos::subview(refEdgesVec, ie, Kokkos::ALL());
598 //auto edgeVecHost = Kokkos::create_mirror_view(edgeVec);
599
600 if(isHCurlBasis) {
602 } else if(isHDivBasis) {
604 } else {
605 deep_copy(edgeVec, 1.0);
606 }
607
608 ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
609 ordinal_type numBasisEPoints = range_size(basisEPointsRange(edgeDim, ie));
610 ordinal_type numTargetEPoints = range_size(targetEPointsRange(edgeDim, ie));
611
612 ScalarViewType basisDofAtBasisEPoints("BasisDofAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
613 ScalarViewType tragetDofAtTargetEPoints("TargetDofAtTargetEPoints",numCells, numTargetEPoints);
614 ScalarViewType weightedBasisAtBasisEPoints("weightedTanBasisAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
615 ScalarViewType weightedBasisAtTargetEPoints("weightedTanBasisAtTargetEPoints",numCells,edgeCardinality, numTargetEPoints);
616 ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints);
617
618 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(edgeDim,ie));
619 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(edgeDim,ie));
620
621 //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
622 ordinal_type offsetBasis = basisEPointsRange(edgeDim, ie).first;
623 ordinal_type offsetTarget = targetEPointsRange(edgeDim, ie).first;
624
625
626 typedef ComputeBasisCoeffsOnEdges_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
627 decltype(computedDofs), decltype(tagToOrdinal), decltype(targetAtTargetEPoints)> functorTypeEdge;
628
629 Kokkos::parallel_for(policy, functorTypeEdge(basisCoeffs, negPartialProj, basisDofAtBasisEPoints,
630 basisAtBasisEPoints, basisEWeights, weightedBasisAtBasisEPoints, targetEWeights,
631 basisAtTargetEPoints, weightedBasisAtTargetEPoints, computedDofs, tagToOrdinal,
632 targetAtTargetEPoints,tragetDofAtTargetEPoints, refEdgesVec, fieldDim,
633 edgeCardinality, offsetBasis, offsetTarget, numVertexDofs, edgeDim, ie));
634
635
636 ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality, edgeCardinality),
637 edgeRhsMat_("rhsMat_", numCells, edgeCardinality);
638
639 FunctionSpaceTools<DeviceType >::integrate(edgeMassMat_, basisDofAtBasisEPoints, weightedBasisAtBasisEPoints);
640 FunctionSpaceTools<DeviceType >::integrate(edgeRhsMat_, tragetDofAtTargetEPoints, weightedBasisAtTargetEPoints);
641 FunctionSpaceTools<DeviceType >::integrate(edgeRhsMat_, negPartialProj, weightedBasisAtBasisEPoints,true);
642
643
644 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
645 ScalarViewType t_("t",numCells, edgeCardinality);
646 WorkArrayViewType w_("w",numCells, edgeCardinality);
647
648 auto edgeDof = Kokkos::subview(tagToOrdinal, edgeDim, ie, Kokkos::ALL());
649
650 ElemSystem edgeSystem("edgeSystem", false);
651 edgeSystem.solve(basisCoeffs, edgeMassMat_, edgeRhsMat_, t_, w_, edgeDof, edgeCardinality);
652 }
653
654 typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamFace;
655 if(numFaces>0)
656 subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellBasis->getBaseCellTopology().getKey());
657
658 for(ordinal_type iface=0; iface<numFaces; ++iface) {
659 const auto topoKey = refTopologyKey(faceDim,iface);
660 ordinal_type faceCardinality = cellBasis->getDofCount(faceDim,iface);
661
662 ordinal_type numTargetEPoints = range_size(targetEPointsRange(faceDim, iface));
663 ordinal_type numBasisEPoints = range_size(basisEPointsRange(faceDim, iface));
664
665 ScalarViewType faceBasisDofAtBasisEPoints("faceBasisDofAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim);
666 ScalarViewType wBasisDofAtBasisEPoints("weightedBasisDofAtBasisEPoints",numCells,faceCardinality, numBasisEPoints,faceDofDim);
667
668 ScalarViewType faceBasisAtTargetEPoints("faceBasisDofAtTargetEPoints",numCells,faceCardinality, numTargetEPoints,faceDofDim);
669 ScalarViewType wBasisDofAtTargetEPoints("weightedBasisDofAtTargetEPoints",numCells,faceCardinality, numTargetEPoints,faceDofDim);
670
671 ScalarViewType targetDofAtTargetEPoints("targetDofAtTargetEPoints",numCells, numTargetEPoints,faceDofDim);
672 ScalarViewType negPartialProj("negComputedProjection", numCells,numBasisEPoints,faceDofDim);
673
674 ordinal_type offsetBasis = basisEPointsRange(faceDim, iface).first;
675 ordinal_type offsetTarget = targetEPointsRange(faceDim, iface).first;
676 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(faceDim,iface));
677 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(faceDim,iface));
678
679
680 typedef ComputeBasisCoeffsOnFaces_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
681 decltype(computedDofs), decltype(tagToOrdinal), decltype(orts), decltype(targetAtTargetEPoints), decltype(subcellParamFace)> functorTypeFace;
682
683 Kokkos::parallel_for(policy, functorTypeFace(basisCoeffs, negPartialProj,faceBasisDofAtBasisEPoints,
684 basisAtBasisEPoints, basisEWeights, wBasisDofAtBasisEPoints, targetEWeights,
685 basisAtTargetEPoints, wBasisDofAtTargetEPoints, computedDofs, tagToOrdinal,
686 orts, targetAtTargetEPoints,targetDofAtTargetEPoints,
687 subcellParamFace, fieldDim, faceCardinality, offsetBasis,
688 offsetTarget, numVertexDofs+numEdgeDofs, numFaces, faceDim,
689 dim, iface, topoKey, isHCurlBasis, isHDivBasis));
690
691 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
692 ScalarViewType faceMassMat_("faceMassMat_", numCells, faceCardinality, faceCardinality),
693 faceRhsMat_("rhsMat_", numCells, faceCardinality);
694
695 FunctionSpaceTools<DeviceType >::integrate(faceMassMat_, faceBasisDofAtBasisEPoints, wBasisDofAtBasisEPoints);
696 FunctionSpaceTools<DeviceType >::integrate(faceRhsMat_, targetDofAtTargetEPoints, wBasisDofAtTargetEPoints);
697 FunctionSpaceTools<DeviceType >::integrate(faceRhsMat_, negPartialProj, wBasisDofAtBasisEPoints,true);
698
699 ScalarViewType t_("t",numCells, faceCardinality);
700 WorkArrayViewType w_("w",numCells,faceCardinality);
701
702 auto faceDof = Kokkos::subview(tagToOrdinal, faceDim, iface, Kokkos::ALL());
703
704 ElemSystem faceSystem("faceSystem", false);
705 faceSystem.solve(basisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDof, faceCardinality);
706 }
707
708 ordinal_type numElemDofs = cellBasis->getDofCount(dim,0);
709
710
711 if(numElemDofs>0) {
712
713 auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
714
715 range_type cellPointsRange = targetEPointsRange(dim, 0);
716
717 ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
718 ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
719
720 ScalarViewType internalBasisAtBasisEPoints("internalBasisAtBasisEPoints",numCells,numElemDofs, numBasisEPoints, fieldDim);
721 ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints, fieldDim);
722
723 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
724 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
725 ordinal_type offsetBasis = basisEPointsRange(dim, 0).first;
726 ordinal_type offsetTarget = targetEPointsRange(dim, 0).first;
727
728
729 ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",numCells,numElemDofs, numBasisEPoints,fieldDim);
730 ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTargetEPoints,fieldDim);
731
732 typedef ComputeBasisCoeffsOnCells_L2<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights), decltype(computedDofs), decltype(cellDofs)> functorType;
733 Kokkos::parallel_for(policy, functorType( basisCoeffs, negPartialProj, internalBasisAtBasisEPoints,
734 basisAtBasisEPoints, basisEWeights, wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints,
735 computedDofs, cellDofs, fieldDim, numElemDofs, offsetBasis, offsetTarget, numVertexDofs+numEdgeDofs+numFaceDofs));
736
737 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
738 ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs),
739 cellRhsMat_("rhsMat_", numCells, numElemDofs);
740
741 FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, internalBasisAtBasisEPoints, wBasisAtBasisEPoints);
742 if(fieldDim==1)
743 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()),
744 Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
745 else
746 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), wBasisDofAtTargetEPoints);
747 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, negPartialProj, wBasisAtBasisEPoints, true);
748
749 ScalarViewType t_("t",numCells, numElemDofs);
750 WorkArrayViewType w_("w",numCells,numElemDofs);
751 ElemSystem cellSystem("cellSystem", false);
752 cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
753 }
754}
755
756template<typename ViewType1, typename ViewType2>
758 const ViewType1 basisAtBasisEPoints_;
759 const ViewType2 basisEWeights_;
760 const ViewType1 wBasisAtBasisEPoints_;
761 const ViewType2 targetEWeights_;
762 const ViewType1 basisAtTargetEPoints_;
763 const ViewType1 wBasisDofAtTargetEPoints_;
764 ordinal_type fieldDim_;
765 ordinal_type numElemDofs_;
766
767 MultiplyBasisByWeights(const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisAtBasisEPoints, const ViewType2 targetEWeights,
768 const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEPoints,
769 ordinal_type fieldDim, ordinal_type numElemDofs) :
770 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
771 basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEPoints),
772 fieldDim_(fieldDim), numElemDofs_(numElemDofs) {}
773
774 void
775 KOKKOS_INLINE_FUNCTION
776 operator()(const ordinal_type ic) const {
777
778 for(ordinal_type j=0; j <numElemDofs_; ++j) {
779 for(ordinal_type d=0; d <fieldDim_; ++d) {
780 for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
781 wBasisAtBasisEPoints_(ic,j,iq,d) = basisAtBasisEPoints_(ic,j,iq,d) * basisEWeights_(iq);
782 }
783 for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
784 wBasisDofAtTargetEPoints_(ic,j,iq,d) = basisAtTargetEPoints_(ic,j,iq,d)* targetEWeights_(iq);
785 }
786 }
787 }
788 }
789};
790
791template<typename DeviceType>
792template<typename BasisType>
793void
794ProjectionTools<DeviceType>::getL2DGEvaluationPoints(typename BasisType::ScalarViewType ePoints,
795 const BasisType* cellBasis,
797 const EvalPointsType ePointType) {
798
799 ordinal_type dim = cellBasis->getBaseCellTopology().getDimension();
800 auto cellEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,ePointType));
801 RealSpaceTools<DeviceType>::clone(ePoints, cellEPoints);
802}
803
804template<typename DeviceType>
805template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
806typename funValsValueType, class ...funValsProperties,
807typename BasisType,
808typename ortValueType,class ...ortProperties>
809void
810ProjectionTools<DeviceType>::getL2DGBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
811 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtTargetEPoints,
812 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
813 const BasisType* cellBasis,
815
816 typedef typename BasisType::scalarType scalarType;
817 typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
818 const auto cellTopo = cellBasis->getBaseCellTopology();
819
820 ordinal_type dim = cellTopo.getDimension();
821
822 auto basisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
823 projStruct->getEvalPoints(dim,0,EvalPointsType::BASIS));
824 auto targetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
825 projStruct->getEvalPoints(dim,0,EvalPointsType::TARGET));
826
827
828 ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
829 ordinal_type basisCardinality = cellBasis->getCardinality();
830 ordinal_type numCells = targetAtTargetEPoints.extent(0);
831 const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
832
833 ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
834 ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, fieldDim);
835 ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
836 {
837 if(fieldDim == 1) {
838 ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints);
839 ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints);
840 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
841 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
842
843 RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtTargetEPoints, Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL()));
844 RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtBasisEPoints, Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL()));
845 OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtBasisEPoints, Kokkos::ALL(), Kokkos::ALL(),
846 Kokkos::ALL(),0), nonOrientedBasisAtBasisEPoints, orts, cellBasis);
847 OrientationTools<DeviceType>::modifyBasisByOrientation(Kokkos::subview(basisAtTargetEPoints, Kokkos::ALL(),
848 Kokkos::ALL(), Kokkos::ALL(),0), nonOrientedBasisAtTargetEPoints, orts, cellBasis);
849 }
850 else {
851 ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints,fieldDim);
852 ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints,fieldDim);
853 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
854 cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
855
856 RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtTargetEPoints, Kokkos::subview(nonOrientedBasisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
857 RealSpaceTools<DeviceType>::clone(nonOrientedBasisAtBasisEPoints, Kokkos::subview(nonOrientedBasisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
858 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
859 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
860 }
861 }
862
863 const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
864 ordinal_type numElemDofs = cellBasis->getCardinality();
865
866 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
867 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
868
869 ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",numCells,numElemDofs, numTotalBasisEPoints,fieldDim);
870 ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTotalTargetEPoints,fieldDim);
871
872 typedef MultiplyBasisByWeights<decltype(basisAtBasisEPoints), decltype(basisEWeights)> functorType;
873 Kokkos::parallel_for( "Multiply basis by weights", policy, functorType(basisAtBasisEPoints, basisEWeights,
874 wBasisAtBasisEPoints, targetEWeights, basisAtTargetEPoints, wBasisDofAtTargetEPoints, fieldDim, numElemDofs));// )){
875
876 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
877 ScalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs, numElemDofs),
878 cellRhsMat_("rhsMat_", numCells, numElemDofs);
879
880 FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, basisAtBasisEPoints, wBasisAtBasisEPoints);
881 if(fieldDim==1)
882 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints,
883 Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
884 else
885 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints, wBasisDofAtTargetEPoints);
886
887 Kokkos::DynRankView<ordinal_type,Kokkos::HostSpace> hCellDofs("cellDoFs", numElemDofs);
888 for(ordinal_type i=0; i<numElemDofs; ++i) hCellDofs(i) = i;
889 auto cellDofs = Kokkos::create_mirror_view_and_copy(MemSpaceType(),hCellDofs);
890
891 ScalarViewType t_("t",numCells, numElemDofs);
892 WorkArrayViewType w_("w",numCells,numElemDofs);
893 ElemSystem cellSystem("cellSystem", false);
894 cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
895}
896
897template<typename DeviceType>
898template<typename basisViewType, typename targetViewType, typename BasisType>
899void
901 const targetViewType targetAtTargetEPoints,
902 const BasisType* cellBasis,
904
905 typedef typename BasisType::scalarType scalarType;
906 typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
907 const auto cellTopo = cellBasis->getBaseCellTopology();
908
909 ordinal_type dim = cellTopo.getDimension();
910
911 auto basisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
912 projStruct->getEvalPoints(dim,0,EvalPointsType::BASIS));
913 auto targetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),
914 projStruct->getEvalPoints(dim,0,EvalPointsType::TARGET));
915
916 ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1));
917 ordinal_type basisCardinality = cellBasis->getCardinality();
918 ordinal_type numCells = targetAtTargetEPoints.extent(0);
919 const ordinal_type fieldDim = (targetAtTargetEPoints.rank()==2) ? 1 : targetAtTargetEPoints.extent(2);
920
921 ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints();
922 ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",1,basisCardinality, numTotalBasisEPoints, fieldDim);
923 ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, fieldDim);
924 {
925 if(fieldDim == 1) {
926 cellBasis->getValues(Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),0), targetEPoints);
927 cellBasis->getValues(Kokkos::subview(basisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),0), basisEPoints);
928
929 RealSpaceTools<DeviceType>::clone(basisAtTargetEPoints, Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
930 }
931 else {
932 cellBasis->getValues(Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), targetEPoints);
933 cellBasis->getValues(Kokkos::subview(basisAtBasisEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), basisEPoints);
934
935 RealSpaceTools<DeviceType>::clone(basisAtTargetEPoints, Kokkos::subview(basisAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
936 }
937 }
938
939 const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
940 ordinal_type numElemDofs = cellBasis->getCardinality();
941
942 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
943 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
944
945 ScalarViewType wBasisAtBasisEPoints("weightedBasisAtBasisEPoints",1,numElemDofs, numTotalBasisEPoints,fieldDim);
946 ScalarViewType wBasisDofAtTargetEPoints("weightedBasisAtTargetEPoints",numCells,numElemDofs, numTotalTargetEPoints,fieldDim);
947 Kokkos::DynRankView<ordinal_type, DeviceType> cellDofs("cellDoFs", numElemDofs);
948
949 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,numElemDofs),
950 KOKKOS_LAMBDA (const int &j) {
951 for(ordinal_type d=0; d <fieldDim; ++d) {
952 for(ordinal_type iq=0; iq < numTotalBasisEPoints; ++iq)
953 wBasisAtBasisEPoints(0,j,iq,d) = basisAtBasisEPoints(0,j,iq,d) * basisEWeights(iq);
954 for(ordinal_type iq=0; iq <numTotalTargetEPoints; ++iq) {
955 wBasisDofAtTargetEPoints(0,j,iq,d) = basisAtTargetEPoints(0,j,iq,d)* targetEWeights(iq);
956 }
957 }
958 cellDofs(j) = j;
959 });
960 RealSpaceTools<DeviceType>::clone(wBasisDofAtTargetEPoints, Kokkos::subview(wBasisDofAtTargetEPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
961
962 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
963 ScalarViewType cellMassMat_("cellMassMat_", 1, numElemDofs, numElemDofs),
964 cellRhsMat_("rhsMat_", numCells, numElemDofs);
965
966 FunctionSpaceTools<DeviceType >::integrate(cellMassMat_, basisAtBasisEPoints, wBasisAtBasisEPoints);
967 if(fieldDim==1)
968 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints,
969 Kokkos::subview(wBasisDofAtTargetEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
970 else
971 FunctionSpaceTools<DeviceType >::integrate(cellRhsMat_, targetAtTargetEPoints, wBasisDofAtTargetEPoints);
972
973 ScalarViewType t_("t",1, numElemDofs);
974 WorkArrayViewType w_("w",numCells,numElemDofs);
975 ElemSystem cellSystem("cellSystem", true);
976 cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numElemDofs);
977}
978
979}
980}
981
982#endif
983
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 getRefSubcellTangents(TanViewType tangents, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) subCell tangents.
static void mapToModifiedReference(outPointViewType outPoints, const refPointViewType refPoints, const shards::CellTopology cellTopo, const ordinal_type cellOrt=0)
Computes modified parameterization maps of 1- and 2-subcells with orientation.
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...