17 #ifndef __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
18 #define __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
27 template<
typename po
intValueType>
28 KOKKOS_INLINE_FUNCTION
31 const pointValueType y,
32 const pointValueType z ) {
43 if ( (0.0 <= xyz && xyz <= 0.5) &&
44 (0.0 <= x && x <= 0.5) &&
45 (0.0 <= y && y <= 0.5) &&
46 (0.0 <= z && z <= 0.5) )
50 if ( (0.5 <= xyz && xyz <= 1.0) &&
51 (0.5 <= x && x <= 1.0) &&
52 (0.0 <= y && y <= 0.5) &&
53 (0.0 <= z && z <= 0.5) )
57 if ( (0.5 <= xyz && xyz <= 1.0) &&
58 (0.0 <= x && x <= 0.5) &&
59 (0.5 <= y && y <= 1.0) &&
60 (0.0 <= z && z <= 0.5) )
64 if ( (0.5 <= xyz && xyz <= 1.0) &&
65 (0.0 <= x && x <= 0.5) &&
66 (0.0 <= y && y <= 0.5) &&
67 (0.5 <= z && z <= 1.0) )
71 if ( (0.0 <= yz && yz <= 0.5) &&
72 (0.5 <= xy && xy <= 1.0) &&
73 (0.5 <= xz && xz <= 1.0) &&
74 (0.0 <= x && x <= 0.5) )
78 if ( (0.5 <= xy && xy <= 1.0) &&
79 (0.5 <= yz && yz <= 1.0) &&
80 (0.5 <= xz && xz <= 1.0) &&
81 (0.75 <= xyz && xyz <= 1.0) )
85 if ( (0.5 <= yz && yz <= 1.0) &&
86 (0.0 <= xy && xy <= 0.5) &&
87 (0.5 <= xz && xz <= 1.0) &&
88 (0.0 <= z && z <= 0.5) )
92 if ( (0.0 <= yz && yz <= 0.5) &&
93 (0.0 <= xy && xy <= 0.5) &&
94 (0.5 <= xz && xz <= 1.0) &&
95 (0.0 <= y && y <= 0.25) )
99 if ( (0.0 <= xz && xz <= 0.5) &&
100 (0.0 <= yz && yz <= 0.5) &&
101 (0.5 <= xy && xy <= 1.0) &&
102 (0.0 <= z && z <= 0.25) )
106 if ( (0.0 <= xz && xz <= 0.5) &&
107 (0.5 <= xy && xy <= 1.0) &&
108 (0.5 <= yz && yz <= 1.0) &&
109 (0.0 <= y && y <= 0.5) )
113 if ( (0.0 <= xz && xz <= 0.5) &&
114 (0.5 <= yz && yz <= 1.0) &&
115 (0.0 <= xy && xy <= 0.5) &&
116 (0.0 <= x && x <= 0.25) )
120 if ( (0.5 <= xyz && xyz <= 0.75) &&
121 (0.0 <= xz && xz <= 0.5) &&
122 (0.0 <= yz && yz <= 0.5) &&
123 (0.0 <= xy && xy <= 0.5) )
129 template<EOperator opType>
130 template<
typename outputValueViewType,
131 typename inputPointViewType>
132 KOKKOS_INLINE_FUNCTION
136 const inputPointViewType input ) {
138 case OPERATOR_VALUE: {
139 const typename inputPointViewType::value_type r = input(0);
140 const typename inputPointViewType::value_type s = input(1);
141 const typename inputPointViewType::value_type t = input(2);
144 for (
auto i=0;i<10;++i)
145 output.access(i) = 0.0;
147 const auto subtet = getLocalSubTet( r, s, t );
153 typename inputPointViewType::value_type aux = 0.0;
156 output.access(0) = 1. - 2. * (r + s + t);
157 output.access(4) = 2. * r;
158 output.access(6) = 2. * s;
159 output.access(7) = 2. * t;
162 output.access(1) = 2. * r - 1.;
163 output.access(4) = 2. - 2. * (r + s + t);
164 output.access(5) = 2. * s;
165 output.access(8) = 2. * t;
168 output.access(2) = 2. * s - 1.;
169 output.access(5) = 2. * r;
170 output.access(6) = 2. - 2. * (r + s + t);
171 output.access(9) = 2. * t;
174 output.access(3) = 2. * t - 1.;
175 output.access(7) = 2. - 2. * (r + s + t);
176 output.access(8) = 2. * r;
177 output.access(9) = 2. * s;
180 output.access(4) = 1. - 2. * (s + t);
181 output.access(5) = 2. * (r + s) - 1.;
182 output.access(8) = 2. * (r + t) - 1.;
186 output.access(5) = 2. * (r + s) - 1.;
187 output.access(8) = 2. * (r + t) - 1.;
188 output.access(9) = 2. * (s + t) - 1.;
189 aux = 4. - 4. * (r + s + t);
192 output.access(7) = 1. - 2. * (r + s);
193 output.access(8) = 2. * (r + t) - 1.;
194 output.access(9) = 2. * (s + t) - 1.;
198 output.access(4) = 1. - 2. * (s + t);
199 output.access(7) = 1. - 2. * (r + s);
200 output.access(8) = 2. * (r + t) - 1.;
204 output.access(4) = 1. - 2. * (s + t);
205 output.access(5) = 2. * (r + s) - 1.;
206 output.access(6) = 1. - 2. * (r + t);
210 output.access(5) = 2. * (r + s) - 1.;
211 output.access(6) = 1. - 2. * (r + t);
212 output.access(9) = 2. * (s + t) - 1.;
216 output.access(6) = 1. - 2. * (r + t);
217 output.access(7) = 1. - 2. * (r + s);
218 output.access(9) = 2. * (s + t) - 1.;
222 output.access(4) = 1. - 2. * (s + t);
223 output.access(6) = 1. - 2. * (r + t);
224 output.access(7) = 1. - 2. * (r + s);
225 aux = 4. * (r + s + t) - 2.;
228 for (
auto i=4;i<10;++i)
229 output.access(i) += aux/6.0;
233 case OPERATOR_GRAD: {
234 const typename inputPointViewType::value_type r = input(0);
235 const typename inputPointViewType::value_type s = input(1);
236 const typename inputPointViewType::value_type t = input(2);
238 output.access(0,0) = (-17 + 20*r + 20*s + 20*t)/8.;
239 output.access(0,1) = (-17 + 20*r + 20*s + 20*t)/8.;
240 output.access(0,2) = (-17 + 20*r + 20*s + 20*t)/8.;
241 output.access(1,0) = -0.375 + (5*r)/2.;
242 output.access(1,1) = 0.;
243 output.access(1,2) = 0.;
244 output.access(2,0) = 0.;
245 output.access(2,1) = -0.375 + (5*s)/2.;
246 output.access(2,2) = 0.;
247 output.access(3,0) = 0.;
248 output.access(3,1) = 0.;
249 output.access(3,2) = -0.375 + (5*t)/2.;
250 output.access(4,0) = (-35*(-1 + 2*r + s + t))/12.;
251 output.access(4,1) = (-4 - 35*r + 5*s + 10*t)/12.;
252 output.access(4,2) = (-4 - 35*r + 10*s + 5*t)/12.;
253 output.access(5,0) = (-1 + 5*r + 40*s - 5*t)/12.;
254 output.access(5,1) = (-1 + 40*r + 5*s - 5*t)/12.;
255 output.access(5,2) = (-5*(-1 + r + s + 2*t))/12.;
256 output.access(6,0) = (-4 + 5*r - 35*s + 10*t)/12.;
257 output.access(6,1) = (-35*(-1 + r + 2*s + t))/12.;
258 output.access(6,2) = (-4 + 10*r - 35*s + 5*t)/12.;
259 output.access(7,0) = (-4 + 5*r + 10*s - 35*t)/12.;
260 output.access(7,1) = (-4 + 10*r + 5*s - 35*t)/12.;
261 output.access(7,2) = (-35*(-1 + r + s + 2*t))/12.;
262 output.access(8,0) = (-1 + 5*r - 5*s + 40*t)/12.;
263 output.access(8,1) = (-5*(-1 + r + 2*s + t))/12.;
264 output.access(8,2) = (-1 + 40*r - 5*s + 5*t)/12.;
265 output.access(9,0) = (-5*(-1 + 2*r + s + t))/12.;
266 output.access(9,1) = (-1 - 5*r + 5*s + 40*t)/12.;
267 output.access(9,2) = (-1 - 5*r + 40*s + 5*t)/12.;
271 const ordinal_type jend = output.extent(1);
272 const ordinal_type iend = output.extent(0);
274 for (ordinal_type j=0;j<jend;++j)
275 for (
auto i=0;i<iend;++i)
276 output.access(i, j) = 0.0;
280 INTREPID2_TEST_FOR_ABORT(
true ,
281 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
286 template<
typename DT,
287 typename outputValueValueType,
class ...outputValueProperties,
288 typename inputPointValueType,
class ...inputPointProperties>
290 Basis_HGRAD_TET_COMP12_FEM::
291 getValues(
const typename DT::execution_space& space,
292 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
293 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
294 const EOperator operatorType ) {
295 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
296 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
297 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
300 const auto loopSize = inputPoints.extent(0);
301 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
303 switch (operatorType) {
305 case OPERATOR_VALUE: {
306 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
307 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
312 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
313 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
316 case OPERATOR_CURL: {
317 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
318 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
322 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
323 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
335 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
336 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
340 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
341 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
349 template<
typename DT,
typename OT,
typename PT>
352 const ordinal_type spaceDim = 3;
353 this->basisCardinality_ = 10;
354 this->basisDegree_ = 1;
355 this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
356 this->basisType_ = BASIS_FEM_DEFAULT;
357 this->basisCoordinates_ = COORDINATES_CARTESIAN;
358 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
362 const ordinal_type tagSize = 4;
363 const ordinal_type posScDim = 0;
364 const ordinal_type posScOrd = 1;
365 const ordinal_type posDfOrd = 2;
368 ordinal_type tags[] = { 0, 0, 0, 1,
383 this->setOrdinalTagData(this->tagToOrdinal_,
386 this->basisCardinality_,
394 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
395 dofCoords(
"dofCoordsHost", this->basisCardinality_, spaceDim);
397 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
398 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
399 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
400 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
401 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
402 dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
403 dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
404 dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
405 dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
406 dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
408 this->dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoords);
409 Kokkos::deep_copy(this->dofCoords_, dofCoords);
412 template<
typename DT,
typename OT,
typename PT>
415 ordinal_type& perTeamSpaceSize,
416 ordinal_type& perThreadSpaceSize,
418 const EOperator operatorType)
const {
419 perTeamSpaceSize = 0;
420 perThreadSpaceSize = 0;
423 template<
typename DT,
typename OT,
typename PT>
424 KOKKOS_INLINE_FUNCTION
427 OutputViewType outputValues,
428 const PointViewType inputPoints,
429 const EOperator operatorType,
430 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
431 const typename DT::execution_space::scratch_memory_space & scratchStorage,
432 const ordinal_type subcellDim,
433 const ordinal_type subcellOrdinal)
const {
435 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
436 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_COMP12_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
438 (void) scratchStorage;
440 const int numPoints = inputPoints.extent(0);
442 switch(operatorType) {
444 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
445 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
446 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
451 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
452 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
453 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
454 Impl::Basis_HGRAD_TET_COMP12_FEM::Serial<OPERATOR_GRAD>::getValues( output, input);
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
FEM basis evaluation on a reference Tetrahedron cell.
Basis_HGRAD_TET_COMP12_FEM()
Constructor.
static KOKKOS_INLINE_FUNCTION ordinal_type getLocalSubTet(const pointValueType x, const pointValueType y, const pointValueType z)
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.
virtual void getScratchSpaceSize(ordinal_type &perTeamSpaceSize, ordinal_type &perThreadSpaceSize, const PointViewType inputPointsconst, const EOperator operatorType=OPERATOR_VALUE) const override
Return the size of the scratch space, in bytes, needed for using the team-level implementation of get...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.