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