Intrepid2
Intrepid2_HGRAD_QUAD_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_QUAD_C1_FEM_DEF_HPP__
17#define __INTREPID2_HGRAD_QUAD_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_QUAD_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
38 // output with dimensions (basisCardinality_)
39 output.access(0) = (1.0 - x)*(1.0 - y)/4.0;
40 output.access(1) = (1.0 + x)*(1.0 - y)/4.0;
41 output.access(2) = (1.0 + x)*(1.0 + y)/4.0;
42 output.access(3) = (1.0 - x)*(1.0 + y)/4.0;
43 break;
44 }
45 case OPERATOR_GRAD : {
46 const auto x = input(0);
47 const auto y = input(1);
48
49 // output with dimensions (basisCardinality_, spaceDim)
50 output.access(0, 0) = -(1.0 - y)/4.0;
51 output.access(0, 1) = -(1.0 - x)/4.0;
52
53 output.access(1, 0) = (1.0 - y)/4.0;
54 output.access(1, 1) = -(1.0 + x)/4.0;
55
56 output.access(2, 0) = (1.0 + y)/4.0;
57 output.access(2, 1) = (1.0 + x)/4.0;
58
59 output.access(3, 0) = -(1.0 + y)/4.0;
60 output.access(3, 1) = (1.0 - x)/4.0;
61 break;
62 }
63 case OPERATOR_CURL : {
64 const auto x = input(0);
65 const auto y = input(1);
66
67 // output with dimensions (basisCardinality_, spaceDim)
68 output.access(0, 0) = -(1.0 - x)/4.0;
69 output.access(0, 1) = (1.0 - y)/4.0;
70
71 output.access(1, 0) = -(1.0 + x)/4.0;
72 output.access(1, 1) = -(1.0 - y)/4.0;
73
74 output.access(2, 0) = (1.0 + x)/4.0;
75 output.access(2, 1) = -(1.0 + y)/4.0;
76
77 output.access(3, 0) = (1.0 - x)/4.0;
78 output.access(3, 1) = (1.0 + y)/4.0;
79 break;
80 }
81 case OPERATOR_D2 : {
82 // output with dimensions (basisCardinality_, D2Cardinality=3)
83 output.access(0, 0) = 0.0;
84 output.access(0, 1) = 0.25;
85 output.access(0, 2) = 0.0;
86
87 output.access(1, 0) = 0.0;
88 output.access(1, 1) = -0.25;
89 output.access(1, 2) = 0.0;
90
91 output.access(2, 0) = 0.0;
92 output.access(2, 1) = 0.25;
93 output.access(2, 2) = 0.0;
94
95 output.access(3, 0) = 0.0;
96 output.access(3, 1) = -0.25;
97 output.access(3, 2) = 0.0;
98 break;
99 }
100 case OPERATOR_MAX : {
101 const ordinal_type jend = output.extent(1);
102 const ordinal_type iend = output.extent(0);
103
104 for (ordinal_type j=0;j<jend;++j)
105 for (ordinal_type i=0;i<iend;++i)
106 output.access(i, j) = 0.0;
107 break;
108 }
109 default: {
110 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
111 opType != OPERATOR_GRAD &&
112 opType != OPERATOR_CURL &&
113 opType != OPERATOR_D2 &&
114 opType != OPERATOR_MAX,
115 ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C1_FEM::Serial::getValues) operator is not supported");
116
117 }
118 }
119 }
120
121 template<typename DT,
122 typename outputValueValueType, class ...outputValueProperties,
123 typename inputPointValueType, class ...inputPointProperties>
124 void
125 Basis_HGRAD_QUAD_C1_FEM::
126 getValues( const typename DT::execution_space& space,
127 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
128 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
129 const EOperator operatorType ) {
130 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
131 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
132 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
133
134 // Number of evaluation points = dim 0 of inputPoints
135 const auto loopSize = inputPoints.extent(0);
136 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
137
138 switch (operatorType) {
139
140 case OPERATOR_VALUE: {
141 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
142 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
143 break;
144 }
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 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
154 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
155 break;
156 }
157 case OPERATOR_DIV: {
158 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV, std::invalid_argument,
159 ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
160 break;
161 }
162 case OPERATOR_D2: {
163 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
164 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
165 break;
166 }
167 case OPERATOR_D3:
168 case OPERATOR_D4:
169 case OPERATOR_D5:
170 case OPERATOR_D6:
171 case OPERATOR_D7:
172 case OPERATOR_D8:
173 case OPERATOR_D9:
174 case OPERATOR_D10: {
175 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
176 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
177 break;
178 }
179 default: {
180 INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
181 ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): Invalid operator type");
182 }
183 }
184 }
185 }
186 // -------------------------------------------------------------------------------------
187
188 template<typename DT, typename OT, typename PT>
191 const ordinal_type spaceDim = 2;
192 this->basisCardinality_ = 4;
193 this->basisDegree_ = 1;
194 this->basisCellTopologyKey_ = shards::Quadrilateral<4>::key;
195 this->basisType_ = BASIS_FEM_DEFAULT;
196 this->basisCoordinates_ = COORDINATES_CARTESIAN;
197 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
198
199 // initialize tags
200 {
201 // Basis-dependent intializations
202 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
203 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
204 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
205 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
206
207 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
208 ordinal_type tags[16] = { 0, 0, 0, 1,
209 0, 1, 0, 1,
210 0, 2, 0, 1,
211 0, 3, 0, 1 };
212
213 // host tags
214 OrdinalTypeArray1DHost tagView(&tags[0], 16);
215
216 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
217 //OrdinalTypeArray2DHost ordinalToTag;
218 //OrdinalTypeArray3DHost tagToOrdinal;
219 this->setOrdinalTagData(this->tagToOrdinal_,
220 this->ordinalToTag_,
221 tagView,
222 this->basisCardinality_,
223 tagSize,
224 posScDim,
225 posScOrd,
226 posDfOrd);
227
228 //this->tagToOrdinal_ = Kokkos::create_mirror_view(typename DT::memory_space(), tagToOrdinal);
229 //Kokkos::deep_copy(this->tagToOrdinal_, tagToOrdinal);
230
231 //this->ordinalToTag_ = Kokkos::create_mirror_view(typename DT::memory_space(), ordinalToTag);
232 //Kokkos::deep_copy(this->ordinalToTag_, ordinalToTag);
233 }
234
235 // dofCoords on host and create its mirror view to device
236 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
237 dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
238
239 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0;
240 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0;
241 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0;
242 dofCoords(3,0) = -1.0; dofCoords(3,1) = 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
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
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_QUAD_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() );
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_QUAD_C1_FEM::Serial<OPERATOR_GRAD>::getValues( output, input);
291 });
292 break;
293 case OPERATOR_CURL:
294 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
295 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
296 const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
297 Impl::Basis_HGRAD_QUAD_C1_FEM::Serial<OPERATOR_CURL>::getValues( output, input);
298 });
299 break;
300 default: {
301 INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C1_FEM::getValues), Operator Type not supported.");
302 }
303 }
304 }
305
306}// namespace Intrepid2
307#endif
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...
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.
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.