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>
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 ViewType1 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 ViewType1 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 SpT>
246 template<typename BasisType,
247 typename ortValueType, class ...ortProperties>
248 void
249 ProjectionTools<SpT>::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  typedef typename BasisType::scalarType scalarType;
256  typedef Kokkos::DynRankView<scalarType,SpT> ScalarViewType;
257  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
258  auto refTopologyKey = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTopologyKey());
259  //const auto cellTopoKey = cellBasis->getBaseCellTopology().getKey();
260  ordinal_type dim = cellBasis->getBaseCellTopology().getDimension();
261  ordinal_type sideDim = dim-1;
262 
263  ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
264 
265  ordinal_type numCells = orts.extent(0);
266 
268  typename CellTools<SpT>::subcellParamViewType subcellParamSide;
269  if(numSides>0)
270  CellTools<SpT>::getSubcellParametrization(subcellParamSide, sideDim, cellBasis->getBaseCellTopology());
271 
272  auto evalPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getPointsRange(evalPointType));
273 
274  for(ordinal_type is=0; is<numSides; ++is) {
275  ScalarViewType sideWorkview("sideWorkview", numCells, projStruct->getMaxNumEvalPoints(evalPointType), sideDim);
276 
277  const auto topoKey = refTopologyKey(sideDim,is);
278  auto sideBasisEPoints = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getEvalPoints(sideDim,is,evalPointType));
279 
280  Kokkos::parallel_for
281  ("Evaluate Points Sides ",
282  Kokkos::RangePolicy<SpT, int> (0, numCells),
283  KOKKOS_LAMBDA (const size_t ic) {
284  auto sidePointsRange = evalPointsRange(sideDim, is);
285  auto sideRefPointsRange = range_type(0, range_size(sidePointsRange));
286  auto orientedBasisEPoints = Kokkos::subview(sideWorkview, ic, sideRefPointsRange, range_type(0,sideDim));
287 
288  ordinal_type sOrt[6];
289  if(dim == 3)
290  orts(ic).getFaceOrientation(sOrt, numSides);
291  else
292  orts(ic).getEdgeOrientation(sOrt, numSides);
293  ordinal_type ort = sOrt[is];
294 
295  Impl::OrientationTools::mapToModifiedReference(orientedBasisEPoints,sideBasisEPoints,topoKey,ort);
296  CellTools<SpT>::mapToReferenceSubcell(Kokkos::subview(targetEPoints,ic,sidePointsRange,Kokkos::ALL()), orientedBasisEPoints, subcellParamSide, sideDim, is, dim);
297  });
298  }
299 
300  if(cellBasis->getDofCount(dim,0) <= 0)
301  return;
302 
303  auto evalDivPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getDerivPointsRange(evalPointType));
304  auto cellDivPointsRange = evalDivPointsRange(dim, 0);
305  auto cellBasisDivEPoints = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getDerivEvalPoints(dim,0,evalPointType));
306 
307  RealSpaceTools<SpT>::clone(Kokkos::subview(targetDivEPoints, Kokkos::ALL(), cellDivPointsRange, Kokkos::ALL()), cellBasisDivEPoints);
308 
309 
310  if(projStruct->getTargetEvalPoints(dim, 0).data() != NULL)
311  {
312  auto cellPointsRange = evalPointsRange(dim, 0);
313  auto cellBasisEPoints = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getEvalPoints(dim,0,evalPointType));
314  RealSpaceTools<SpT>::clone(Kokkos::subview(targetEPoints, Kokkos::ALL(), cellPointsRange, Kokkos::ALL()), cellBasisEPoints);
315  }
316 }
317 
318 
319 template<typename SpT>
320 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
321 typename funValsValueType, class ...funValsProperties,
322 typename BasisType,
323 typename ortValueType,class ...ortProperties>
324 void
325 ProjectionTools<SpT>::getHDivBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
326  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEPoints,
327  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetDivAtDivEPoints,
328  const typename BasisType::ScalarViewType targetEPoints,
329  const typename BasisType::ScalarViewType targetDivEPoints,
330  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
331  const BasisType* cellBasis,
332  ProjectionStruct<SpT, typename BasisType::scalarType> * projStruct){
333  typedef typename BasisType::scalarType scalarType;
334  typedef Kokkos::DynRankView<scalarType,SpT> ScalarViewType;
335  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
336  const auto cellTopo = cellBasis->getBaseCellTopology();
337  ordinal_type dim = cellTopo.getDimension();
338  ordinal_type numTotalEvaluationPoints(targetAtEPoints.extent(1)),
339  numTotalDivEvaluationPoints(targetDivAtDivEPoints.extent(1));
340  ordinal_type basisCardinality = cellBasis->getCardinality();
341  ordinal_type numCells = targetAtEPoints.extent(0);
342  const ordinal_type sideDim = dim-1;
343 
344  const std::string& name = cellBasis->getName();
345 
346  ordinal_type numSides = cellBasis->getBaseCellTopology().getSideCount();
347  ScalarViewType refSideNormal("refSideNormal", dim);
348 
349  ordinal_type numSideDofs(0);
350  for(ordinal_type is=0; is<numSides; ++is)
351  numSideDofs += cellBasis->getDofCount(sideDim,is);
352 
353  Kokkos::View<ordinal_type*, SpT> computedDofs("computedDofs",numSideDofs);
354 
355  const Kokkos::RangePolicy<SpT> policy(0, numCells);
356 
357  auto targetEPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTargetPointsRange());
358  auto basisEPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getBasisPointsRange());
359 
360  auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(), cellBasis->getAllDofOrdinal());
361 
362  ordinal_type numTotalBasisEPoints = projStruct->getNumBasisEvalPoints(), numTotalBasisDivEPoints = projStruct->getNumBasisDerivEvalPoints();
363  ScalarViewType basisEPoints_("basisEPoints",numCells,numTotalBasisEPoints, dim);
364  ScalarViewType basisDivEPoints("basisDivEPoints",numCells,numTotalBasisDivEPoints, dim);
365  getHDivEvaluationPoints(basisEPoints_, basisDivEPoints, orts, cellBasis, projStruct, EvalPointsType::BASIS);
366 
367  ScalarViewType basisAtBasisEPoints("basisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
368  ScalarViewType basisAtTargetEPoints("basisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
369  {
370  ScalarViewType nonOrientedBasisAtBasisEPoints("nonOrientedBasisAtBasisEPoints",numCells,basisCardinality, numTotalBasisEPoints, dim);
371  ScalarViewType nonOrientedBasisAtTargetEPoints("nonOrientedBasisAtTargetEPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
372  for(ordinal_type ic=0; ic<numCells; ++ic) {
373  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetEPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
374  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtBasisEPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisEPoints_, ic, Kokkos::ALL(), Kokkos::ALL()));
375  }
376 
377  OrientationTools<SpT>::modifyBasisByOrientation(basisAtBasisEPoints, nonOrientedBasisAtBasisEPoints, orts, cellBasis);
378  OrientationTools<SpT>::modifyBasisByOrientation(basisAtTargetEPoints, nonOrientedBasisAtTargetEPoints, orts, cellBasis);
379  }
380 
381  ScalarViewType basisDivAtBasisDivEPoints;
382  ScalarViewType basisDivAtTargetDivEPoints;
383  if(numTotalDivEvaluationPoints>0) {
384  ScalarViewType nonOrientedBasisDivAtTargetDivEPoints, nonOrientedBasisDivAtBasisDivEPoints;
385  basisDivAtBasisDivEPoints = ScalarViewType ("basisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
386  nonOrientedBasisDivAtBasisDivEPoints = ScalarViewType ("nonOrientedBasisDivAtBasisDivEPoints",numCells,basisCardinality, numTotalBasisDivEPoints);
387  basisDivAtTargetDivEPoints = ScalarViewType("basisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
388  nonOrientedBasisDivAtTargetDivEPoints = ScalarViewType("nonOrientedBasisDivAtTargetDivEPoints",numCells,basisCardinality, numTotalDivEvaluationPoints);
389 
390  for(ordinal_type ic=0; ic<numCells; ++ic) {
391  cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtBasisDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(basisDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
392  cellBasis->getValues(Kokkos::subview(nonOrientedBasisDivAtTargetDivEPoints,ic,Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(targetDivEPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_DIV);
393  }
394  OrientationTools<SpT>::modifyBasisByOrientation(basisDivAtBasisDivEPoints, nonOrientedBasisDivAtBasisDivEPoints, orts, cellBasis);
395  OrientationTools<SpT>::modifyBasisByOrientation(basisDivAtTargetDivEPoints, nonOrientedBasisDivAtTargetDivEPoints, orts, cellBasis);
396  }
397 
398  ScalarViewType refSidesNormal("refSidesNormal", numSides, dim);
399  ordinal_type computedDofsCount = 0;
400  for(ordinal_type is=0; is<numSides; ++is) {
401 
402  ordinal_type sideCardinality = cellBasis->getDofCount(sideDim,is);
403 
404  ordinal_type numTargetEPoints = range_size(targetEPointsRange(sideDim,is));
405  ordinal_type numBasisEPoints = range_size(basisEPointsRange(sideDim,is));
406 
407  for(ordinal_type i=0; i<sideCardinality; ++i)
408  computedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(sideDim, is, i);
409 
410  auto sideNormal = Kokkos::subview(refSidesNormal,is,Kokkos::ALL());
411  auto sideNormalHost = Kokkos::create_mirror_view(sideNormal);
412  CellTools<SpT>::getReferenceSideNormal(sideNormalHost, is, cellTopo);
413  Kokkos::deep_copy(sideNormal, sideNormalHost);
414 
415  ScalarViewType basisNormalAtBasisEPoints("normalBasisAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
416  ScalarViewType wBasisNormalAtBasisEPoints("wBasisNormalAtBasisEPoints",numCells,sideCardinality, numBasisEPoints);
417  ScalarViewType wBasisNormalAtTargetEPoints("wBasisNormalAtTargetEPoints",numCells,sideCardinality, numTargetEPoints);
418  ScalarViewType targetNormalAtTargetEPoints("targetNormalAtTargetEPoints",numCells, numTargetEPoints);
419 
420  ordinal_type offsetBasis = basisEPointsRange(sideDim,is).first;
421  ordinal_type offsetTarget = targetEPointsRange(sideDim,is).first;
422  auto targetEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTargetEvalWeights(sideDim,is));
423  auto basisEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getBasisEvalWeights(sideDim,is));
424 
425 
426  typedef ComputeBasisCoeffsOnSides_HDiv<ScalarViewType, decltype(basisEWeights), decltype(tagToOrdinal)> functorTypeSide;
427  Kokkos::parallel_for(policy, functorTypeSide(basisNormalAtBasisEPoints, basisAtBasisEPoints,
428  basisEWeights, wBasisNormalAtBasisEPoints, targetEWeights,
429  basisAtTargetEPoints, wBasisNormalAtTargetEPoints, tagToOrdinal,
430  targetAtEPoints, targetNormalAtTargetEPoints,
431  refSidesNormal, sideCardinality, offsetBasis,
432  offsetTarget, sideDim,
433  dim, is));
434 
435  ScalarViewType sideMassMat_("sideMassMat_", numCells, sideCardinality+1, sideCardinality+1),
436  sideRhsMat_("rhsMat_", numCells, sideCardinality+1);
437 
438  ScalarViewType targetEWeights_("targetEWeights", numCells, 1, targetEWeights.extent(0));
439  RealSpaceTools<SpT>::clone(targetEWeights_, targetEWeights);
440 
441  range_type range_H(0, sideCardinality);
442  range_type range_B(sideCardinality, sideCardinality+1);
443  ScalarViewType ones("ones",numCells,1,numBasisEPoints);
444  Kokkos::deep_copy(ones,1);
445 
446  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(sideMassMat_, Kokkos::ALL(), range_H, range_H), basisNormalAtBasisEPoints, wBasisNormalAtBasisEPoints);
447  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(sideMassMat_, Kokkos::ALL(), range_H, range_B), wBasisNormalAtBasisEPoints, ones);
448 
449  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(sideRhsMat_, Kokkos::ALL(), range_H), targetNormalAtTargetEPoints, wBasisNormalAtTargetEPoints);
450  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(sideRhsMat_, Kokkos::ALL(), range_B), targetNormalAtTargetEPoints, targetEWeights_);
451 
452  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, SpT> WorkArrayViewType;
453  ScalarViewType t_("t",numCells, sideCardinality+1);
454  WorkArrayViewType w_("w",numCells, sideCardinality+1);
455 
456  auto sideDof = Kokkos::subview(tagToOrdinal, sideDim, is, Kokkos::ALL());
457 
458  ElemSystem sideSystem("sideSystem", false);
459  sideSystem.solve(basisCoeffs, sideMassMat_, sideRhsMat_, t_, w_, sideDof, sideCardinality, 1);
460  }
461 
462 
463  //Cell
464  ordinal_type numCellDofs = cellBasis->getDofCount(dim,0);
465  if(numCellDofs==0)
466  return;
467 
468  Basis<SpT,scalarType,scalarType> *hcurlBasis = NULL;
469  if(cellTopo.getKey() == shards::getCellTopologyData<shards::Hexahedron<8> >()->key)
470  hcurlBasis = new Basis_HCURL_HEX_In_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
471  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Tetrahedron<4> >()->key)
472  hcurlBasis = new Basis_HCURL_TET_In_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
473  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Quadrilateral<4> >()->key)
474  hcurlBasis = new Basis_HGRAD_QUAD_Cn_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
475  else if(cellTopo.getKey() == shards::getCellTopologyData<shards::Triangle<3> >()->key)
476  hcurlBasis = new Basis_HGRAD_TRI_Cn_FEM<SpT,scalarType,scalarType>(cellBasis->getDegree());
477  else {
478  std::stringstream ss;
479  ss << ">>> ERROR (Intrepid2::ProjectionTools::getHDivEvaluationPoints): "
480  << "Method not implemented for basis " << name;
481  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
482  }
483  if(hcurlBasis == NULL) return;
484 
485 
486  auto targetDivEPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTargetDerivPointsRange());
487  auto basisDivEPointsRange = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getBasisDerivPointsRange());
488 
489  ordinal_type numTargetDivEPoints = range_size(targetDivEPointsRange(dim,0));
490  ordinal_type numBasisDivEPoints = range_size(basisDivEPointsRange(dim,0));
491 
492  auto targetDivEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTargetDerivEvalWeights(dim,0));
493  auto divEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getBasisDerivEvalWeights(dim,0));
494 
495  ordinal_type offsetBasisDiv = basisDivEPointsRange(dim,0).first;
496  ordinal_type offsetTargetDiv = targetDivEPointsRange(dim,0).first;
497 
498  ScalarViewType weightedBasisDivAtBasisEPoints("weightedBasisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
499  ScalarViewType weightedBasisDivAtTargetEPoints("weightedBasisDivAtTargetEPoints",numCells, numCellDofs, numTargetDivEPoints);
500  ScalarViewType basisDivAtBasisEPoints("basisDivAtBasisEPoints",numCells,numCellDofs, numBasisDivEPoints);
501  ScalarViewType targetSideDivAtBasisEPoints("targetSideDivAtBasisEPoints",numCells, numBasisDivEPoints);
502 
503  auto cellDofs = Kokkos::subview(tagToOrdinal, dim, 0, Kokkos::ALL());
504  typedef ComputeBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights), decltype(computedDofs), decltype(cellDofs)> functorType;
505  Kokkos::parallel_for(policy, functorType( basisCoeffs, targetSideDivAtBasisEPoints, basisDivAtBasisEPoints,
506  basisDivAtBasisDivEPoints, divEWeights, weightedBasisDivAtBasisEPoints, targetDivEWeights, basisDivAtTargetDivEPoints, weightedBasisDivAtTargetEPoints,
507  computedDofs, cellDofs, numCellDofs, offsetBasisDiv, offsetTargetDiv, numSideDofs));
508 
509 
510  ordinal_type hcurlBasisCardinality = hcurlBasis->getCardinality();
511  ordinal_type numCurlInteriorDOFs = hcurlBasis->getDofCount(dim,0);
512 
513  range_type range_H(0, numCellDofs);
514  range_type range_B(numCellDofs, numCellDofs+numCurlInteriorDOFs);
515 
516 
517  ScalarViewType massMat_("massMat_",numCells,numCellDofs+numCurlInteriorDOFs,numCellDofs+numCurlInteriorDOFs),
518  rhsMatTrans("rhsMatTrans",numCells,numCellDofs+numCurlInteriorDOFs);
519 
520  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_H), basisDivAtBasisEPoints, weightedBasisDivAtBasisEPoints);
521  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetDivAtDivEPoints, weightedBasisDivAtTargetEPoints);
522  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_H), targetSideDivAtBasisEPoints, weightedBasisDivAtBasisEPoints,true);
523 
524  if(numCurlInteriorDOFs>0){
525  ordinal_type numTargetEPoints = range_size(targetEPointsRange(dim,0));
526  ordinal_type numBasisEPoints = range_size(basisEPointsRange(dim,0));
527 
528  auto targetEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getTargetEvalWeights(dim,0));
529  auto basisEWeights = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(),projStruct->getBasisEvalWeights(dim,0));
530 
531  ordinal_type offsetBasis = basisEPointsRange(dim,0).first;
532 
533  auto basisEPoints = Kokkos::subview(basisEPoints_, 0, basisEPointsRange(dim,0), Kokkos::ALL());
534 
535  ScalarViewType negPartialProjAtBasisEPoints("targetSideAtBasisEPoints",numCells, numBasisEPoints, dim);
536  ScalarViewType nonWeightedBasisAtBasisEPoints("basisAtBasisEPoints",numCells,numCellDofs, numBasisEPoints, dim);
537  ScalarViewType hcurlBasisCurlAtBasisEPoints("hcurlBasisCurlAtBasisEPoints",hcurlBasisCardinality, numBasisEPoints,dim);
538  ScalarViewType hcurlBasisCurlAtTargetEPoints("hcurlBasisCurlAtTargetEPoints", hcurlBasisCardinality,numTargetEPoints, dim);
539  ScalarViewType wHcurlBasisCurlAtBasisEPoints("wHcurlBasisHcurlAtBasisEPoints", numCells, numCurlInteriorDOFs, numBasisEPoints,dim);
540  ScalarViewType wHcurlBasisCurlAtTargetEPoints("wHcurlBasisHcurlAtTargetEPoints",numCells, numCurlInteriorDOFs, numTargetEPoints,dim);
541 
542  hcurlBasis->getValues(hcurlBasisCurlAtBasisEPoints, basisEPoints, OPERATOR_CURL);
543  hcurlBasis->getValues(hcurlBasisCurlAtTargetEPoints, Kokkos::subview(targetEPoints,0,targetEPointsRange(dim,0),Kokkos::ALL()), OPERATOR_CURL);
544 
545  auto hCurlTagToOrdinal = Kokkos::create_mirror_view_and_copy(typename SpT::memory_space(), hcurlBasis->getAllDofOrdinal());
546  auto cellHCurlDof = Kokkos::subview(hCurlTagToOrdinal, dim, 0, range_type(0, numCurlInteriorDOFs));
547 
548  typedef ComputeHCurlBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights),
549  decltype(tagToOrdinal), decltype(computedDofs), decltype(cellDofs)> functorTypeHCurlCells;
550  Kokkos::parallel_for(policy, functorTypeHCurlCells(basisCoeffs, negPartialProjAtBasisEPoints, nonWeightedBasisAtBasisEPoints,
551  basisAtBasisEPoints, hcurlBasisCurlAtBasisEPoints, basisEWeights, wHcurlBasisCurlAtBasisEPoints, targetEWeights,
552  hcurlBasisCurlAtTargetEPoints, wHcurlBasisCurlAtTargetEPoints, tagToOrdinal, computedDofs, cellHCurlDof,
553  numCellDofs, offsetBasis, numSideDofs, dim));
554 
555  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(massMat_, Kokkos::ALL(), range_H,range_B), nonWeightedBasisAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints);
556  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), Kokkos::subview(targetAtEPoints, Kokkos::ALL(), targetEPointsRange(dim,0), Kokkos::ALL()), wHcurlBasisCurlAtTargetEPoints);
557  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(rhsMatTrans, Kokkos::ALL(), range_B), negPartialProjAtBasisEPoints, wHcurlBasisCurlAtBasisEPoints,true);
558  }
559  delete hcurlBasis;
560 
561  typedef Kokkos::DynRankView<scalarType, Kokkos::LayoutRight, SpT> WorkArrayViewType;
562  ScalarViewType t_("t",numCells, numCellDofs+numCurlInteriorDOFs);
563  WorkArrayViewType w_("w",numCells, numCellDofs+numCurlInteriorDOFs);
564 
565  ElemSystem cellSystem("cellSystem", true);
566  cellSystem.solve(basisCoeffs, massMat_, rhsMatTrans, t_, w_, cellDofs, numCellDofs, numCurlInteriorDOFs);
567 }
568 
569 
570 
571 }
572 }
573 
574 #endif
575 
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< ExecSpaceType, typename BasisType::scalarType > *projStruct)
Computes the basis coefficients of the HDiv projection of the target function.
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
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.
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
static void modifyBasisByOrientation(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input, const Kokkos::DynRankView< ortValueType, ortProperties...> orts, const BasisType *basis)
Modify basis due to orientation.
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 mapToReferenceSubcell(Kokkos::DynRankView< refSubcellPointValueType, refSubcellPointProperties...> refSubcellPoints, const Kokkos::DynRankView< paramPointValueType, paramPointProperties...> paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
ordinal_type getMaxNumEvalPoints(const EvalPointsType type) const
Returns the maximum number of evaluation points across all the subcells.
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 getHDivEvaluationPoints(typename BasisType::ScalarViewType evaluationPoints, typename BasisType::ScalarViewType divEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties...> cellOrientations, const BasisType *cellBasis, ProjectionStruct< ExecSpaceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=EvalPointsType::TARGET)
Computes evaluation points for HDiv projection.
Header file for the Intrepid2::FunctionSpaceTools class.
static void clone(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input)
Clone input array.
Header file for Intrepid2::ArrayTools class providing utilities for array operations.
view_type getTargetEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the points where to evaluate the target function on a subcell.
static void mapToModifiedReference(outPointViewType outPoints, const refPointViewType refPoints, const shards::CellTopology cellTopo, const ordinal_type cellOrt=0)
Computes modified parameterization maps of 1- and 2-subcells with orientation.
static void setSubcellParametrization()
Defines orientation-preserving parametrizations of reference edges and faces of cell topologies with ...
An helper class to compute the evaluation points and weights needed for performing projections...
static void getSubcellParametrization(subcellParamViewType &subcellParam, const ordinal_type subcellDim, const shards::CellTopology parentCell)
Returns array with the coefficients of the parametrization maps for the edges or faces of a reference...
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.