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