Intrepid2
Intrepid2_HGRAD_TET_C2_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_HGRAD_TET_C2_FEM_DEF_HPP__
17#define __INTREPID2_HGRAD_TET_C2_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22
23 namespace Impl {
24
25 template<EOperator opType>
26 template<typename OutputViewType,
27 typename inputViewType>
28 KOKKOS_INLINE_FUNCTION
29 void
30 Basis_HGRAD_TET_C2_FEM::Serial<opType>::
31 getValues( OutputViewType output,
32 const inputViewType input ) {
33 switch (opType) {
34 case OPERATOR_VALUE: {
35 const auto x = input(0);
36 const auto y = input(1);
37 const auto z = input(2);
38
39 // output is a rank-2 array with dimensions (basisCardinality_, dim0)
40 output.access(0) = (-1. + x + y + z)*(-1. + 2.*x + 2.*y + 2.*z);
41 output.access(1) = x*(-1. + 2.*x);
42 output.access(2) = y*(-1. + 2.*y);
43 output.access(3) = z*(-1. + 2.*z);
44
45 output.access(4) = -4.*x*(-1. + x + y + z);
46 output.access(5) = 4.*x*y;
47 output.access(6) = -4.*y*(-1. + x + y + z);
48 output.access(7) = -4.*z*(-1. + x + y + z);
49 output.access(8) = 4.*x*z;
50 output.access(9) = 4.*y*z;
51 break;
52 }
53 case OPERATOR_D1:
54 case OPERATOR_GRAD: {
55 const auto x = input(0);
56 const auto y = input(1);
57 const auto z = input(2);
58
59 // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
60 output.access(0, 0) = -3.+ 4.*x + 4.*y + 4.*z;
61 output.access(0, 1) = -3.+ 4.*x + 4.*y + 4.*z;
62 output.access(0, 2) = -3.+ 4.*x + 4.*y + 4.*z;
63
64 output.access(1, 0) = -1.+ 4.*x;
65 output.access(1, 1) = 0.;
66 output.access(1, 2) = 0.;
67
68 output.access(2, 0) = 0.;
69 output.access(2, 1) = -1.+ 4.*y;
70 output.access(2, 2) = 0.;
71
72 output.access(3, 0) = 0.;
73 output.access(3, 1) = 0.;
74 output.access(3, 2) = -1.+ 4.*z;
75
76
77 output.access(4, 0) = -4.*(-1.+ 2*x + y + z);
78 output.access(4, 1) = -4.*x;
79 output.access(4, 2) = -4.*x;
80
81 output.access(5, 0) = 4.*y;
82 output.access(5, 1) = 4.*x;
83 output.access(5, 2) = 0.;
84
85 output.access(6, 0) = -4.*y;
86 output.access(6, 1) = -4.*(-1.+ x + 2*y + z);
87 output.access(6, 2) = -4.*y;
88
89 output.access(7, 0) = -4.*z;
90 output.access(7, 1) = -4.*z;
91 output.access(7, 2) = -4.*(-1.+ x + y + 2*z);
92
93 output.access(8, 0) = 4.*z;
94 output.access(8, 1) = 0.;
95 output.access(8, 2) = 4.*x;
96
97 output.access(9, 0) = 0.;
98 output.access(9, 1) = 4.*z;
99 output.access(9, 2) = 4.*y;
100 break;
101 }
102 case OPERATOR_D2: {
103 output.access(0, 0) = 4.;
104 output.access(0, 1) = 4.;
105 output.access(0, 2) = 4.;
106 output.access(0, 3) = 4.;
107 output.access(0, 4) = 4.;
108 output.access(0, 5) = 4.;
109
110 output.access(1, 0) = 4.;
111 output.access(1, 1) = 0.;
112 output.access(1, 2) = 0.;
113 output.access(1, 3) = 0.;
114 output.access(1, 4) = 0.;
115 output.access(1, 5) = 0.;
116
117 output.access(2, 0) = 0.;
118 output.access(2, 1) = 0.;
119 output.access(2, 2) = 0.;
120 output.access(2, 3) = 4.;
121 output.access(2, 4) = 0.;
122 output.access(2, 5) = 0.;
123
124 output.access(3, 0) = 0.;
125 output.access(3, 1) = 0.;
126 output.access(3, 2) = 0.;
127 output.access(3, 3) = 0.;
128 output.access(3, 4) = 0.;
129 output.access(3, 5) = 4.;
130
131 output.access(4, 0) = -8.;
132 output.access(4, 1) = -4.;
133 output.access(4, 2) = -4.;
134 output.access(4, 3) = 0.;
135 output.access(4, 4) = 0.;
136 output.access(4, 5) = 0.;
137
138 output.access(5, 0) = 0.;
139 output.access(5, 1) = 4.;
140 output.access(5, 2) = 0.;
141 output.access(5, 3) = 0.;
142 output.access(5, 4) = 0.;
143 output.access(5, 5) = 0.;
144
145 output.access(6, 0) = 0.;
146 output.access(6, 1) = -4.;
147 output.access(6, 2) = 0.;
148 output.access(6, 3) = -8.;
149 output.access(6, 4) = -4.;
150 output.access(6, 5) = 0;
151
152 output.access(7, 0) = 0.;
153 output.access(7, 1) = 0.;
154 output.access(7, 2) = -4.;
155 output.access(7, 3) = 0.;
156 output.access(7, 4) = -4.;
157 output.access(7, 5) = -8.;
158
159 output.access(8, 0) = 0.;
160 output.access(8, 1) = 0.;
161 output.access(8, 2) = 4.;
162 output.access(8, 3) = 0.;
163 output.access(8, 4) = 0.;
164 output.access(8, 5) = 0.;
165
166 output.access(9, 0) = 0.;
167 output.access(9, 1) = 0.;
168 output.access(9, 2) = 0.;
169 output.access(9, 3) = 0.;
170 output.access(9, 4) = 4.;
171 output.access(9, 5) = 0.;
172 break;
173 }
174 case OPERATOR_MAX: {
175 const ordinal_type jend = output.extent(1);
176 const ordinal_type iend = output.extent(0);
177
178 for (ordinal_type j=0;j<jend;++j)
179 for (ordinal_type i=0;i<iend;++i)
180 output.access(i, j) = 0.0;
181 break;
182 }
183 default: {
184 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
185 opType != OPERATOR_GRAD &&
186 opType != OPERATOR_D1 &&
187 opType != OPERATOR_D2 &&
188 opType != OPERATOR_MAX,
189 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_C2_FEM::Serial::getValues) operator is not supported");
190 }
191 }
192 }
193
194 template<typename DT,
195 typename outputValueValueType, class ...outputValueProperties,
196 typename inputPointValueType, class ...inputPointProperties>
197 void
198 Basis_HGRAD_TET_C2_FEM::
199 getValues( const typename DT::execution_space& space,
200 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
201 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
202 const EOperator operatorType ) {
203 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
204 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
205 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
206
207 // Number of evaluation points = dim 0 of inputPoints
208 const auto loopSize = inputPoints.extent(0);
209 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
210
211 switch (operatorType) {
212
213 case OPERATOR_VALUE: {
214 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
215 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
216 break;
217 }
218 case OPERATOR_GRAD:
219 case OPERATOR_D1: {
220 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
221 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
222 break;
223 }
224 case OPERATOR_CURL: {
225 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
226 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
227 break;
228 }
229 case OPERATOR_DIV: {
230 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
231 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
232 break;
233 }
234 case OPERATOR_D2: {
235 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
236 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
237 break;
238 }
239 case OPERATOR_D3:
240 case OPERATOR_D4:
241 case OPERATOR_D5:
242 case OPERATOR_D6:
243 case OPERATOR_D7:
244 case OPERATOR_D8:
245 case OPERATOR_D9:
246 case OPERATOR_D10: {
247 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
248 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
249 break;
250 }
251 default: {
252 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
253 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
254 }
255 }
256 }
257 }
258 // -------------------------------------------------------------------------------------
259
260 template<typename DT, typename OT, typename PT>
261 Basis_HGRAD_TET_C2_FEM<DT,OT,PT>::
262 Basis_HGRAD_TET_C2_FEM() {
263 const ordinal_type spaceDim = 3;
264 this->basisCardinality_ = 10;
265 this->basisDegree_ = 2;
266 this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
267 this->basisType_ = BASIS_FEM_DEFAULT;
268 this->basisCoordinates_ = COORDINATES_CARTESIAN;
269 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
270
271 // intialize tags
272 {
273 // Basis-dependent intializations
274 const ordinal_type tagSize = 4; // size of DoF tag
275 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
276 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
277 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
278
279 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
280 ordinal_type tags[40] = { 0, 0, 0, 1,
281 0, 1, 0, 1,
282 0, 2, 0, 1,
283 0, 3, 0, 1,
284 1, 0, 0, 1,
285 1, 1, 0, 1,
286 1, 2, 0, 1,
287 1, 3, 0, 1,
288 1, 4, 0, 1,
289 1, 5, 0, 1,
290 };
291
292 // host tags
293 OrdinalTypeArray1DHost tagView(&tags[0], 40);
294
295 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
296 this->setOrdinalTagData(this->tagToOrdinal_,
297 this->ordinalToTag_,
298 tagView,
299 this->basisCardinality_,
300 tagSize,
301 posScDim,
302 posScOrd,
303 posDfOrd);
304 }
305
306 // dofCoords on host and create its mirror view to device
307 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
308 dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
309
310 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
311 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
312 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
313 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
314
315 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
316 dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
317 dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
318 dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
319 dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
320 dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
321
322 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
323 Kokkos::deep_copy(this->dofCoords_, dofCoords);
324 }
325
326 template<typename DT, typename OT, typename PT>
327 void
328 Basis_HGRAD_TET_C2_FEM<DT,OT,PT>::getScratchSpaceSize(
329 ordinal_type& perTeamSpaceSize,
330 ordinal_type& perThreadSpaceSize,
331 const PointViewType inputPoints,
332 const EOperator operatorType) const {
333 perTeamSpaceSize = 0;
334 perThreadSpaceSize = 0;
335 }
336
337 template<typename DT, typename OT, typename PT>
338 KOKKOS_INLINE_FUNCTION
339 void
340 Basis_HGRAD_TET_C2_FEM<DT,OT,PT>::getValues(
341 OutputViewType outputValues,
342 const PointViewType inputPoints,
343 const EOperator operatorType,
344 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
345 const typename DT::execution_space::scratch_memory_space & scratchStorage,
346 const ordinal_type subcellDim,
347 const ordinal_type subcellOrdinal) const {
348
349 INTREPID2_TEST_FOR_ABORT( !((subcellDim <= 0) && (subcellOrdinal == -1)),
350 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_C2_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
351
352 (void) scratchStorage; //avoid unused variable warning
353
354 const int numPoints = inputPoints.extent(0);
355
356 switch(operatorType) {
357 case OPERATOR_VALUE:
358 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
359 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
360 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
361 Impl::Basis_HGRAD_TET_C2_FEM::Serial<OPERATOR_VALUE>::getValues( output, input);
362 });
363 break;
364 case OPERATOR_GRAD:
365 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
366 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
367 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
368 Impl::Basis_HGRAD_TET_C2_FEM::Serial<OPERATOR_GRAD>::getValues( output, input);
369 });
370 break;
371 default: {}
372 }
373 }
374
375}// namespace Intrepid2
376#endif
Intrepid2::Basis_HGRAD_TET_C2_FEM::Basis_HGRAD_TET_C2_FEM
Basis_HGRAD_TET_C2_FEM()
Constructor.
Definition Intrepid2_HGRAD_TET_C2_FEMDef.hpp:262
Intrepid2::Basis_HGRAD_TET_C2_FEM::getScratchSpaceSize
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...
Definition Intrepid2_HGRAD_TET_C2_FEMDef.hpp:328
Intrepid2::Basis_HGRAD_TET_C2_FEM::getValues
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.
Definition Intrepid2_HGRAD_TET_C2_FEM.hpp:167
Intrepid2::Basis::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition Intrepid2_Basis.hpp:338
Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125
Intrepid2::Impl::Basis_HGRAD_TET_C2_FEM::Serial
See Intrepid2::Basis_HGRAD_TET_C2_FEM.
Definition Intrepid2_HGRAD_TET_C2_FEM.hpp:80
Generated by doxygen 1.9.8