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
55
56
57namespace Intrepid2 {
58namespace Experimental {
59
60
61template<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
125template<typename ViewType1, typename ViewType2, typename ViewType3,
126typename 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
177template<typename ViewType1, typename ViewType2, typename ViewType3,
178typename 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
245template<typename DeviceType>
246template<typename BasisType,
247typename ortValueType, class ...ortProperties>
248void
249ProjectionTools<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
311template<typename DeviceType>
312template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
313typename funValsValueType, class ...funValsProperties,
314typename BasisType,
315typename ortValueType,class ...ortProperties>
316void
317ProjectionTools<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
418 typedef ComputeBasisCoeffsOnSides_HDiv<ScalarViewType, decltype(basisEWeights), decltype(tagToOrdinal), decltype(targetAtEPoints)> functorTypeSide;
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
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());
496 typedef ComputeBasisCoeffsOnCells_HDiv<decltype(basisCoeffs), ScalarViewType, decltype(divEWeights), decltype(computedDofs), decltype(cellDofs)> functorType;
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 Intrepid2::ArrayTools class providing utilities for array operations.
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
Header file for the Intrepid2::FunctionSpaceTools class.
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
Implementation of the default H(grad)-compatible FEM basis of degree n on Quadrilateral cell Implemen...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
ordinal_type getCardinality() const
Returns cardinality of the basis.
const OrdinalTypeArray3DHost getAllDofOrdinal() const
DoF tag to ordinal data structure.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
static void 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.
An helper class to compute the evaluation points and weights needed for performing projections.
ordinal_type getNumBasisEvalPoints()
Returns number of basis evaluation points.
view_type getDerivEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the evaluation points for basis/target derivatives on a subcell.
const range_tag getDerivPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target derivative evaluation points on subcells.
const range_tag getBasisPointsRange() const
Returns the range tag of the basis evaluation points subcells.
const range_tag getBasisDerivPointsRange() const
Returns the range tag of the derivative evaluation points on subcell.
const range_tag getTargetPointsRange() const
Returns the range of the target function evaluation points on subcells.
view_type getEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId, EvalPointsType type) const
Returns the basis/target evaluation points on a subcell.
const key_tag getTopologyKey() const
Returns the key tag for subcells.
view_type getTargetDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function derivatives evaluation weights on a subcell.
view_type getBasisEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis evaluation weights on a subcell.
const range_tag getTargetDerivPointsRange() const
Returns the range tag of the target function derivative evaluation points on subcells.
ordinal_type getNumBasisDerivEvalPoints()
Returns number of evaluation points for basis derivatives.
view_type getTargetEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the function evaluation weights on a subcell.
view_type getBasisDerivEvalWeights(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the basis derivatives evaluation weights on a subcell.
view_type getTargetEvalPoints(const ordinal_type subCellDim, const ordinal_type subCellId)
Returns the points where to evaluate the target function on a subcell.
const range_tag getPointsRange(const EvalPointsType type) const
Returns the range tag of the basis/target evaluation points in subcells.
static void 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.
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.
static void integrate(Kokkos::DynRankView< outputValueValueType, outputValueProperties... > outputValues, const Kokkos::DynRankView< leftValueValueType, leftValueProperties... > leftValues, const Kokkos::DynRankView< rightValueValueType, rightValueProperties... > rightValues, const bool sumInto=false)
Contracts leftValues and rightValues arrays on the point and possibly space dimensions and stores the...
static KOKKOS_INLINE_FUNCTION void mapSubcellCoordsToRefCell(coordsViewType cellCoords, const subcellCoordsViewType subCellCoords, const ParamViewType subcellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Maps points defined on the subCell manifold into the parent Cell.
static void modifyBasisByOrientation(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input, const OrientationViewType orts, const BasisType *basis)
Modify basis due to orientation.
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
Class to solve a square system A x = b on each cell A is expected to be saddle a point (KKT) matrix o...
void solve(ViewType1 basisCoeffs, ViewType2 elemMat, ViewType2 elemRhs, ViewType2 tau, ViewType3 w, const ViewType4 elemDof, ordinal_type n, ordinal_type m=0)
Solve the system and returns the basis coefficients solve the system either using Kokkos Kernel QR or...