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 
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_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::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.