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