Intrepid2
Intrepid2_ProjectionToolsDefHCURL.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
48
49#ifndef __INTREPID2_PROJECTIONTOOLSDEFHCURL_HPP__
50#define __INTREPID2_PROJECTIONTOOLSDEFHCURL_HPP__
51
55
56
57namespace Intrepid2 {
58namespace Experimental {
59
60template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4>
61struct ComputeBasisCoeffsOnEdges_HCurl {
62 const ViewType1 basisTanAtBasisEPoints_;
63 const ViewType1 basisAtBasisEPoints_;
64 const ViewType2 basisEWeights_;
65 const ViewType1 wTanBasisAtBasisEPoints_;
66 const ViewType2 targetEWeights_;
67 const ViewType1 basisAtTargetEPoints_;
68 const ViewType1 wTanBasisAtTargetEPoints_;
69 const ViewType3 tagToOrdinal_;
70 const ViewType4 targetAtTargetEPoints_;
71 const ViewType1 targetTanAtTargetEPoints_;
72 const ViewType1 refEdgesTangent_;
73 ordinal_type edgeCardinality_;
74 ordinal_type offsetBasis_;
75 ordinal_type offsetTarget_;
76 ordinal_type edgeDim_;
77 ordinal_type dim_;
78 ordinal_type iedge_;
79
80 ComputeBasisCoeffsOnEdges_HCurl(const ViewType1 basisTanAtBasisEPoints,
81 const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wTanBasisAtBasisEPoints, const ViewType2 targetEWeights,
82 const ViewType1 basisAtTargetEPoints, const ViewType1 wTanBasisAtTargetEPoints, const ViewType3 tagToOrdinal,
83 const ViewType4 targetAtTargetEPoints, const ViewType1 targetTanAtTargetEPoints,
84 const ViewType1 refEdgesTangent, ordinal_type edgeCardinality, ordinal_type offsetBasis,
85 ordinal_type offsetTarget, ordinal_type edgeDim,
86 ordinal_type dim, ordinal_type iedge) :
87 basisTanAtBasisEPoints_(basisTanAtBasisEPoints),
88 basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wTanBasisAtBasisEPoints_(wTanBasisAtBasisEPoints), targetEWeights_(targetEWeights),
89 basisAtTargetEPoints_(basisAtTargetEPoints), wTanBasisAtTargetEPoints_(wTanBasisAtTargetEPoints),
90 tagToOrdinal_(tagToOrdinal), targetAtTargetEPoints_(targetAtTargetEPoints),
91 targetTanAtTargetEPoints_(targetTanAtTargetEPoints),
92 refEdgesTangent_(refEdgesTangent), edgeCardinality_(edgeCardinality), offsetBasis_(offsetBasis),
93 offsetTarget_(offsetTarget), edgeDim_(edgeDim), dim_(dim), iedge_(iedge)
94 {}
95
96 void
97 KOKKOS_INLINE_FUNCTION
98 operator()(const ordinal_type ic) const {
99
100 ordinal_type numBasisEPoints = basisEWeights_.extent(0);
101 ordinal_type numTargetEPoints = targetEWeights_.extent(0);
102 for(ordinal_type j=0; j <edgeCardinality_; ++j) {
103 ordinal_type jdof = tagToOrdinal_(edgeDim_, iedge_, j);
104 for(ordinal_type iq=0; iq <numBasisEPoints; ++iq) {
105 for(ordinal_type d=0; d <dim_; ++d)
106 basisTanAtBasisEPoints_(ic,j,iq) += refEdgesTangent_(iedge_,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
107 wTanBasisAtBasisEPoints_(ic,j,iq) = basisTanAtBasisEPoints_(ic,j,iq)*basisEWeights_(iq);
108 }
109
110 for(ordinal_type iq=0; iq <numTargetEPoints; ++iq) {
111 typename ViewType2::value_type tmp = 0;
112 for(ordinal_type d=0; d <dim_; ++d)
113 tmp += refEdgesTangent_(iedge_,d)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
114 wTanBasisAtTargetEPoints_(ic,j,iq) = tmp*targetEWeights_(iq);
115 }
116 }
117 for(ordinal_type iq=0; iq <numTargetEPoints; ++iq)
118 for(ordinal_type d=0; d <dim_; ++d)
119 targetTanAtTargetEPoints_(ic,iq) += refEdgesTangent_(iedge_,d)*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
120 }
121};
122
123
124template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4,
125typename ViewType5, typename ViewType6, typename ViewType7, typename ViewType8>
126struct ComputeBasisCoeffsOnFaces_HCurl {
127 const ViewType1 basisCoeffs_;
128 const ViewType2 orts_;
129 const ViewType3 negPartialProjTan_;
130 const ViewType3 negPartialProjCurlNormal_;
131 const ViewType3 hgradBasisGradAtBasisEPoints_;
132 const ViewType3 wHgradBasisGradAtBasisEPoints_;
133 const ViewType3 basisCurlAtBasisCurlEPoints_;
134 const ViewType3 basisCurlNormalAtBasisCurlEPoints_;
135 const ViewType3 basisAtBasisEPoints_;
136 const ViewType3 normalTargetCurlAtTargetEPoints_;
137 const ViewType3 basisTanAtBasisEPoints_;
138 const ViewType3 hgradBasisGradAtTargetEPoints_;
139 const ViewType3 wHgradBasisGradAtTargetEPoints_;
140 const ViewType3 wNormalBasisCurlAtBasisCurlEPoints_;
141 const ViewType3 basisCurlAtTargetCurlEPoints_;
142 const ViewType3 wNormalBasisCurlBasisAtTargetCurlEPoints_;
143 const ViewType4 targetAtTargetEPoints_;
144 const ViewType3 targetTanAtTargetEPoints_;
145 const ViewType4 targetCurlAtTargetCurlEPoints_;
146 const ViewType5 basisEWeights_;
147 const ViewType5 targetEWeights_;
148 const ViewType5 basisCurlEWeights_;
149 const ViewType5 targetCurlEWeights_;
150 const ViewType6 tagToOrdinal_;
151 const ViewType6 hGradTagToOrdinal_;
152 const ViewType7 faceParametrization_;
153 const ViewType8 computedDofs_;
154 unsigned refTopologyKey_;
155 ordinal_type offsetBasis_;
156 ordinal_type offsetBasisCurl_;
157 ordinal_type offsetTarget_;
158 ordinal_type offsetTargetCurl_;
159 ordinal_type iface_;
160 ordinal_type hgradCardinality_;
161 ordinal_type numFaces_;
162 ordinal_type numFaceDofs_;
163 ordinal_type numEdgeDofs_;
164 ordinal_type faceDim_;
165 ordinal_type dim_;
166
167
168
169
170 ComputeBasisCoeffsOnFaces_HCurl(const ViewType1 basisCoeffs,
171 const ViewType2 orts, const ViewType3 negPartialProjTan, const ViewType3 negPartialProjCurlNormal,
172 const ViewType3 hgradBasisGradAtBasisEPoints, const ViewType3 wHgradBasisGradAtBasisEPoints,
173 const ViewType3 basisCurlAtBasisCurlEPoints, const ViewType3 basisCurlNormalAtBasisCurlEPoints,
174 const ViewType3 basisAtBasisEPoints,
175 const ViewType3 normalTargetCurlAtTargetEPoints,
176 const ViewType3 basisTanAtBasisEPoints,
177 const ViewType3 hgradBasisGradAtTargetEPoints, const ViewType3 wHgradBasisGradAtTargetEPoints,
178 const ViewType3 wNormalBasisCurlAtBasisCurlEPoints, const ViewType3 basisCurlAtTargetCurlEPoints,
179 const ViewType3 wNormalBasisCurlBasisAtTargetCurlEPoints, const ViewType4 targetAtTargetEPoints,
180 const ViewType3 targetTanAtTargetEPoints, const ViewType4 targetCurlAtTargetCurlEPoints,
181 const ViewType5 basisEWeights, const ViewType5 targetEWeights,
182 const ViewType5 basisCurlEWeights, const ViewType5 targetCurlEWeights, const ViewType6 tagToOrdinal,
183 const ViewType6 hGradTagToOrdinal, const ViewType7 faceParametrization,
184 const ViewType8 computedDofs, unsigned refTopologyKey, ordinal_type offsetBasis,
185 ordinal_type offsetBasisCurl, ordinal_type offsetTarget,
186 ordinal_type offsetTargetCurl, ordinal_type iface,
187 ordinal_type hgradCardinality, ordinal_type numFaces,
188 ordinal_type numFaceDofs, ordinal_type numEdgeDofs,
189 ordinal_type faceDim, ordinal_type dim):
190 basisCoeffs_(basisCoeffs),
191 orts_(orts), negPartialProjTan_(negPartialProjTan), negPartialProjCurlNormal_(negPartialProjCurlNormal),
192 hgradBasisGradAtBasisEPoints_(hgradBasisGradAtBasisEPoints), wHgradBasisGradAtBasisEPoints_(wHgradBasisGradAtBasisEPoints),
193 basisCurlAtBasisCurlEPoints_(basisCurlAtBasisCurlEPoints), basisCurlNormalAtBasisCurlEPoints_(basisCurlNormalAtBasisCurlEPoints),
194 basisAtBasisEPoints_(basisAtBasisEPoints),
195 normalTargetCurlAtTargetEPoints_(normalTargetCurlAtTargetEPoints), basisTanAtBasisEPoints_(basisTanAtBasisEPoints),
196 hgradBasisGradAtTargetEPoints_(hgradBasisGradAtTargetEPoints), wHgradBasisGradAtTargetEPoints_(wHgradBasisGradAtTargetEPoints),
197 wNormalBasisCurlAtBasisCurlEPoints_(wNormalBasisCurlAtBasisCurlEPoints), basisCurlAtTargetCurlEPoints_(basisCurlAtTargetCurlEPoints),
198 wNormalBasisCurlBasisAtTargetCurlEPoints_(wNormalBasisCurlBasisAtTargetCurlEPoints), targetAtTargetEPoints_(targetAtTargetEPoints),
199 targetTanAtTargetEPoints_(targetTanAtTargetEPoints), targetCurlAtTargetCurlEPoints_(targetCurlAtTargetCurlEPoints),
200 basisEWeights_(basisEWeights), targetEWeights_(targetEWeights),
201 basisCurlEWeights_(basisCurlEWeights), targetCurlEWeights_(targetCurlEWeights), tagToOrdinal_(tagToOrdinal),
202 hGradTagToOrdinal_(hGradTagToOrdinal), faceParametrization_(faceParametrization),
203 computedDofs_(computedDofs), refTopologyKey_(refTopologyKey), offsetBasis_(offsetBasis),
204 offsetBasisCurl_(offsetBasisCurl), offsetTarget_(offsetTarget),
205 offsetTargetCurl_(offsetTargetCurl), iface_(iface),
206 hgradCardinality_(hgradCardinality), numFaces_(numFaces),
207 numFaceDofs_(numFaceDofs), numEdgeDofs_(numEdgeDofs),
208 faceDim_(faceDim), dim_(dim){}
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
217 ordinal_type ort = fOrt[iface_];
218
219 typename ViewType3::value_type data[3*3];
220 auto tangentsAndNormal = ViewType3(data, dim_, dim_);
221
222 Impl::OrientationTools::getRefSideTangentsAndNormal(tangentsAndNormal, faceParametrization_,refTopologyKey_, iface_, ort);
223
224 ordinal_type numBasisEPoints = basisEWeights_.extent(0);
225 ordinal_type numTargetEPoints = targetEWeights_.extent(0);
226 for(ordinal_type j=0; j <hgradCardinality_; ++j) {
227 ordinal_type face_dof = hGradTagToOrdinal_(faceDim_, 0, j);
228 for(ordinal_type d=0; d <faceDim_; ++d) {
229 for(ordinal_type iq=0; iq <numBasisEPoints; ++iq)
230 wHgradBasisGradAtBasisEPoints_(ic, j, iq, d) = hgradBasisGradAtBasisEPoints_(face_dof,iq,d) * basisEWeights_(iq);
231 for(ordinal_type iq=0; iq <numTargetEPoints; ++iq)
232 wHgradBasisGradAtTargetEPoints_(ic,j,iq,d)= hgradBasisGradAtTargetEPoints_(face_dof,iq,d) * targetEWeights_(iq);
233 }
234 }
235
236 //Note: we are not considering the jacobian of the orientation map for normals since it is simply a scalar term for the integrals and it does not affect the projection
237 ordinal_type numBasisCurlEPoints = basisCurlEWeights_.extent(0);
238 for(ordinal_type j=0; j <numFaceDofs_; ++j) {
239 ordinal_type jdof = tagToOrdinal_(faceDim_, iface_, j);
240 for(ordinal_type iq=0; iq <numBasisEPoints; ++iq) {
241 for(ordinal_type d=0; d <dim_; ++d) {
242 for(ordinal_type itan=0; itan <dim_-1; ++itan)
243 basisTanAtBasisEPoints_(ic,j,iq,itan) += tangentsAndNormal(itan,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
244 }
245 }
246 for(ordinal_type iq=0; iq <numBasisCurlEPoints; ++iq) {
247 for(ordinal_type d=0; d <dim_; ++d)
248 basisCurlNormalAtBasisCurlEPoints_(ic,j,iq) += tangentsAndNormal(dim_-1,d)*basisCurlAtBasisCurlEPoints_(ic,jdof,offsetBasisCurl_+iq,d);
249 wNormalBasisCurlAtBasisCurlEPoints_(ic,j,iq) = basisCurlNormalAtBasisCurlEPoints_(ic,j,iq) * basisCurlEWeights_(iq);
250 }
251
252 ordinal_type numTargetCurlEPoints = targetCurlEWeights_.extent(0);
253 for(ordinal_type iq=0; iq <numTargetCurlEPoints; ++iq) {
254 typename ViewType3::value_type tmp=0;
255 for(ordinal_type d=0; d <dim_; ++d)
256 tmp += tangentsAndNormal(dim_-1,d)*basisCurlAtTargetCurlEPoints_(ic,jdof,offsetTargetCurl_+iq,d);
257 wNormalBasisCurlBasisAtTargetCurlEPoints_(ic,j,iq) = tmp*targetCurlEWeights_(iq);
258 }
259 }
260
261 for(ordinal_type j=0; j <numEdgeDofs_; ++j) {
262 ordinal_type jdof = computedDofs_(j);
263 for(ordinal_type iq=0; iq <numBasisEPoints; ++iq)
264 for(ordinal_type d=0; d <dim_; ++d) {
265 negPartialProjCurlNormal_(ic,iq) -= tangentsAndNormal(dim_-1,d)*basisCoeffs_(ic,jdof)*basisCurlAtBasisCurlEPoints_(ic,jdof,offsetBasisCurl_+iq,d);
266 for(ordinal_type itan=0; itan <dim_-1; ++itan)
267 negPartialProjTan_(ic,iq,itan) -= tangentsAndNormal(itan,d)*basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
268 }
269 }
270
271 ordinal_type numTargetCurlEPoints = targetCurlEWeights_.extent(0);
272 for(ordinal_type iq=0; iq <numTargetEPoints; ++iq)
273 for(ordinal_type d=0; d <dim_; ++d)
274 for(ordinal_type itan=0; itan <dim_-1; ++itan)
275 targetTanAtTargetEPoints_(ic,iq,itan) += tangentsAndNormal(itan,d)*targetAtTargetEPoints_(ic,offsetTarget_+iq,d);
276
277 for(ordinal_type iq=0; iq <numTargetCurlEPoints; ++iq)
278 for(ordinal_type d=0; d <dim_; ++d)
279 normalTargetCurlAtTargetEPoints_(ic,iq) += tangentsAndNormal(dim_-1,d)*targetCurlAtTargetCurlEPoints_(ic,offsetTargetCurl_+iq,d);
280 }
281};
282
283
284template<typename ViewType1, typename ViewType2, typename ViewType3,
285typename ViewType4, typename ViewType5>
286struct ComputeBasisCoeffsOnCell_HCurl {
287 const ViewType1 basisCoeffs_;
288 const ViewType2 negPartialProj_;
289 const ViewType2 negPartialProjCurl_;
290 const ViewType2 cellBasisAtBasisEPoints_;
291 const ViewType2 cellBasisCurlAtBasisCurlEPoints_;
292 const ViewType2 basisAtBasisEPoints_;
293 const ViewType2 hgradBasisGradAtBasisEPoints_;
294 const ViewType2 basisCurlAtBasisCurlEPoints_;
295 const ViewType2 hgradBasisGradAtTargetEPoints_;
296 const ViewType2 basisCurlAtTargetCurlEPoints_;
297 const ViewType3 basisEWeights_;
298 const ViewType3 basisCurlEWeights_;
299 const ViewType2 wHgradBasisGradAtBasisEPoints_;
300 const ViewType2 wBasisCurlAtBasisCurlEPoints_;
301 const ViewType3 targetEWeights_;
302 const ViewType3 targetCurlEWeights_;
303 const ViewType2 wHgradBasisGradAtTargetEPoints_;
304 const ViewType2 wBasisCurlAtTargetCurlEPoints_;
305 const ViewType4 computedDofs_;
306 const ViewType5 tagToOrdinal_;
307 const ViewType5 hGradTagToOrdinal_;
308 ordinal_type numCellDofs_;
309 ordinal_type hgradCardinality_;
310 ordinal_type offsetBasis_;
311 ordinal_type offsetBasisCurl_;
312 ordinal_type offsetTargetCurl_;
313 ordinal_type numEdgeFaceDofs_;
314 ordinal_type dim_;
315 ordinal_type derDim_;
316
317 ComputeBasisCoeffsOnCell_HCurl(const ViewType1 basisCoeffs, ViewType2 negPartialProj, ViewType2 negPartialProjCurl,
318 const ViewType2 cellBasisAtBasisEPoints, const ViewType2 cellBasisCurlAtBasisCurlEPoints,
319 const ViewType2 basisAtBasisEPoints, const ViewType2 hgradBasisGradAtBasisEPoints, const ViewType2 basisCurlAtBasisCurlEPoints,
320 const ViewType2 hgradBasisGradAtTargetEPoints, const ViewType2 basisCurlAtTargetCurlEPoints,
321 const ViewType3 basisEWeights, const ViewType3 basisCurlEWeights,
322 const ViewType2 wHgradBasisGradAtBasisEPoints, const ViewType2 wBasisCurlAtBasisCurlEPoints,
323 const ViewType3 targetEWeights, const ViewType3 targetCurlEWeights,
324 const ViewType2 wHgradBasisGradAtTargetEPoints,
325 const ViewType2 wBasisCurlAtTargetCurlEPoints, const ViewType4 computedDofs,
326 const ViewType5 tagToOrdinal, const ViewType5 hGradTagToOrdinal,
327 ordinal_type numCellDofs, ordinal_type hgradCardinality,
328 ordinal_type offsetBasis, ordinal_type offsetBasisCurl, ordinal_type offsetTargetCurl,
329 ordinal_type numEdgeFaceDofs, ordinal_type dim, ordinal_type derDim) :
330 basisCoeffs_(basisCoeffs), negPartialProj_(negPartialProj), negPartialProjCurl_(negPartialProjCurl),
331 cellBasisAtBasisEPoints_(cellBasisAtBasisEPoints), cellBasisCurlAtBasisCurlEPoints_(cellBasisCurlAtBasisCurlEPoints),
332 basisAtBasisEPoints_(basisAtBasisEPoints), hgradBasisGradAtBasisEPoints_(hgradBasisGradAtBasisEPoints),
333 basisCurlAtBasisCurlEPoints_(basisCurlAtBasisCurlEPoints),
334 hgradBasisGradAtTargetEPoints_(hgradBasisGradAtTargetEPoints),
335 basisCurlAtTargetCurlEPoints_(basisCurlAtTargetCurlEPoints),
336 basisEWeights_(basisEWeights), basisCurlEWeights_(basisCurlEWeights),
337 wHgradBasisGradAtBasisEPoints_(wHgradBasisGradAtBasisEPoints),
338 wBasisCurlAtBasisCurlEPoints_(wBasisCurlAtBasisCurlEPoints),
339 targetEWeights_(targetEWeights), targetCurlEWeights_(targetCurlEWeights),
340 wHgradBasisGradAtTargetEPoints_(wHgradBasisGradAtTargetEPoints),
341 wBasisCurlAtTargetCurlEPoints_(wBasisCurlAtTargetCurlEPoints),
342 computedDofs_(computedDofs), tagToOrdinal_(tagToOrdinal), hGradTagToOrdinal_(hGradTagToOrdinal),
343 numCellDofs_(numCellDofs), hgradCardinality_(hgradCardinality),
344 offsetBasis_(offsetBasis), offsetBasisCurl_(offsetBasisCurl), offsetTargetCurl_(offsetTargetCurl),
345 numEdgeFaceDofs_(numEdgeFaceDofs), dim_(dim), derDim_(derDim) {}
346
347 void
348 KOKKOS_INLINE_FUNCTION
349 operator()(const ordinal_type ic) const {
350
351 ordinal_type numBasisPoints = basisEWeights_.extent(0);
352 ordinal_type numBasisCurlPoints = basisCurlEWeights_.extent(0);
353 ordinal_type numTargetPoints = targetEWeights_.extent(0);
354 ordinal_type numTargetCurlPoints = targetCurlEWeights_.extent(0);
355 for(ordinal_type j=0; j <hgradCardinality_; ++j) {
356 ordinal_type idof = hGradTagToOrdinal_(dim_, 0, j);
357 for(ordinal_type d=0; d <dim_; ++d) {
358 for(ordinal_type iq=0; iq <numBasisPoints; ++iq)
359 wHgradBasisGradAtBasisEPoints_(ic,j,iq,d) = hgradBasisGradAtBasisEPoints_(idof,iq,d)*basisEWeights_(iq);
360 for(ordinal_type iq=0; iq <numTargetPoints; ++iq)
361 wHgradBasisGradAtTargetEPoints_(ic,j,iq,d) = hgradBasisGradAtTargetEPoints_(idof,iq,d)*targetEWeights_(iq);
362 }
363 }
364 for(ordinal_type j=0; j <numCellDofs_; ++j) {
365 ordinal_type idof = tagToOrdinal_(dim_, 0, j);
366 for(ordinal_type d=0; d <dim_; ++d)
367 for(ordinal_type iq=0; iq <numBasisPoints; ++iq)
368 cellBasisAtBasisEPoints_(ic,j,iq,d)=basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
369
370 for(ordinal_type d=0; d <derDim_; ++d) {
371 for(ordinal_type iq=0; iq <numBasisCurlPoints; ++iq) {
372 cellBasisCurlAtBasisCurlEPoints_(ic,j,iq,d)=basisCurlAtBasisCurlEPoints_(ic,idof,offsetBasisCurl_+iq,d);
373 wBasisCurlAtBasisCurlEPoints_(ic,j,iq,d)=cellBasisCurlAtBasisCurlEPoints_(ic,j,iq,d)*basisCurlEWeights_(iq);
374 }
375 for(ordinal_type iq=0; iq <numTargetCurlPoints; ++iq)
376 wBasisCurlAtTargetCurlEPoints_(ic,j,iq,d) = basisCurlAtTargetCurlEPoints_(ic,idof,offsetTargetCurl_+iq,d)*targetCurlEWeights_(iq);
377 }
378 }
379 for(ordinal_type j=0; j < numEdgeFaceDofs_; ++j) {
380 ordinal_type jdof = computedDofs_(j);
381 for(ordinal_type d=0; d <derDim_; ++d)
382 for(ordinal_type iq=0; iq <numBasisCurlPoints; ++iq)
383 negPartialProjCurl_(ic,iq,d) -= basisCoeffs_(ic,jdof)*basisCurlAtBasisCurlEPoints_(ic,jdof,offsetBasisCurl_+iq,d);
384 for(ordinal_type d=0; d <dim_; ++d)
385 for(ordinal_type iq=0; iq <numBasisPoints; ++iq)
386 negPartialProj_(ic,iq,d) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
387 }
388 }
389};
390
391
392template<typename DeviceType>
393template<typename BasisType,
394typename ortValueType, class ...ortProperties>
395void
396ProjectionTools<DeviceType>::getHCurlEvaluationPoints(typename BasisType::ScalarViewType targetEPoints,
397 typename BasisType::ScalarViewType targetCurlEPoints,
398 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
399 const BasisType* cellBasis,
401 const EvalPointsType evalPointType) {
402 const auto cellTopo = cellBasis->getBaseCellTopology();
403 ordinal_type dim = cellTopo.getDimension();
404 ordinal_type numCells = targetEPoints.extent(0);
405 const ordinal_type edgeDim = 1;
406 const ordinal_type faceDim = 2;
407
408 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
409 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
410
411 typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamEdge, subcellParamFace;
412 if(numEdges>0)
413 subcellParamEdge = RefSubcellParametrization<DeviceType>::get(edgeDim, cellTopo.getKey());
414 if(numFaces>0)
415 subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellTopo.getKey());
416
417 auto refTopologyKey = projStruct->getTopologyKey();
418
419 auto evalPointsRange = projStruct->getPointsRange(evalPointType);
420 auto curlEPointsRange = projStruct->getDerivPointsRange(evalPointType);
421
422 for(ordinal_type ie=0; ie<numEdges; ++ie) {
423
424 auto edgePointsRange = evalPointsRange(edgeDim, ie);
425 auto edgeEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(edgeDim,ie,evalPointType));
426 const auto topoKey = refTopologyKey(edgeDim, ie);
427 Kokkos::parallel_for
428 ("Evaluate Points Edges ",
429 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
430 KOKKOS_LAMBDA (const size_t ic) {
431
432 ordinal_type eOrt[12];
433 orts(ic).getEdgeOrientation(eOrt, numEdges);
434 ordinal_type ort = eOrt[ie];
435
436 Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(targetEPoints,ic,edgePointsRange,Kokkos::ALL()),
437 edgeEPoints, subcellParamEdge, topoKey, ie, ort);
438 });
439 }
440
441
442 for(ordinal_type iface=0; iface<numFaces; ++iface) {
443 auto faceCurlPointsRange = curlEPointsRange(faceDim, iface);
444 auto faceCurlEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getDerivEvalPoints(faceDim,iface,evalPointType));
445
446 auto facePointsRange = evalPointsRange(faceDim, iface);
447 auto faceEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(faceDim,iface,evalPointType));
448
449 const auto topoKey = refTopologyKey(faceDim, iface);
450 Kokkos::parallel_for
451 ("Evaluate Points Faces ",
452 Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
453 KOKKOS_LAMBDA (const size_t ic) {
454
455 ordinal_type fOrt[6];
456 orts(ic).getFaceOrientation(fOrt, numFaces);
457 ordinal_type ort = fOrt[iface];
458
459 Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(targetEPoints, ic, facePointsRange, Kokkos::ALL()),
460 faceEPoints, subcellParamFace, topoKey, iface, ort);
461
462 Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(targetCurlEPoints, ic, faceCurlPointsRange, Kokkos::ALL()),
463 faceCurlEPoints, subcellParamFace, topoKey, iface, ort);
464 });
465 }
466
467
468 if(cellBasis->getDofCount(dim,0)>0) {
469 auto cellEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,evalPointType));
470 RealSpaceTools<DeviceType>::clone(Kokkos::subview(targetEPoints, Kokkos::ALL(), evalPointsRange(dim, 0), Kokkos::ALL()), cellEPoints);
471
472 auto cellCurlEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getDerivEvalPoints(dim,0,evalPointType));
473 RealSpaceTools<DeviceType>::clone(Kokkos::subview(targetCurlEPoints, Kokkos::ALL(), curlEPointsRange(dim, 0), Kokkos::ALL()), cellCurlEPoints);
474 }
475}
476
477
478template<typename DeviceType>
479template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
480typename funValsValueType, class ...funValsProperties,
481typename BasisType,
482typename ortValueType,class ...ortProperties>
483void
484ProjectionTools<DeviceType>::getHCurlBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
485 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtTargetEPoints,
486 const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetCurlAtTargetCurlEPoints,
487 const typename BasisType::ScalarViewType targetEPoints,
488 const typename BasisType::ScalarViewType targetCurlEPoints,
489 const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
490 const BasisType* cellBasis,
492
493 typedef typename BasisType::scalarType scalarType;
494 typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
495 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
496 const auto cellTopo = cellBasis->getBaseCellTopology();
497 ordinal_type dim = cellTopo.getDimension();
498 ordinal_type numTotalTargetEPoints(targetAtTargetEPoints.extent(1)),
499 numTotalTargetCurlEPoints(targetCurlAtTargetCurlEPoints.extent(1));
500 ordinal_type basisCardinality = cellBasis->getCardinality();
501 ordinal_type numCells = targetAtTargetEPoints.extent(0);
502 const ordinal_type edgeDim = 1;
503 const ordinal_type faceDim = 2;
504 const ordinal_type derDim = dim == 3 ? dim : 1;
505
506 const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
507
508 const std::string& name = cellBasis->getName();
509
510 ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
511 ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
512
513 ScalarViewType refEdgesTangent("refEdgesTangent", numEdges, dim);
514 ScalarViewType refFacesTangents("refFaceTangents", numFaces, dim, 2);
515 ScalarViewType refFacesNormal("refFaceNormal", numFaces, dim);
516
517 ordinal_type numEdgeDofs(0);
518 for(ordinal_type ie=0; ie<numEdges; ++ie)
519 numEdgeDofs += cellBasis->getDofCount(edgeDim,ie);
520
521 ordinal_type numTotalFaceDofs(0);
522 for(ordinal_type iface=0; iface<numFaces; ++iface)
523 numTotalFaceDofs += cellBasis->getDofCount(faceDim,iface);
524
525 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal());
526
527 Kokkos::View<ordinal_type*, DeviceType> computedDofs("computedDofs",numEdgeDofs+numTotalFaceDofs);
528
529 auto targetEPointsRange = projStruct->getTargetPointsRange();
530 auto targetCurlEPointsRange = projStruct->getTargetDerivPointsRange();
531
532 auto basisEPointsRange = projStruct->getBasisPointsRange();
533 auto basisCurlEPointsRange = projStruct->getBasisDerivPointsRange();
534
535 auto refTopologyKey = projStruct->getTopologyKey();
536
537 ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(), numTotalBasisCurlEPoints = projStruct->getNumBasisDerivEvalPoints();
538
539 ScalarViewType basisEPoints("basisEPoints",numCells,numTotalBasisEPoints, dim);
540 ScalarViewType basisCurlEPoints("basisCurlEPoints",numCells,numTotalBasisCurlEPoints, dim);
541 getHCurlEvaluationPoints(basisEPoints, basisCurlEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
542
543 ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
544 ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, dim);
545 {
546 ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
547 ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalTargetEPoints, dim);
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
553 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
554 OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
555 }
556
557 ScalarViewType basisCurlAtBasisCurlEPoints;
558 ScalarViewType basisCurlAtTargetCurlEPoints;
559 if(numTotalBasisCurlEPoints>0) {
560 ScalarViewType nonOrientedBasisCurlAtTargetCurlEPoints, nonOrientedBasisCurlAtBasisCurlEPoints;
561 if (dim == 3) {
562 basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints, dim);
563 nonOrientedBasisCurlAtBasisCurlEPoints = ScalarViewType ("nonOrientedBasisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints, dim);
564 basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints, dim);
565 nonOrientedBasisCurlAtTargetCurlEPoints = ScalarViewType("nonOrientedBasisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints, dim);
566 } else {
567 basisCurlAtBasisCurlEPoints = ScalarViewType ("basisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints);
568 nonOrientedBasisCurlAtBasisCurlEPoints = ScalarViewType ("nonOrientedBasisCurlAtBasisCurlEPoints",numCells,basisCardinality, numTotalBasisCurlEPoints);
569 basisCurlAtTargetCurlEPoints = ScalarViewType("basisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints);
570 nonOrientedBasisCurlAtTargetCurlEPoints = ScalarViewType("nonOrientedBasisCurlAtTargetCurlEPoints",numCells,basisCardinality, numTotalTargetCurlEPoints);
571 }
572 for(ordinal_type ic=0; ic<numCells; ++ic) {
573 cellBasis->getValues(Kokkos::subview(nonOrientedBasisCurlAtBasisCurlEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisCurlEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_CURL);
574 cellBasis->getValues(Kokkos::subview(nonOrientedBasisCurlAtTargetCurlEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetCurlEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_CURL);
575 }
576 OrientationTools<DeviceType>::modifyBasisByOrientation(basisCurlAtBasisCurlEPoints, nonOrientedBasisCurlAtBasisCurlEPoints, orts, cellBasis);
577 OrientationTools<DeviceType>::modifyBasisByOrientation(basisCurlAtTargetCurlEPoints, nonOrientedBasisCurlAtTargetCurlEPoints, orts, cellBasis);
578 }
579
580 ordinal_type computedDofsCount = 0;
581 for(ordinal_type ie=0; ie<numEdges; ++ie) {
582
583 ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
584 ordinal_type numBasisEPoints = range_size(basisEPointsRange(edgeDim, ie));
585 ordinal_type numTargetEPoints = range_size(targetEPointsRange(edgeDim, ie));
586
587 {
588 auto refEdgeTan = Kokkos::subview(refEdgesTangent, ie, Kokkos::ALL());
589 CellTools<DeviceType>::getReferenceEdgeTangent(refEdgeTan, ie, cellTopo);
590 }
591
592 ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
593 ScalarViewType basisTanAtTargetEPoints("basisTanAtTargetEPoints",numCells,edgeCardinality, numTargetEPoints);
594 ScalarViewType weightedTanBasisAtBasisEPoints("weightedTanBasisAtBasisEPoints",numCells,edgeCardinality, numBasisEPoints);
595 ScalarViewType weightedTanBasisAtTargetEPoints("weightedTanBasisAtTargetEPoints",numCells,edgeCardinality, numTargetEPoints);
596 ScalarViewType targetTanAtTargetEPoints("normalTargetAtTargetEPoints",numCells, numTargetEPoints);
597
598 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(edgeDim,ie));
599 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(edgeDim,ie));
600
601 //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
602 ordinal_type offsetBasis = basisEPointsRange(edgeDim, ie).first;
603 ordinal_type offsetTarget = targetEPointsRange(edgeDim, ie).first;
604
605 typedef ComputeBasisCoeffsOnEdges_HCurl<ScalarViewType, decltype(basisEWeights), decltype(tagToOrdinal), decltype(targetAtTargetEPoints)> functorTypeEdge;
606 Kokkos::parallel_for(policy, functorTypeEdge(basisTanAtBasisEPoints,basisAtBasisEPoints,basisEWeights,
607 weightedTanBasisAtBasisEPoints, targetEWeights,
608 basisAtTargetEPoints, weightedTanBasisAtTargetEPoints, tagToOrdinal,
609 targetAtTargetEPoints, targetTanAtTargetEPoints,
610 refEdgesTangent, edgeCardinality, offsetBasis,
611 offsetTarget, edgeDim,
612 dim, ie));
613
614 ScalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality+1, edgeCardinality+1),
615 edgeRhsMat_("rhsMat_", numCells, edgeCardinality+1);
616
617 ScalarViewType eWeights_("eWeights_", numCells, 1, basisEWeights.extent(0)), targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0));
618 RealSpaceTools<DeviceType>::clone(eWeights_, basisEWeights);
619 RealSpaceTools<DeviceType>::clone(targetEWeights_, targetEWeights);
620
621 range_type range_H(0, edgeCardinality);
622 range_type range_B(edgeCardinality, edgeCardinality+1);
623 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(edgeMassMat_,Kokkos::ALL(),range_H,range_H), basisTanAtBasisEPoints, weightedTanBasisAtBasisEPoints);
624 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(edgeMassMat_,Kokkos::ALL(),range_H,range_B), basisTanAtBasisEPoints, eWeights_);
625 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(edgeRhsMat_,Kokkos::ALL(),range_H), targetTanAtTargetEPoints, weightedTanBasisAtTargetEPoints);
626 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(edgeRhsMat_,Kokkos::ALL(),range_B), targetTanAtTargetEPoints, targetEWeights_);
627
628 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
629 ScalarViewType t_("t",numCells, edgeCardinality+1);
630 WorkArrayViewType w_("w",numCells, edgeCardinality+1);
631
632 auto edgeDofs = Kokkos::subview(tagToOrdinal, edgeDim, ie, range_type(0,edgeCardinality));
633 ElemSystem edgeSystem("edgeSystem", false);
634 edgeSystem.solve(basisCoeffs, edgeMassMat_, edgeRhsMat_, t_, w_, edgeDofs, edgeCardinality, 1);
635
636 auto computedEdgeDofs = Kokkos::subview(computedDofs, range_type(computedDofsCount,computedDofsCount+edgeCardinality));
637 deep_copy(computedEdgeDofs, edgeDofs);
638 computedDofsCount += edgeCardinality;
639
640 }
641
642 typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamFace;
643 if(numFaces>0)
644 subcellParamFace = RefSubcellParametrization<DeviceType>::get(faceDim, cellBasis->getBaseCellTopology().getKey());
645
647 for(ordinal_type iface=0; iface<numFaces; ++iface) {
648
649 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
650 hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
651 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
652 hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
653 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Wedge<6> >()->key) {
654 if(iface < 3)
655 hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
656 else
657 hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
658 }
659 else {
660 std::stringstream ss;
661 ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
662 << "Method not implemented for basis " << name;
663 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
664 }
665
666 ordinal_type numTargetEPoints = range_size(targetEPointsRange(faceDim, iface));
667 ordinal_type numTargetCurlEPoints = range_size(targetCurlEPointsRange(faceDim, iface));
668 ordinal_type numBasisEPoints = range_size(basisEPointsRange(faceDim, iface));
669 ordinal_type numBasisCurlEPoints = range_size(basisCurlEPointsRange(faceDim, iface));
670
671 ordinal_type numFaceDofs = cellBasis->getDofCount(faceDim,iface);
672
673 ScalarViewType hgradBasisGradAtBasisEPoints("hgradBasisGradAtBasisEPoints",hgradBasis->getCardinality(), numBasisEPoints, faceDim);
674 ScalarViewType hgradBasisGradAtTargetEPoints("hgradBasisGradAtTargetEPoints",hgradBasis->getCardinality(), numTargetEPoints, faceDim);
675
676 ordinal_type hgradCardinality = hgradBasis->getDofCount(faceDim,0);
677
678 auto refBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalPoints(faceDim, iface));
679 auto refTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalPoints(faceDim, iface));
680 hgradBasis->getValues(hgradBasisGradAtBasisEPoints, refBasisEPoints, OPERATOR_GRAD);
681 hgradBasis->getValues(hgradBasisGradAtTargetEPoints, refTargetEPoints, OPERATOR_GRAD);
682
683 ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",numCells,numFaceDofs, numBasisEPoints,dim-1);
684 ScalarViewType basisTanAtTargetEPoints("basisTanAtTargetEPoints",numCells,numFaceDofs, numTargetEPoints,dim-1);
685 ScalarViewType basisCurlNormalAtBasisCurlEPoints("normaBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints);
686 ScalarViewType wNormalBasisCurlAtBasisCurlEPoints("weightedNormalBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints);
687
688 ScalarViewType targetTanAtTargetEPoints("targetTanAtTargetEPoints",numCells, numTargetEPoints, dim-1);
689 ScalarViewType normalTargetCurlAtTargetEPoints("normalTargetCurlAtTargetEPoints",numCells, numTargetCurlEPoints);
690 ScalarViewType wNormalBasisCurlBasisAtTargetCurlEPoints("weightedNormalBasisCurlAtTargetCurlEPoints",numCells,numFaceDofs, numTargetCurlEPoints);
691
692 ScalarViewType wHgradBasisGradAtBasisEPoints("wHgradBasisGradAtBasisEPoints",numCells, hgradCardinality, numBasisEPoints, faceDim);
693 ScalarViewType wHgradBasisGradAtTargetEPoints("wHgradBasisGradAtTargetEPoints",numCells, hgradCardinality, numTargetEPoints, faceDim);
694
695 ScalarViewType negPartialProjCurlNormal("mNormalComputedProjection", numCells,numBasisEPoints);
696 ScalarViewType negPartialProjTan("negPartialProjTan", numCells,numBasisEPoints,dim-1);
697
698
699 ordinal_type offsetBasis = basisEPointsRange(faceDim, iface).first;
700 ordinal_type offsetBasisCurl = basisCurlEPointsRange(faceDim, iface).first;
701 ordinal_type offsetTarget = targetEPointsRange(faceDim, iface).first;
702 ordinal_type offsetTargetCurl = targetCurlEPointsRange(faceDim, iface).first;
703
704 auto hGradTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hgradBasis->getAllDofOrdinal());
705
706 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(faceDim,iface));
707 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(faceDim,iface));
708 auto targetCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(faceDim,iface));
709 auto basisCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(faceDim,iface));
710 const auto topoKey = refTopologyKey(faceDim, iface);
711 typedef ComputeBasisCoeffsOnFaces_HCurl<decltype(basisCoeffs), decltype(orts), ScalarViewType, decltype(targetAtTargetEPoints), decltype(basisEWeights),
712 decltype(tagToOrdinal), decltype(subcellParamFace), decltype(computedDofs)> functorTypeFaces;
713 Kokkos::parallel_for(policy, functorTypeFaces(basisCoeffs,
714 orts, negPartialProjTan, negPartialProjCurlNormal,
715 hgradBasisGradAtBasisEPoints, wHgradBasisGradAtBasisEPoints,
716 basisCurlAtBasisCurlEPoints, basisCurlNormalAtBasisCurlEPoints,
717 basisAtBasisEPoints,
718 normalTargetCurlAtTargetEPoints, basisTanAtBasisEPoints,
719 hgradBasisGradAtTargetEPoints, wHgradBasisGradAtTargetEPoints,
720 wNormalBasisCurlAtBasisCurlEPoints, basisCurlAtTargetCurlEPoints,
721 wNormalBasisCurlBasisAtTargetCurlEPoints, targetAtTargetEPoints,
722 targetTanAtTargetEPoints, targetCurlAtTargetCurlEPoints,
723 basisEWeights, targetEWeights,
724 basisCurlEWeights, targetCurlEWeights, tagToOrdinal,
725 hGradTagToOrdinal, subcellParamFace,
726 computedDofs, topoKey, offsetBasis,
727 offsetBasisCurl, offsetTarget,
728 offsetTargetCurl, iface,
729 hgradCardinality, numFaces,
730 numFaceDofs, numEdgeDofs,
731 faceDim, dim));
732
733
734 ScalarViewType faceMassMat_("faceMassMat_", numCells, numFaceDofs+hgradCardinality, numFaceDofs+hgradCardinality),
735 faceRhsMat_("rhsMat_", numCells, numFaceDofs+hgradCardinality);
736 range_type range_H(0, numFaceDofs);
737 range_type range_B(numFaceDofs, numFaceDofs+hgradCardinality);
738 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_H), basisCurlNormalAtBasisCurlEPoints, wNormalBasisCurlAtBasisCurlEPoints);
739 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_B), basisTanAtBasisEPoints, wHgradBasisGradAtBasisEPoints);
740
741 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), normalTargetCurlAtTargetEPoints, wNormalBasisCurlBasisAtTargetCurlEPoints);
742 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), negPartialProjCurlNormal, wNormalBasisCurlAtBasisCurlEPoints,true);
743
744 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), targetTanAtTargetEPoints, wHgradBasisGradAtTargetEPoints);
745 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), negPartialProjTan, wHgradBasisGradAtBasisEPoints,true);
746
747
748 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
749 ScalarViewType t_("t",numCells, numFaceDofs+hgradCardinality);
750 WorkArrayViewType w_("w",numCells, numFaceDofs+hgradCardinality);
751
752 auto faceDofs = Kokkos::subview(tagToOrdinal, faceDim, iface, range_type(0,numFaceDofs));
753 ElemSystem faceSystem( "faceSystem", false);
754 faceSystem.solve(basisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDofs, numFaceDofs, hgradCardinality);
755
756 auto computedFaceDofs = Kokkos::subview(computedDofs, range_type(computedDofsCount,computedDofsCount+numFaceDofs));
757 deep_copy(computedFaceDofs, faceDofs);
758 computedDofsCount += numFaceDofs;
759
760 delete hgradBasis;
761 }
762
763 ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
764 if(numCellDofs>0) {
765 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
766 hgradBasis = new Basis_HGRAD_HEX_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
767 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
768 hgradBasis = new Basis_HGRAD_TET_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
769 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Wedge<6> >()->key)
770 hgradBasis = new typename DerivedNodalBasisFamily<DeviceType,scalarType,scalarType>::HGRAD_WEDGE(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
771 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
772 hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
773 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
774 hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
775 else {
776 std::stringstream ss;
777 ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
778 << "Method not implemented for basis " << name;
779 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
780 }
781
782 range_type cellPointsRange = targetEPointsRange(dim, 0);
783 range_type cellCurlPointsRange = targetCurlEPointsRange(dim, 0);
784
785 ordinal_type numTargetCurlEPoints = range_size(targetCurlEPointsRange(dim,0));
786 ordinal_type numBasisCurlEPoints = range_size(basisCurlEPointsRange(dim,0));
787 ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
788 ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
789
790 ScalarViewType hgradBasisGradAtBasisEPoints("hgradBasisGradAtBasisEPoints",hgradBasis->getCardinality(), numBasisEPoints, dim);
791 ScalarViewType hgradBasisGradAtTargetEPoints("hgradBasisGradAtTargetEPoints",hgradBasis->getCardinality(), numTargetEPoints, dim);
792
793 ordinal_type hgradCardinality = hgradBasis->getDofCount(dim,0);
794 ScalarViewType wHgradBasisGradAtTargetEPoints("wHgradBasisGradAtTargetEPoints",numCells, hgradCardinality, numTargetEPoints, dim);
795 ScalarViewType wHgradBasisGradAtBasisEPoints("wHgradBasisGradAtBasisEPoints",numCells, hgradCardinality, numBasisEPoints, dim);
796
797 hgradBasis->getValues(hgradBasisGradAtBasisEPoints,Kokkos::subview(basisEPoints, 0, basisEPointsRange(dim, 0), Kokkos::ALL()), OPERATOR_GRAD);
798 hgradBasis->getValues(hgradBasisGradAtTargetEPoints,Kokkos::subview(targetEPoints, 0, targetEPointsRange(dim, 0), Kokkos::ALL()),OPERATOR_GRAD);
799
800 ScalarViewType cellBasisAtBasisEPoints("basisCellAtEPoints",numCells,numCellDofs, numBasisEPoints, dim);
801 ScalarViewType cellBasisCurlAtCurlEPoints("cellBasisCurlAtCurlEPoints",numCells,numCellDofs, numBasisCurlEPoints, derDim);
802 ScalarViewType negPartialProjCurl("negPartialProjCurl", numCells, numBasisEPoints, derDim);
803 ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints, dim);
804 ScalarViewType wBasisCurlAtCurlEPoints("weightedBasisCurlAtBasisEPoints",numCells,numCellDofs, numBasisCurlEPoints,derDim);
805 ScalarViewType wBasisCurlBasisAtTargetCurlEPoints("weightedBasisCurlAtTargetCurlEPoints",numCells,numCellDofs, numTargetCurlEPoints,derDim);
806
807 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
808 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
809 auto targetCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0));
810 auto basisCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0));
811 ordinal_type offsetBasis = basisEPointsRange(dim, 0).first;
812 ordinal_type offsetBasisCurl = basisCurlEPointsRange(dim, 0).first;
813 ordinal_type offsetTargetCurl = targetCurlEPointsRange(dim, 0).first;
814
815
816 auto hGradTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hgradBasis->getAllDofOrdinal());
817
818 typedef ComputeBasisCoeffsOnCell_HCurl<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
819 decltype(computedDofs), decltype(tagToOrdinal)> functorTypeCell;
820 Kokkos::parallel_for(policy, functorTypeCell(basisCoeffs, negPartialProj, negPartialProjCurl,
821 cellBasisAtBasisEPoints, cellBasisCurlAtCurlEPoints,
822 basisAtBasisEPoints, hgradBasisGradAtBasisEPoints, basisCurlAtBasisCurlEPoints,
823 hgradBasisGradAtTargetEPoints, basisCurlAtTargetCurlEPoints,
824 basisEWeights, basisCurlEWeights, wHgradBasisGradAtBasisEPoints,
825 wBasisCurlAtCurlEPoints, targetEWeights, targetCurlEWeights,
826 wHgradBasisGradAtTargetEPoints,
827 wBasisCurlBasisAtTargetCurlEPoints, computedDofs,
828 tagToOrdinal, hGradTagToOrdinal,
829 numCellDofs, hgradCardinality,
830 offsetBasis, offsetBasisCurl, offsetTargetCurl,
831 numEdgeDofs+numTotalFaceDofs, dim, derDim));
832
833 ScalarViewType cellMassMat_("cellMassMat_", numCells, numCellDofs+hgradCardinality, numCellDofs+hgradCardinality),
834 cellRhsMat_("rhsMat_", numCells, numCellDofs+hgradCardinality);
835
836 range_type range_H(0, numCellDofs);
837 range_type range_B(numCellDofs, numCellDofs+hgradCardinality);
838 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_H), cellBasisCurlAtCurlEPoints, wBasisCurlAtCurlEPoints);
839 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_B), cellBasisAtBasisEPoints, wHgradBasisGradAtBasisEPoints);
840 if(dim==3)
841 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtTargetCurlEPoints,Kokkos::ALL(),cellCurlPointsRange,Kokkos::ALL()), wBasisCurlBasisAtTargetCurlEPoints);
842 else
843 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtTargetCurlEPoints,Kokkos::ALL(),cellCurlPointsRange), Kokkos::subview(wBasisCurlBasisAtTargetCurlEPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
844 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), negPartialProjCurl, wBasisCurlAtCurlEPoints, true);
845
846 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), wHgradBasisGradAtTargetEPoints);
847 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), negPartialProj, wHgradBasisGradAtBasisEPoints, true);
848
849 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
850 ScalarViewType t_("t",numCells, numCellDofs+hgradCardinality);
851 WorkArrayViewType w_("w",numCells, numCellDofs+hgradCardinality);
852
853 auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
854 ElemSystem cellSystem( "cellSystem", true);
855 cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numCellDofs, hgradCardinality);
856
857 delete hgradBasis;
858 }
859}
860}
861}
862
863#endif
864
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.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
Implementation of the default H(grad)-compatible FEM basis of degree n on Quadrilateral cell Implemen...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Tetrahedron ce...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
ordinal_type getCardinality() const
Returns cardinality of the basis.
const OrdinalTypeArray3DHost getAllDofOrdinal() const
DoF tag to ordinal data structure.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation https://trili...
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.
An helper class to compute the evaluation points and weights needed for performing projections.
ordinal_type getNumBasisEvalPoints()
Returns number of basis evaluation points.
view_type getDerivEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the evaluation points for basis/target derivatives on a subcell.
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
const range_tag getBasisPointsRange() const
Returns the range tag of the basis evaluation points subcells.
const range_tag getBasisDerivPointsRange() const
Returns the range tag of the derivative evaluation points on subcell.
const range_tag getTargetPointsRange() const
Returns the range of the target function evaluation points on subcells.
view_type getBasisEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis evaluation points on a subcell.
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 getTargetDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function derivatives evaluation weights on a subcell.
view_type getBasisEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis evaluation weights on a subcell.
const range_tag getTargetDerivPointsRange() const
Returns the range tag of the target function derivative evaluation points on subcells.
ordinal_type getNumBasisDerivEvalPoints()
Returns number of evaluation points for basis derivatives.
view_type getTargetEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function evaluation weights on a subcell.
view_type getBasisDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis derivatives evaluation weights on a subcell.
view_type getTargetEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the points where to evaluate the target function 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 getHCurlBasisCoeffs(Kokkos::DynRankView< basisCoeffsValueType, basisCoeffsProperties... > basisCoeffs, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetAtEvalPoints, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetCurlAtCurlEvalPoints, const typename BasisType::ScalarViewType evaluationPoints, const typename BasisType::ScalarViewType curlEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct)
Computes the basis coefficients of the HCurl projection of the target function.
static void getHCurlEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints, typename BasisType::ScalarViewType curlEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=EvalPointsType::TARGET)
Computes evaluation points for HCurl projection.
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...