Intrepid2
Intrepid2_ProjectionToolsDefHDIV.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_PROJECTIONTOOLSDEFHDIV_HPP__
50 #define __INTREPID2_PROJECTIONTOOLSDEFHDIV_HPP__
51 
53 #include "Intrepid2_ArrayTools.hpp"
55 
56 
57 namespace Intrepid2 {
58 namespace Experimental {
59 
60 
61 template<typename ViewType1, typename ViewType2, typename ViewType3, typename ViewType4>
63  const ViewType1 sideBasisNormalAtBasisEPoints_;
64  const ViewType1 basisAtBasisEPoints_;
65  const ViewType2 basisEWeights_;
66  const ViewType1 wBasisDofAtBasisEPoints_;
67  const ViewType2 targetEWeights_;
68  const ViewType1 basisAtTargetEPoints_;
69  const ViewType1 wBasisDofAtTargetEPoints_;
70  const ViewType3 tagToOrdinal_;
71  const ViewType4 targetAtEPoints_;
72  const ViewType1 targetAtTargetEPoints_;
73  const ViewType1 refSidesNormal_;
74  ordinal_type sideCardinality_;
75  ordinal_type offsetBasis_;
76  ordinal_type offsetTarget_;
77  ordinal_type sideDim_;
78  ordinal_type dim_;
79  ordinal_type iside_;
80 
81  ComputeBasisCoeffsOnSides_HDiv(const ViewType1 sideBasisNormalAtBasisEPoints,
82  const ViewType1 basisAtBasisEPoints, const ViewType2 basisEWeights, const ViewType1 wBasisDofAtBasisEPoints, const ViewType2 targetEWeights,
83  const ViewType1 basisAtTargetEPoints, const ViewType1 wBasisDofAtTargetEvalPoint, const ViewType3 tagToOrdinal,
84  const ViewType4 targetAtEPoints, const ViewType1 targetAtTargetEPoints,
85  const ViewType1 refSidesNormal, ordinal_type sideCardinality, ordinal_type offsetBasis,
86  ordinal_type offsetTarget, ordinal_type sideDim,
87  ordinal_type dim, ordinal_type iside) :
88  sideBasisNormalAtBasisEPoints_(sideBasisNormalAtBasisEPoints),
89  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisDofAtBasisEPoints_(wBasisDofAtBasisEPoints), targetEWeights_(targetEWeights),
90  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisDofAtTargetEPoints_(wBasisDofAtTargetEvalPoint),
91  tagToOrdinal_(tagToOrdinal), targetAtEPoints_(targetAtEPoints),
92  targetAtTargetEPoints_(targetAtTargetEPoints),
93  refSidesNormal_(refSidesNormal), sideCardinality_(sideCardinality), offsetBasis_(offsetBasis),
94  offsetTarget_(offsetTarget), sideDim_(sideDim), dim_(dim), iside_(iside)
95  {}
96 
97  void
98  KOKKOS_INLINE_FUNCTION
99  operator()(const ordinal_type ic) const {
100 
101  //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
102  for(ordinal_type j=0; j <sideCardinality_; ++j) {
103  ordinal_type jdof = tagToOrdinal_(sideDim_, iside_, j);
104 
105  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
106  for(ordinal_type d=0; d <dim_; ++d)
107  sideBasisNormalAtBasisEPoints_(ic,j,iq) += refSidesNormal_(iside_,d)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq,d);
108  wBasisDofAtBasisEPoints_(ic,j,iq) = sideBasisNormalAtBasisEPoints_(ic,j,iq) * basisEWeights_(iq);
109  }
110  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
111  typename ViewType2::value_type sum=0;
112  for(ordinal_type d=0; d <dim_; ++d)
113  sum += refSidesNormal_(iside_,d)*basisAtTargetEPoints_(ic,jdof,offsetTarget_+iq,d);
114  wBasisDofAtTargetEPoints_(ic,j,iq) = sum * targetEWeights_(iq);
115  }
116  }
117 
118  for(ordinal_type d=0; d <dim_; ++d)
119  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq)
120  targetAtTargetEPoints_(ic,iq) += refSidesNormal_(iside_,d)*targetAtEPoints_(ic,offsetTarget_+iq,d);
121  }
122 };
123 
124 
125 template<typename ViewType1, typename ViewType2, typename ViewType3,
126 typename ViewType4, typename ViewType5>
128  const ViewType1 basisCoeffs_;
129  const ViewType2 negPartialProjAtBasisEPoints_;
130  const ViewType2 nonWeightedBasisAtBasisEPoints_;
131  const ViewType2 basisAtBasisEPoints_;
132  const ViewType3 basisEWeights_;
133  const ViewType2 wBasisAtBasisEPoints_;
134  const ViewType3 targetEWeights_;
135  const ViewType2 basisAtTargetEPoints_;
136  const ViewType2 wBasisAtTargetEPoints_;
137  const ViewType4 computedDofs_;
138  const ViewType5 cellDof_;
139  ordinal_type numCellDofs_;
140  ordinal_type offsetBasis_;
141  ordinal_type offsetTarget_;
142  ordinal_type numSideDofs_;
143 
144  ComputeBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType2 negPartialProjAtBasisEPoints, const ViewType2 nonWeightedBasisAtBasisEPoints,
145  const ViewType2 basisAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wBasisAtBasisEPoints, const ViewType3 targetEWeights,
146  const ViewType2 basisAtTargetEPoints, const ViewType2 wBasisAtTargetEPoints, const ViewType4 computedDofs, const ViewType5 cellDof,
147  ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type offsetTarget, ordinal_type numSideDofs) :
148  basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints),
149  basisAtBasisEPoints_(basisAtBasisEPoints), basisEWeights_(basisEWeights), wBasisAtBasisEPoints_(wBasisAtBasisEPoints), targetEWeights_(targetEWeights),
150  basisAtTargetEPoints_(basisAtTargetEPoints), wBasisAtTargetEPoints_(wBasisAtTargetEPoints),
151  computedDofs_(computedDofs), cellDof_(cellDof),numCellDofs_(numCellDofs), offsetBasis_(offsetBasis),
152  offsetTarget_(offsetTarget), numSideDofs_(numSideDofs) {}
153 
154  void
155  KOKKOS_INLINE_FUNCTION
156  operator()(const ordinal_type ic) const {
157 
158  for(ordinal_type j=0; j <numCellDofs_; ++j) {
159  ordinal_type idof = cellDof_(j);
160  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq) {
161  nonWeightedBasisAtBasisEPoints_(ic,j,iq) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq);
162  wBasisAtBasisEPoints_(ic,j,iq) = nonWeightedBasisAtBasisEPoints_(ic,j,iq) * basisEWeights_(iq);
163  }
164  for(ordinal_type iq=0; iq <ordinal_type(targetEWeights_.extent(0)); ++iq) {
165  wBasisAtTargetEPoints_(ic,j,iq) = basisAtTargetEPoints_(ic,idof,offsetTarget_+iq)* targetEWeights_(iq);
166  }
167  }
168  for(ordinal_type j=0; j < numSideDofs_; ++j) {
169  ordinal_type jdof = computedDofs_(j);
170  for(ordinal_type iq=0; iq <ordinal_type(basisEWeights_.extent(0)); ++iq)
171  negPartialProjAtBasisEPoints_(ic,iq) -= basisCoeffs_(ic,jdof)*basisAtBasisEPoints_(ic,jdof,offsetBasis_+iq);
172  }
173  }
174 };
175 
176 
177 template<typename ViewType1, typename ViewType2, typename ViewType3,
178 typename ViewType4, typename ViewType5, typename ViewType6>
180  const ViewType1 basisCoeffs_;
181  const ViewType2 negPartialProjAtBasisEPoints_;
182  const ViewType2 nonWeightedBasisAtBasisEPoints_;
183  const ViewType2 basisAtBasisEPoints_;
184  const ViewType2 hcurlBasisCurlAtBasisEPoints_;
185  const ViewType3 basisEWeights_;
186  const ViewType2 wHCurlBasisAtBasisEPoints_;
187  const ViewType3 targetEWeights_;
188  const ViewType2 hcurlBasisCurlAtTargetEPoints_;
189  const ViewType2 wHCurlBasisAtTargetEPoints_;
190  const ViewType4 tagToOrdinal_;
191  const ViewType5 computedDofs_;
192  const ViewType6 hCurlDof_;
193  ordinal_type numCellDofs_;
194  ordinal_type offsetBasis_;
195  ordinal_type numSideDofs_;
196  ordinal_type dim_;
197 
198  ComputeHCurlBasisCoeffsOnCells_HDiv(const ViewType1 basisCoeffs, ViewType2 negPartialProjAtBasisEPoints, const ViewType2 nonWeightedBasisAtBasisEPoints,
199  const ViewType2 basisAtBasisEPoints, const ViewType2 hcurlBasisCurlAtBasisEPoints, const ViewType3 basisEWeights, const ViewType2 wHCurlBasisAtBasisEPoints, const ViewType3 targetEWeights,
200  const ViewType2 hcurlBasisCurlAtTargetEPoints, const ViewType2 wHCurlBasisAtTargetEPoints, const ViewType4 tagToOrdinal, const ViewType5 computedDofs, const ViewType6 hCurlDof,
201  ordinal_type numCellDofs, ordinal_type offsetBasis, ordinal_type numSideDofs, ordinal_type dim) :
202  basisCoeffs_(basisCoeffs), negPartialProjAtBasisEPoints_(negPartialProjAtBasisEPoints), nonWeightedBasisAtBasisEPoints_(nonWeightedBasisAtBasisEPoints),
203  basisAtBasisEPoints_(basisAtBasisEPoints), hcurlBasisCurlAtBasisEPoints_(hcurlBasisCurlAtBasisEPoints), basisEWeights_(basisEWeights), wHCurlBasisAtBasisEPoints_(wHCurlBasisAtBasisEPoints), targetEWeights_(targetEWeights),
204  hcurlBasisCurlAtTargetEPoints_(hcurlBasisCurlAtTargetEPoints), wHCurlBasisAtTargetEPoints_(wHCurlBasisAtTargetEPoints),
205  tagToOrdinal_(tagToOrdinal), computedDofs_(computedDofs), hCurlDof_(hCurlDof),numCellDofs_(numCellDofs), offsetBasis_(offsetBasis),
206  numSideDofs_(numSideDofs), dim_(dim) {}
207 
208  void
209  KOKKOS_INLINE_FUNCTION
210  operator()(const ordinal_type ic) const {
211 
212  ordinal_type numBasisEPoints = basisEWeights_.extent(0);
213 
214  for(ordinal_type i=0; i<numSideDofs_; ++i) {
215  ordinal_type idof = computedDofs_(i);
216  for(ordinal_type iq=0; iq<numBasisEPoints; ++iq){
217  for(ordinal_type d=0; d <dim_; ++d)
218  negPartialProjAtBasisEPoints_(ic,iq,d) -= basisCoeffs_(ic,idof)*basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
219  }
220  }
221 
222  for(ordinal_type i=0; i<numCellDofs_; ++i) {
223  ordinal_type idof = tagToOrdinal_(dim_, 0, i);
224  for(ordinal_type iq=0; iq<numBasisEPoints; ++iq) {
225  for(ordinal_type d=0; d<dim_; ++d)
226  nonWeightedBasisAtBasisEPoints_(ic,i,iq,d) = basisAtBasisEPoints_(ic,idof,offsetBasis_+iq,d);
227  }
228  }
229 
230  for(ordinal_type i=0; i<static_cast<ordinal_type>(hCurlDof_.extent(0)); ++i) {
231  ordinal_type idof = hCurlDof_(i);
232  for(ordinal_type d=0; d<dim_; ++d) {
233  for(ordinal_type iq=0; iq<numBasisEPoints; ++iq) {
234  wHCurlBasisAtBasisEPoints_(ic,i,iq,d) = hcurlBasisCurlAtBasisEPoints_(idof,iq,d)*basisEWeights_(iq);
235  }
236  for(ordinal_type iq=0; iq<static_cast<ordinal_type>(targetEWeights_.extent(0)); ++iq) {
237  wHCurlBasisAtTargetEPoints_(ic,i,iq,d) = hcurlBasisCurlAtTargetEPoints_(idof,iq,d)*targetEWeights_(iq);
238  }
239  }
240  }
241  }
242 };
243 
244 
245 template<typename DeviceType>
246 template<typename BasisType,
247 typename ortValueType, class ...ortProperties>
248 void
249 ProjectionTools<DeviceType>::getHDivEvaluationPoints(typename BasisType::ScalarViewType targetEPoints,
250  typename BasisType::ScalarViewType targetDivEPoints,
251  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
252  const BasisType* cellBasis,
254  const EvalPointsType evalPointType){
255  auto refTopologyKey = projStruct->getTopologyKey();
256  //const auto cellTopoKey = cellBasis->getBaseCellTopology().getKey();
257  ordinal_type dim = cellBasis->getBaseCellTopology().getDimension();
258  ordinal_type sideDim = dim-1;
259 
260  ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
261 
262  ordinal_type numCells = orts.extent(0);
263 
264  typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParamSide;
265  if(numSides>0)
266  subcellParamSide = RefSubcellParametrization<DeviceType>::get(sideDim, cellBasis->getBaseCellTopology().getKey());
267 
268  auto evalPointsRange = projStruct->getPointsRange(evalPointType);
269 
270  for(ordinal_type is=0; is<numSides; ++is) {
271 
272  const auto topoKey = refTopologyKey(sideDim,is);
273  auto sideBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(sideDim,is,evalPointType));
274 
275  auto sidePointsRange = evalPointsRange(sideDim, is);
276  Kokkos::parallel_for
277  ("Evaluate Points Sides ",
278  Kokkos::RangePolicy<ExecSpaceType, int> (0, numCells),
279  KOKKOS_LAMBDA (const size_t ic) {
280 
281  ordinal_type sOrt[6];
282  if(dim == 3)
283  orts(ic).getFaceOrientation(sOrt, numSides);
284  else
285  orts(ic).getEdgeOrientation(sOrt, numSides);
286  ordinal_type ort = sOrt[is];
287 
288  Impl::OrientationTools::mapSubcellCoordsToRefCell(Kokkos::subview(targetEPoints,ic,sidePointsRange,Kokkos::ALL()),
289  sideBasisEPoints, subcellParamSide, topoKey, is, ort);
290  });
291  }
292 
293  if(cellBasis->getDofCount(dim,0) <= 0)
294  return;
295 
296  auto evalDivPointsRange = projStruct->getDerivPointsRange(evalPointType);
297  auto cellDivPointsRange = evalDivPointsRange(dim, 0);
298  auto cellBasisDivEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getDerivEvalPoints(dim,0,evalPointType));
299 
300  RealSpaceTools<DeviceType>::clone(Kokkos::subview(targetDivEPoints, Kokkos::ALL(), cellDivPointsRange, Kokkos::ALL()), cellBasisDivEPoints);
301 
302  auto cellPointsRange = evalPointsRange(dim, 0);
303  if(projStruct->getTargetEvalPoints(dim, 0).data() != NULL)
304  {
305  auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getEvalPoints(dim,0,evalPointType));
306  RealSpaceTools<DeviceType>::clone(Kokkos::subview(targetEPoints, Kokkos::ALL(), cellPointsRange, Kokkos::ALL()), cellBasisEPoints);
307  }
308 }
309 
310 
311 template<typename DeviceType>
312 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
313 typename funValsValueType, class ...funValsProperties,
314 typename BasisType,
315 typename ortValueType,class ...ortProperties>
316 void
317 ProjectionTools<DeviceType>::getHDivBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
318  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEPoints,
319  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEPoints,
320  const typename BasisType::ScalarViewType targetEPoints,
321  const typename BasisType::ScalarViewType targetDivEPoints,
322  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
323  const BasisType* cellBasis,
325  typedef typename BasisType::scalarType scalarType;
326  typedef Kokkos::DynRankView<scalarType,DeviceType> ScalarViewType;
327  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
328  const auto cellTopo = cellBasis->getBaseCellTopology();
329  ordinal_type dim = cellTopo.getDimension();
330  ordinal_type numTotalEvaluationPoints(targetAtEPoints.extent(1)),
331  numTotalDivEvaluationPoints(targetDivAtDivEPoints.extent(1));
332  ordinal_type basisCardinality = cellBasis->getCardinality();
333  ordinal_type numCells = targetAtEPoints.extent(0);
334  const ordinal_type sideDim = dim-1;
335 
336  const std::string& name = cellBasis->getName();
337 
338  ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
339  ScalarViewType refSideNormal("refSideNormal", dim);
340 
341  ordinal_type numSideDofs(0);
342  for(ordinal_type is=0; is<numSides; ++is)
343  numSideDofs += cellBasis->getDofCount(sideDim,is);
344 
345  Kokkos::View<ordinal_type*, DeviceType> computedDofs("computedDofs",numSideDofs);
346 
347  const Kokkos::RangePolicy<ExecSpaceType> policy(0, numCells);
348 
349  auto targetEPointsRange = projStruct->getTargetPointsRange();
350  auto basisEPointsRange = projStruct->getBasisPointsRange();
351 
352  auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), cellBasis->getAllDofOrdinal());
353 
354  ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(), numTotalBasisDivEPoints = projStruct->getNumBasisDerivEvalPoints();
355  ScalarViewType basisEPoints_("basisEPoints",numCells,numTotalBasisEPoints, dim);
356  ScalarViewType basisDivEPoints("basisDivEPoints",numCells,numTotalBasisDivEPoints, dim);
357  getHDivEvaluationPoints(basisEPoints_, basisDivEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
358 
359  ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
360  ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
361  {
362  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
363  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
364  for(ordinal_type ic=0; ic<numCells; ++ic) {
365  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
366  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints_, ic, Kokkos::ALL(), Kokkos::ALL()));
367  }
368 
369  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
370  OrientationTools<DeviceType>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
371  }
372 
373  ScalarViewType basisDivAtBasisDivEPoints;
374  ScalarViewType basisDivAtTargetDivEPoints;
375  if(numTotalDivEvaluationPoints>0) {
376  ScalarViewType nonOrientedBasisDivAtTargetDivEPoints, nonOrientedBasisDivAtBasisDivEPoints;
377  basisDivAtBasisDivEPoints = ScalarViewType ("basisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
378  nonOrientedBasisDivAtBasisDivEPoints = ScalarViewType ("nonOrientedBasisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
379  basisDivAtTargetDivEPoints = ScalarViewType("basisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
380  nonOrientedBasisDivAtTargetDivEPoints = ScalarViewType("nonOrientedBasisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
381 
382  for(ordinal_type ic=0; ic<numCells; ++ic) {
383  cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtBasisDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
384  cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtTargetDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
385  }
386  OrientationTools<DeviceType>::modifyBasisByOrientation(basisDivAtBasisDivEPoints, nonOrientedBasisDivAtBasisDivEPoints, orts, cellBasis);
387  OrientationTools<DeviceType>::modifyBasisByOrientation(basisDivAtTargetDivEPoints, nonOrientedBasisDivAtTargetDivEPoints, orts, cellBasis);
388  }
389 
390  ScalarViewType refSidesNormal("refSidesNormal", numSides, dim);
391  ordinal_type computedDofsCount = 0;
392  for(ordinal_type is=0; is<numSides; ++is) {
393 
394  ordinal_type sideCardinality = cellBasis->getDofCount(sideDim,is);
395 
396  ordinal_type numTargetEPoints = range_size(targetEPointsRange(sideDim,is));
397  ordinal_type numBasisEPoints = range_size(basisEPointsRange(sideDim,is));
398 
399  auto sideDofs = Kokkos::subview(tagToOrdinal, sideDim, is, range_type(0,sideCardinality));
400  auto computedSideDofs = Kokkos::subview(computedDofs, range_type(computedDofsCount,computedDofsCount+sideCardinality));
401  deep_copy(computedSideDofs, sideDofs);
402  computedDofsCount += sideCardinality;
403 
404  auto sideNormal = Kokkos::subview(refSidesNormal,is,Kokkos::ALL());
405  CellTools<DeviceType>::getReferenceSideNormal(sideNormal, is, cellTopo);
406 
407  ScalarViewType basisNormalAtBasisEPoints("normalBasisAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
408  ScalarViewType wBasisNormalAtBasisEPoints("wBasisNormalAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
409  ScalarViewType wBasisNormalAtTargetEPoints("wBasisNormalAtTargetEPoints",numCells,sideCardinality, numTargetEPoints);
410  ScalarViewType targetNormalAtTargetEPoints("targetNormalAtTargetEPoints",numCells, numTargetEPoints);
411 
412  ordinal_type offsetBasis = basisEPointsRange(sideDim,is).first;
413  ordinal_type offsetTarget = targetEPointsRange(sideDim,is).first;
414  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(sideDim,is));
415  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(sideDim,is));
416 
417 
419  Kokkos::parallel_for(policy, functorTypeSide(basisNormalAtBasisEPoints, basisAtBasisEPoints,
420  basisEWeights, wBasisNormalAtBasisEPoints, targetEWeights,
421  basisAtTargetEPoints, wBasisNormalAtTargetEPoints, tagToOrdinal,
422  targetAtEPoints, targetNormalAtTargetEPoints,
423  refSidesNormal, sideCardinality, offsetBasis,
424  offsetTarget, sideDim,
425  dim, is));
426 
427  ScalarViewType sideMassMat_("sideMassMat_", numCells, sideCardinality+1, sideCardinality+1),
428  sideRhsMat_("rhsMat_", numCells, sideCardinality+1);
429 
430  ScalarViewType targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0));
431  RealSpaceTools<DeviceType>::clone(targetEWeights_, targetEWeights);
432 
433  range_type range_H(0, sideCardinality);
434  range_type range_B(sideCardinality, sideCardinality+1);
435  ScalarViewType ones("ones",numCells,1,numBasisEPoints);
436  Kokkos::deep_copy(ones,1);
437 
438  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(sideMassMat_, Kokkos::ALL(), range_H, range_H), basisNormalAtBasisEPoints, wBasisNormalAtBasisEPoints);
439  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(sideMassMat_, Kokkos::ALL(), range_H, range_B), wBasisNormalAtBasisEPoints, ones);
440 
441  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(sideRhsMat_, Kokkos::ALL(), range_H), targetNormalAtTargetEPoints, wBasisNormalAtTargetEPoints);
442  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(sideRhsMat_, Kokkos::ALL(), range_B), targetNormalAtTargetEPoints, targetEWeights_);
443 
444  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
445  ScalarViewType t_("t",numCells, sideCardinality+1);
446  WorkArrayViewType w_("w",numCells, sideCardinality+1);
447 
448  auto sideDof = Kokkos::subview(tagToOrdinal, sideDim, is, Kokkos::ALL());
449 
450  ElemSystem sideSystem("sideSystem", false);
451  sideSystem.solve(basisCoeffs, sideMassMat_, sideRhsMat_, t_, w_, sideDof, sideCardinality, 1);
452  }
453 
454 
455  //Cell
456  ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
457  if(numCellDofs==0)
458  return;
459 
460  Basis<DeviceType,scalarType,scalarType> *hcurlBasis = NULL;
461  if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
462  hcurlBasis = new Basis_HCURL_HEX_In_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
463  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
464  hcurlBasis = new Basis_HCURL_TET_In_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
465  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
466  hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
467  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
468  hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM<DeviceType,scalarType,scalarType>(cellBasis->getDegree());
469  else {
470  std::stringstream ss;
471  ss << ">>> ERROR (Intrepid2::ProjectionTools::getHDivEvaluationPoints): "
472  << "Method not implemented for basis " << name;
473  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
474  }
475  if(hcurlBasis == NULL) return;
476 
477 
478  auto targetDivEPointsRange = projStruct->getTargetDerivPointsRange();
479  auto basisDivEPointsRange = projStruct->getBasisDerivPointsRange();
480 
481  ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0));
482  ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0));
483 
484  auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetDerivEvalWeights(dim,0));
485  auto divEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisDerivEvalWeights(dim,0));
486 
487  ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first;
488  ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first;
489 
490  ScalarViewType weightedBasisDivAtBasisEPoints("weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
491  ScalarViewType weightedBasisDivAtTargetEPoints("weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints);
492  ScalarViewType basisDivAtBasisEPoints("basisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
493  ScalarViewType targetSideDivAtBasisEPoints("targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints);
494 
495  auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
497  Kokkos::parallel_for(policy, functorType( basisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints,
498  basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints,
499  computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs));
500 
501 
502  ordinal_type hcurlBasisCardinality = hcurlBasis->getCardinality();
503  ordinal_type numCurlInteriorDOFs = hcurlBasis->getDofCount(dim,0);
504 
505  range_type range_H(0, numCellDofs);
506  range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs);
507 
508 
509  ScalarViewType massMat_("massMat_",numCells,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs),
510  rhsMatTrans("rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs);
511 
512  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_H), basisDivAtBasisEPoints, weightedBasisDivAtBasisEPoints);
513  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetDivAtDivEPoints, weightedBasisDivAtTargetEPoints);
514  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetSideDivAtBasisEPoints, weightedBasisDivAtBasisEPoints,true);
515 
516  if(numCurlInteriorDOFs>0){
517  ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
518  ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
519 
520  auto targetEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getTargetEvalWeights(dim,0));
521  auto basisEWeights = Kokkos::create_mirror_view_and_copy(MemSpaceType(),projStruct->getBasisEvalWeights(dim,0));
522 
523  ordinal_type offsetBasis = basisEPointsRange(dim,0).first;
524 
525  auto basisEPoints = Kokkos::subview(basisEPoints_, 0, basisEPointsRange(dim,0), Kokkos::ALL());
526 
527  ScalarViewType negPartialProjAtBasisEPoints("targetSideAtBasisEPoints",numCells, numBasisEPoints, dim);
528  ScalarViewType nonWeightedBasisAtBasisEPoints("basisAtBasisEPoints",numCells,numCellDofs, numBasisEPoints, dim);
529  ScalarViewType hcurlBasisCurlAtBasisEPoints("hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim);
530  ScalarViewType hcurlBasisCurlAtTargetEPoints("hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim);
531  ScalarViewType wHcurlBasisCurlAtBasisEPoints("wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim);
532  ScalarViewType wHcurlBasisCurlAtTargetEPoints("wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim);
533 
534  hcurlBasis->getValues(hcurlBasisCurlAtBasisEPoints, basisEPoints, OPERATOR_CURL);
535  hcurlBasis->getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,0,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL);
536 
537  auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(MemSpaceType(), hcurlBasis->getAllDofOrdinal());
538  auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs));
539 
540  typedef ComputeHCurlBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights),
541  decltype(tagToOrdinal), decltype(computedDofs), decltype(cellDofs)> functorTypeHCurlCells;
542  Kokkos::parallel_for(policy, functorTypeHCurlCells(basisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints,
543  basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights,
544  hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof,
545  numCellDofs, offsetBasis, numSideDofs, dim));
546 
547  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_B), nonWeightedBasisAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints);
548  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints);
549  FunctionSpaceTools<DeviceType >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), negPartialProjAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints,true);
550  }
551  delete hcurlBasis;
552 
553  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, DeviceType> WorkArrayViewType;
554  ScalarViewType t_("t",numCells, numCellDofs+numCurlInteriorDOFs);
555  WorkArrayViewType w_("w",numCells, numCellDofs+numCurlInteriorDOFs);
556 
557  ElemSystem cellSystem("cellSystem", true);
558  cellSystem.solve(basisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs);
559 }
560 
561 
562 
563 }
564 }
565 
566 #endif
567 
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.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
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...
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.
const range_tag getBasisDerivPointsRange() const
Returns the range tag of the derivative evaluation points on subcell.
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
view_type getTargetDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function derivatives evaluation weights on a subcell.
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
ordinal_type getCardinality() const
Returns cardinality of the basis.
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.
const OrdinalTypeArray3DHost getAllDofOrdinal() const
DoF tag to ordinal data structure.
static void getReferenceSideNormal(Kokkos::DynRankView< refSideNormalValueType, refSideNormalProperties... > refSideNormal, const ordinal_type sideOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
Header file for the Intrepid2::FunctionSpaceTools class.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
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.
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
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...
static void getHDivEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints, typename BasisType::ScalarViewType divEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=EvalPointsType::TARGET)
Computes evaluation points for HDiv projection.
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...
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...
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.
static void getHDivBasisCoeffs(Kokkos::DynRankView< basisCoeffsValueType, basisCoeffsProperties... > basisCoeffs, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetAtEvalPoints, const Kokkos::DynRankView< funValsValueType, funValsProperties... > targetDivAtDivEvalPoints, const typename BasisType::ScalarViewType evaluationPoints, const typename BasisType::ScalarViewType divEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties... > cellOrientations, const BasisType *cellBasis, ProjectionStruct< DeviceType, typename BasisType::scalarType > *projStruct)
Computes the basis coefficients of the HDiv projection of the target function.
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.