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 
53 #include "Intrepid2_ArrayTools.hpp"
55 
56 
57 namespace Intrepid2 {
58 namespace Experimental {
59 
60 template<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 
124 template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4,
125 typename 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 
284 template<typename ViewType1, typename ViewType2, typename ViewType3,
285 typename 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 
392 template<typename DeviceType>
393 template<typename BasisType,
394 typename ortValueType, class ...ortProperties>
395 void
396 ProjectionTools<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 
478 template<typename DeviceType>
479 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
480 typename funValsValueType, class ...funValsProperties,
481 typename BasisType,
482 typename ortValueType,class ...ortProperties>
483 void
484 ProjectionTools<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 
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 
646  Basis<DeviceType,scalarType,scalarType> *hgradBasis = NULL;
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 the abstract base class Intrepid2::DefaultCubatureFactory.
ordinal_type getNumBasisDerivEvalPoints()
Returns number of evaluation points for basis derivatives.
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 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.
const range_tag getTargetDerivPointsRange() const
Returns the range tag of the target function derivative evaluation points on subcells.
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties... > refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
const range_tag getBasisPointsRange() const
Returns the range tag of the basis evaluation points subcells.
view_type getTargetEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the points where to evaluate the target function on a subcell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
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.
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.
const range_tag getBasisDerivPointsRange() const
Returns the range tag of the derivative evaluation points on subcell.
view_type getTargetDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function derivatives evaluation weights on a subcell.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
view_type getBasisEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis 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.
Header file for the Intrepid2::FunctionSpaceTools class.
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Tetrahedron ce...
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
view_type getBasisEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis evaluation points on a subcell.
Class to solve a square system A x = b on each cell A is expected to be saddle a point (KKT) matrix o...
ordinal_type getNumBasisEvalPoints()
Returns number of basis evaluation points.
Header file for Intrepid2::ArrayTools class providing utilities for array operations.
static void integrate(Kokkos::DynRankView< outputValueValueType, outputValueProperties... > outputValues, const Kokkos::DynRankView< leftValueValueType, leftValueProperties... > leftValues, const Kokkos::DynRankView< rightValueValueType, rightValueProperties... > rightValues, const bool sumInto=false)
Contracts leftValues and rightValues arrays on the point and possibly space dimensions and stores the...
const range_tag getTargetPointsRange() const
Returns the range of the target function evaluation points on subcells.
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell...
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.
void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau, ViewType3 w, const ViewType4 elemDof, ordinal_type n, ordinal_type m=0)
Solve the system and returns the basis coefficients solve the system either using Kokkos Kernel QR or...
An helper class to compute the evaluation points and weights needed for performing projections...
static KOKKOS_INLINE_FUNCTION void getRefSideTangentsAndNormal(TanNormViewType tangentsAndNormal, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) tangents and normal of the side subCell The normals are only defined for side...
Implementation of the default H(grad)-compatible FEM basis of degree n on Quadrilateral cell Implemen...
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
const key_tag getTopologyKey() const
Returns the key tag for subcells.
view_type getEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the basis/target evaluation points on a subcell.
view_type getTargetEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function evaluation weights on a subcell.