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