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
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>
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>
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>
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 {
654 std::stringstream ss;
655 ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
656 << "Method not implemented for basis " << name;
657 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
658 }
659
660 ordinal_type numTargetEPoints = range_size(targetEPointsRange(faceDim, iface));
661 ordinal_type numTargetCurlEPoints = range_size(targetCurlEPointsRange(faceDim, iface));
662 ordinal_type numBasisEPoints = range_size(basisEPointsRange(faceDim, iface));
663 ordinal_type numBasisCurlEPoints = range_size(basisCurlEPointsRange(faceDim, iface));
664
665 ordinal_type numFaceDofs = cellBasis->getDofCount(faceDim,iface);
666
667 ScalarViewType hgradBasisGradAtBasisEPoints("hgradBasisGradAtBasisEPoints",hgradBasis->getCardinality(), numBasisEPoints, faceDim);
668 ScalarViewType hgradBasisGradAtTargetEPoints("hgradBasisGradAtTargetEPoints",hgradBasis->getCardinality(), numTargetEPoints, faceDim);
669
670 ordinal_type hgradCardinality = hgradBasis->getDofCount(faceDim,0);
671
672 auto refBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalPoints(faceDim, iface));
673 auto refTargetEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalPoints(faceDim, iface));
674 hgradBasis->getValues(hgradBasisGradAtBasisEPoints, refBasisEPoints, OPERATOR_GRAD);
675 hgradBasis->getValues(hgradBasisGradAtTargetEPoints, refTargetEPoints, OPERATOR_GRAD);
676
677 ScalarViewType basisTanAtBasisEPoints("basisTanAtBasisEPoints",numCells,numFaceDofs, numBasisEPoints,dim-1);
678 ScalarViewType basisTanAtTargetEPoints("basisTanAtTargetEPoints",numCells,numFaceDofs, numTargetEPoints,dim-1);
679 ScalarViewType basisCurlNormalAtBasisCurlEPoints("normaBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints);
680 ScalarViewType wNormalBasisCurlAtBasisCurlEPoints("weightedNormalBasisCurlAtBasisEPoints",numCells,numFaceDofs, numBasisCurlEPoints);
681
682 ScalarViewType targetTanAtTargetEPoints("targetTanAtTargetEPoints",numCells, numTargetEPoints, dim-1);
683 ScalarViewType normalTargetCurlAtTargetEPoints("normalTargetCurlAtTargetEPoints",numCells, numTargetCurlEPoints);
684 ScalarViewType wNormalBasisCurlBasisAtTargetCurlEPoints("weightedNormalBasisCurlAtTargetCurlEPoints",numCells,numFaceDofs, numTargetCurlEPoints);
685
686 ScalarViewType wHgradBasisGradAtBasisEPoints("wHgradBasisGradAtBasisEPoints",numCells, hgradCardinality, numBasisEPoints, faceDim);
687 ScalarViewType wHgradBasisGradAtTargetEPoints("wHgradBasisGradAtTargetEPoints",numCells, hgradCardinality, numTargetEPoints, faceDim);
688
689 ScalarViewType negPartialProjCurlNormal("mNormalComputedProjection", numCells,numBasisEPoints);
690 ScalarViewType negPartialProjTan("negPartialProjTan", numCells,numBasisEPoints,dim-1);
691
692
693 ordinal_type offsetBasis = basisEPointsRange(faceDim, iface).first;
694 ordinal_type offsetBasisCurl = basisCurlEPointsRange(faceDim, iface).first;
695 ordinal_type offsetTarget = targetEPointsRange(faceDim, iface).first;
696 ordinal_type offsetTargetCurl = targetCurlEPointsRange(faceDim, iface).first;
697
698 auto hGradTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hgradBasis->getAllDofOrdinal());
699
700 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(faceDim,iface));
701 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(faceDim,iface));
702 auto targetCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(faceDim,iface));
703 auto basisCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(faceDim,iface));
704 const auto topoKey = refTopologyKey(faceDim, iface);
705 typedef ComputeBasisCoeffsOnFaces_HCurl<decltype(basisCoeffs), decltype(orts), ScalarViewType, decltype(targetAtTargetEPoints), decltype(basisEWeights),
706 decltype(tagToOrdinal), decltype(subcellParamFace), decltype(computedDofs)> functorTypeFaces;
707 Kokkos::parallel_for(policy, functorTypeFaces(basisCoeffs,
708 orts, negPartialProjTan, negPartialProjCurlNormal,
709 hgradBasisGradAtBasisEPoints, wHgradBasisGradAtBasisEPoints,
710 basisCurlAtBasisCurlEPoints, basisCurlNormalAtBasisCurlEPoints,
711 basisAtBasisEPoints,
712 normalTargetCurlAtTargetEPoints, basisTanAtBasisEPoints,
713 hgradBasisGradAtTargetEPoints, wHgradBasisGradAtTargetEPoints,
714 wNormalBasisCurlAtBasisCurlEPoints, basisCurlAtTargetCurlEPoints,
715 wNormalBasisCurlBasisAtTargetCurlEPoints, targetAtTargetEPoints,
716 targetTanAtTargetEPoints, targetCurlAtTargetCurlEPoints,
717 basisEWeights, targetEWeights,
718 basisCurlEWeights, targetCurlEWeights, tagToOrdinal,
719 hGradTagToOrdinal, subcellParamFace,
720 computedDofs, topoKey, offsetBasis,
721 offsetBasisCurl, offsetTarget,
722 offsetTargetCurl, iface,
723 hgradCardinality, numFaces,
724 numFaceDofs, numEdgeDofs,
725 faceDim, dim));
726
727
728 ScalarViewType faceMassMat_("faceMassMat_", numCells, numFaceDofs+hgradCardinality, numFaceDofs+hgradCardinality),
729 faceRhsMat_("rhsMat_", numCells, numFaceDofs+hgradCardinality);
730 range_type range_H(0, numFaceDofs);
731 range_type range_B(numFaceDofs, numFaceDofs+hgradCardinality);
732 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_H), basisCurlNormalAtBasisCurlEPoints, wNormalBasisCurlAtBasisCurlEPoints);
733 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_B), basisTanAtBasisEPoints, wHgradBasisGradAtBasisEPoints);
734
735 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), normalTargetCurlAtTargetEPoints, wNormalBasisCurlBasisAtTargetCurlEPoints);
736 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), negPartialProjCurlNormal, wNormalBasisCurlAtBasisCurlEPoints,true);
737
738 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), targetTanAtTargetEPoints, wHgradBasisGradAtTargetEPoints);
739 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), negPartialProjTan, wHgradBasisGradAtBasisEPoints,true);
740
741
742 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
743 ScalarViewType t_("t",numCells, numFaceDofs+hgradCardinality);
744 WorkArrayViewType w_("w",numCells, numFaceDofs+hgradCardinality);
745
746 auto faceDofs = Kokkos::subview(tagToOrdinal, faceDim, iface, range_type(0,numFaceDofs));
747 ElemSystem faceSystem( "faceSystem", false);
748 faceSystem.solve(basisCoeffs, faceMassMat_, faceRhsMat_, t_, w_, faceDofs, numFaceDofs, hgradCardinality);
749
750 auto computedFaceDofs = Kokkos::subview(computedDofs, range_type(computedDofsCount,computedDofsCount+numFaceDofs));
751 deep_copy(computedFaceDofs, faceDofs);
752 computedDofsCount += numFaceDofs;
753
754 delete hgradBasis;
755 }
756
757 ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
758 if(numCellDofs>0) {
759 if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
760 hgradBasis = new Basis_HGRAD_HEX_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
761 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
762 hgradBasis = new Basis_HGRAD_TET_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
763 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
764 hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
765 else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
766 hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
767 else {
768 std::stringstream ss;
769 ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
770 << "Method not implemented for basis " << name;
771 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
772 }
773
774 range_type cellPointsRange = targetEPointsRange(dim, 0);
775 range_type cellCurlPointsRange = targetCurlEPointsRange(dim, 0);
776
777 ordinal_type numTargetCurlEPoints = range_size(targetCurlEPointsRange(dim,0));
778 ordinal_type numBasisCurlEPoints = range_size(basisCurlEPointsRange(dim,0));
779 ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
780 ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
781
782 ScalarViewType hgradBasisGradAtBasisEPoints("hgradBasisGradAtBasisEPoints",hgradBasis->getCardinality(), numBasisEPoints, dim);
783 ScalarViewType hgradBasisGradAtTargetEPoints("hgradBasisGradAtTargetEPoints",hgradBasis->getCardinality(), numTargetEPoints, dim);
784
785 ordinal_type hgradCardinality = hgradBasis->getDofCount(dim,0);
786 ScalarViewType wHgradBasisGradAtTargetEPoints("wHgradBasisGradAtTargetEPoints",numCells, hgradCardinality, numTargetEPoints, dim);
787 ScalarViewType wHgradBasisGradAtBasisEPoints("wHgradBasisGradAtBasisEPoints",numCells, hgradCardinality, numBasisEPoints, dim);
788
789 hgradBasis->getValues(hgradBasisGradAtBasisEPoints,Kokkos::subview(basisEPoints, 0, basisEPointsRange(dim, 0), Kokkos::ALL()), OPERATOR_GRAD);
790 hgradBasis->getValues(hgradBasisGradAtTargetEPoints,Kokkos::subview(targetEPoints, 0, targetEPointsRange(dim, 0), Kokkos::ALL()),OPERATOR_GRAD);
791
792 ScalarViewType cellBasisAtBasisEPoints("basisCellAtEPoints",numCells,numCellDofs, numBasisEPoints, dim);
793 ScalarViewType cellBasisCurlAtCurlEPoints("cellBasisCurlAtCurlEPoints",numCells,numCellDofs, numBasisCurlEPoints, derDim);
794 ScalarViewType negPartialProjCurl("negPartialProjCurl", numCells, numBasisEPoints, derDim);
795 ScalarViewType negPartialProj("negPartialProj", numCells, numBasisEPoints, dim);
796 ScalarViewType wBasisCurlAtCurlEPoints("weightedBasisCurlAtBasisEPoints",numCells,numCellDofs, numBasisCurlEPoints,derDim);
797 ScalarViewType wBasisCurlBasisAtTargetCurlEPoints("weightedBasisCurlAtTargetCurlEPoints",numCells,numCellDofs, numTargetCurlEPoints,derDim);
798
799 auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
800 auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
801 auto targetCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0));
802 auto basisCurlEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0));
803 ordinal_type offsetBasis = basisEPointsRange(dim, 0).first;
804 ordinal_type offsetBasisCurl = basisCurlEPointsRange(dim, 0).first;
805 ordinal_type offsetTargetCurl = targetCurlEPointsRange(dim, 0).first;
806
807
808 auto hGradTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hgradBasis->getAllDofOrdinal());
809
810 typedef ComputeBasisCoeffsOnCell_HCurl<decltype(basisCoeffs), ScalarViewType, decltype(basisEWeights),
811 decltype(computedDofs), decltype(tagToOrdinal)> functorTypeCell;
812 Kokkos::parallel_for(policy, functorTypeCell(basisCoeffs, negPartialProj, negPartialProjCurl,
813 cellBasisAtBasisEPoints, cellBasisCurlAtCurlEPoints,
814 basisAtBasisEPoints, hgradBasisGradAtBasisEPoints, basisCurlAtBasisCurlEPoints,
815 hgradBasisGradAtTargetEPoints, basisCurlAtTargetCurlEPoints,
816 basisEWeights, basisCurlEWeights, wHgradBasisGradAtBasisEPoints,
817 wBasisCurlAtCurlEPoints, targetEWeights, targetCurlEWeights,
818 wHgradBasisGradAtTargetEPoints,
819 wBasisCurlBasisAtTargetCurlEPoints, computedDofs,
820 tagToOrdinal, hGradTagToOrdinal,
821 numCellDofs, hgradCardinality,
822 offsetBasis, offsetBasisCurl, offsetTargetCurl,
823 numEdgeDofs+numTotalFaceDofs, dim, derDim));
824
825 ScalarViewType cellMassMat_("cellMassMat_", numCells, numCellDofs+hgradCardinality, numCellDofs+hgradCardinality),
826 cellRhsMat_("rhsMat_", numCells, numCellDofs+hgradCardinality);
827
828 range_type range_H(0, numCellDofs);
829 range_type range_B(numCellDofs, numCellDofs+hgradCardinality);
830 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_H), cellBasisCurlAtCurlEPoints, wBasisCurlAtCurlEPoints);
831 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_B), cellBasisAtBasisEPoints, wHgradBasisGradAtBasisEPoints);
832 if(dim==3)
833 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtTargetCurlEPoints,Kokkos::ALL(),cellCurlPointsRange,Kokkos::ALL()), wBasisCurlBasisAtTargetCurlEPoints);
834 else
835 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));
836 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), negPartialProjCurl, wBasisCurlAtCurlEPoints, true);
837
838 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), Kokkos::subview(targetAtTargetEPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), wHgradBasisGradAtTargetEPoints);
839 FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), negPartialProj, wHgradBasisGradAtBasisEPoints, true);
840
841 typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
842 ScalarViewType t_("t",numCells, numCellDofs+hgradCardinality);
843 WorkArrayViewType w_("w",numCells, numCellDofs+hgradCardinality);
844
845 auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
846 ElemSystem cellSystem( "cellSystem", true);
847 cellSystem.solve(basisCoeffs, cellMassMat_, cellRhsMat_, t_, w_, cellDofs, numCellDofs, hgradCardinality);
848
849 delete hgradBasis;
850 }
851}
852}
853}
854
855#endif
856
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.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
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.
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...