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 
19 namespace 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::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition: Intrepid2_Basis.hpp:338
Intrepid2::Basis_HGRAD_TRI_C2_FEM::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
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
Generated by   doxygen 1.8.5