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