Intrepid2
Intrepid2_HCURL_TET_I1_FEMDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_HCURL_TET_I1_FEM_DEF_HPP__
17#define __INTREPID2_HCURL_TET_I1_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22 namespace Impl {
23
24 template<EOperator opType>
25 template<typename OutputViewType,
26 typename inputViewType>
27 KOKKOS_INLINE_FUNCTION
28 void
29 Basis_HCURL_TET_I1_FEM::Serial<opType>::
30 getValues( OutputViewType output,
31 const inputViewType input ) {
32 switch (opType) {
33 case OPERATOR_VALUE: {
34 const auto x = input(0);
35 const auto y = input(1);
36 const auto z = input(2);
37
38 // output is subview of a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim), dim0 is iteration from range
39 output.access(0, 0) = 2.0*(1.0 - y - z);
40 output.access(0, 1) = 2.0*x;
41 output.access(0, 2) = 2.0*x;
42
43 output.access(1, 0) =-2.0*y;
44 output.access(1, 1) = 2.0*x;
45 output.access(1, 2) = 0.0;
46
47 output.access(2, 0) = -2.0*y;
48 output.access(2, 1) = 2.0*(-1.0 + x + z);
49 output.access(2, 2) = -2.0*y;
50
51 output.access(3, 0) = 2.0*z;
52 output.access(3, 1) = 2.0*z;
53 output.access(3, 2) = 2.0*(1.0 - x - y);
54
55 output.access(4, 0) =-2.0*z;
56 output.access(4, 1) = 0.0;
57 output.access(4, 2) = 2.0*x;
58
59 output.access(5, 0) = 0.0;
60 output.access(5, 1) =-2.0*z;
61 output.access(5, 2) = 2.0*y;
62 break;
63 }
64 case OPERATOR_CURL: {
65 // output is subview of a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim), dim0 is iteration from range
66 output.access(0, 0) = 0.0;
67 output.access(0, 1) =-4.0;
68 output.access(0, 2) = 4.0;
69
70 output.access(1, 0) = 0.0;
71 output.access(1, 1) = 0.0;
72 output.access(1, 2) = 4.0;
73
74 output.access(2, 0) =-4.0;
75 output.access(2, 1) = 0.0;
76 output.access(2, 2) = 4.0;
77
78 output.access(3, 0) =-4.0;
79 output.access(3, 1) = 4.0;
80 output.access(3, 2) = 0.0;
81
82 output.access(4, 0) = 0.0;
83 output.access(4, 1) =-4.0;
84 output.access(4, 2) = 0.0;
85
86 output.access(5, 0) = 4.0;
87 output.access(5, 1) = 0.0;
88 output.access(5, 2) = 0.0;
89 break;
90 }
91 default: {
92 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
93 opType != OPERATOR_CURL,
94 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_I1_FEM::Serial::getValues) operator is not supported");
95 }
96 }
97 }
98
99 template<typename DT,
100 typename outputValueValueType, class ...outputValueProperties,
101 typename inputPointValueType, class ...inputPointProperties>
102 void
103 Basis_HCURL_TET_I1_FEM::
104 getValues( const typename DT::execution_space& space,
105 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
106 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
107 const EOperator operatorType ) {
108 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
109 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
110 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
111
112 // Number of evaluation points = dim 0 of inputPoints
113 const auto loopSize = inputPoints.extent(0);
114 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
115
116 switch (operatorType) {
117
118 case OPERATOR_VALUE: {
119 typedef Functor<outputValueViewType, inputPointViewType, OPERATOR_VALUE> FunctorType;
120 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
121 break;
122 }
123 case OPERATOR_CURL: {
124 typedef Functor<outputValueViewType, inputPointViewType, OPERATOR_CURL> FunctorType;
125 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
126 break;
127 }
128 case OPERATOR_DIV: {
129 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
130 ">>> ERROR (Basis_HCURL_TET_I1_FEM): DIV is invalid operator for HCURL Basis Functions");
131 break;
132 }
133 case OPERATOR_GRAD: {
134 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
135 ">>> ERROR (Basis_HCURL_TET_I1_FEM): GRAD is invalid operator for HCURL Basis Functions");
136 break;
137 }
138 case OPERATOR_D1:
139 case OPERATOR_D2:
140 case OPERATOR_D3:
141 case OPERATOR_D4:
142 case OPERATOR_D5:
143 case OPERATOR_D6:
144 case OPERATOR_D7:
145 case OPERATOR_D8:
146 case OPERATOR_D9:
147 case OPERATOR_D10: {
148 INTREPID2_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
149 (operatorType == OPERATOR_D2) ||
150 (operatorType == OPERATOR_D3) ||
151 (operatorType == OPERATOR_D4) ||
152 (operatorType == OPERATOR_D5) ||
153 (operatorType == OPERATOR_D6) ||
154 (operatorType == OPERATOR_D7) ||
155 (operatorType == OPERATOR_D8) ||
156 (operatorType == OPERATOR_D9) ||
157 (operatorType == OPERATOR_D10) ),
158 std::invalid_argument,
159 ">>> ERROR (Basis_HCURL_TET_I1_FEM): Invalid operator type");
160 break;
161 }
162 default: {
163 INTREPID2_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
164 (operatorType != OPERATOR_GRAD) &&
165 (operatorType != OPERATOR_CURL) &&
166 (operatorType != OPERATOR_DIV) &&
167 (operatorType != OPERATOR_D1) &&
168 (operatorType != OPERATOR_D2) &&
169 (operatorType != OPERATOR_D3) &&
170 (operatorType != OPERATOR_D4) &&
171 (operatorType != OPERATOR_D5) &&
172 (operatorType != OPERATOR_D6) &&
173 (operatorType != OPERATOR_D7) &&
174 (operatorType != OPERATOR_D8) &&
175 (operatorType != OPERATOR_D9) &&
176 (operatorType != OPERATOR_D10) ),
177 std::invalid_argument,
178 ">>> ERROR (Basis_HCURL_TET_I1_FEM): Invalid operator type");
179 }
180 }
181 }
182 }
183
184
185 // -------------------------------------------------------------------------------------
186
187 template< typename DT, typename OT, typename PT >
190 const ordinal_type spaceDim = 3;
191 this->basisCardinality_ = 6;
192 this->basisDegree_ = 1;
193 this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
194 this->basisType_ = BASIS_FEM_DEFAULT;
195 this->basisCoordinates_ = COORDINATES_CARTESIAN;
196 this->functionSpace_ = FUNCTION_SPACE_HCURL;
197
198 // initialize tags
199 {
200 // Basis-dependent intializations
201 const ordinal_type tagSize = 4; // size of DoF tag
202 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
203 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
204 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
205
206 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
207 ordinal_type tags[24] = {
208 1, 0, 0, 1,
209 1, 1, 0, 1,
210 1, 2, 0, 1,
211 1, 3, 0, 1,
212 1, 4, 0, 1,
213 1, 5, 0, 1 };
214
215 //host tags
216 OrdinalTypeArray1DHost tagView(&tags[0], 24);
217
218 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
219 this->setOrdinalTagData(this->tagToOrdinal_,
220 this->ordinalToTag_,
221 tagView,
222 this->basisCardinality_,
223 tagSize,
224 posScDim,
225 posScOrd,
226 posDfOrd);
227 }
228 // dofCoords on host and create its mirror view to device
229 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
230 dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
231
232 dofCoords(0,0) = 0.5; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
233 dofCoords(1,0) = 0.5; dofCoords(1,1) = 0.5; dofCoords(1,2) = 0.0;
234 dofCoords(2,0) = 0.0; dofCoords(2,1) = 0.5; dofCoords(2,2) = 0.0;
235 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 0.5;
236 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.5;
237 dofCoords(5,0) = 0.0; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.5;
238
239 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
240 Kokkos::deep_copy(this->dofCoords_, dofCoords);
241
242 // dofCoeffs on host and create its mirror view to device
243 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
244 dofCoeffs("dofCoeffsHost", this->basisCardinality_,spaceDim);
245
246 dofCoeffs(0,0) = 0.5; dofCoeffs(0,1) = 0.0; dofCoeffs(0,2) = 0.0;
247 dofCoeffs(1,0) = -0.5; dofCoeffs(1,1) = 0.5; dofCoeffs(1,2) = 0.0;
248 dofCoeffs(2,0) = 0.0; dofCoeffs(2,1) = -0.5; dofCoeffs(2,2) = 0.0;
249 dofCoeffs(3,0) = 0.0; dofCoeffs(3,1) = 0.0; dofCoeffs(3,2) = 0.5;
250 dofCoeffs(4,0) = -0.5; dofCoeffs(4,1) = 0.0; dofCoeffs(4,2) = 0.5;
251 dofCoeffs(5,0) = 0.0; dofCoeffs(5,1) = -0.5; dofCoeffs(5,2) = 0.5;
252
253 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffs);
254 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
255
256 }
257
258 template<typename DT, typename OT, typename PT>
259 void
261 ordinal_type& perTeamSpaceSize,
262 ordinal_type& perThreadSpaceSize,
263 const PointViewType inputPoints,
264 const EOperator operatorType) const {
265 perTeamSpaceSize = 0;
266 perThreadSpaceSize = 0;
267 }
268
269 template<typename DT, typename OT, typename PT>
270 KOKKOS_INLINE_FUNCTION
271 void
273 OutputViewType outputValues,
274 const PointViewType inputPoints,
275 const EOperator operatorType,
276 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
277 const typename DT::execution_space::scratch_memory_space & scratchStorage,
278 const ordinal_type subcellDim,
279 const ordinal_type subcellOrdinal) const {
280
281 INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
282 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_I1_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
283
284 (void) scratchStorage; //avoid unused variable warning
285
286 const int numPoints = inputPoints.extent(0);
287
288 switch(operatorType) {
289 case OPERATOR_VALUE:
290 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
291 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
292 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
294 });
295 break;
296 case OPERATOR_CURL:
297 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
298 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
299 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
300 Impl::Basis_HCURL_TET_I1_FEM::Serial<OPERATOR_CURL>::getValues( output, input);
301 });
302 break;
303 default: {
304 INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR: (Intrepid2::Basis_HCURL_TET_I1_FEM::getValues), Operator Type not supported.");
305 }
306 }
307 }
308
309}// namespace Intrepid2
310
311#endif
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
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::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.