Intrepid2
Intrepid2_HCURL_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_HCURL_QUAD_I1_FEM_DEF_HPP__
17 #define __INTREPID2_HCURL_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_HCURL_QUAD_I1_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  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
39  output.access(0, 0) = 0.5*(1.0 - y);
40  output.access(0, 1) = 0.0;
41 
42  output.access(1, 0) = 0.0;
43  output.access(1, 1) = 0.5*(1.0 + x);
44 
45  output.access(2, 0) = -0.5*(1.0 + y);
46  output.access(2, 1) = 0.0;
47 
48  output.access(3, 0) = 0.0;
49  output.access(3, 1) = -0.5*(1.0 - x);
50  break;
51  }
52  case OPERATOR_CURL: {
53  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
54  output.access(0) = 0.5;
55  output.access(1) = 0.5;
56  output.access(2) = 0.5;
57  output.access(3) = 0.5;
58  break;
59  }
60  default: {
61  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
62  opType != OPERATOR_CURL,
63  ">>> ERROR: (Intrepid2::Basis_HCURL_QUAD_I1_FEM::Serial::getValues) operator is not supported");
64  }
65  } //end switch
66  }
67 
68  template<typename DT,
69  typename outputValueValueType, class ...outputValueProperties,
70  typename inputPointValueType, class ...inputPointProperties>
71  void
72  Basis_HCURL_QUAD_I1_FEM::
73  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
74  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
75  const EOperator operatorType ) {
76  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
77  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
78  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
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_CURL: {
91  typedef Functor<outputValueViewType, inputPointViewType, OPERATOR_CURL> FunctorType;
92  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
93  break;
94  }
95  case OPERATOR_DIV: {
96  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV, std::invalid_argument,
97  ">>> ERROR (Basis_HCURL_QUAD_I1_FEM::getValues): DIV is invalid operator for HCURL Basis Functions");
98  break;
99  }
100  case OPERATOR_GRAD: {
101  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_GRAD, std::invalid_argument,
102  ">>> ERROR (Basis_HCURL_QUAD_I1_FEM::getValues): GRAD is invalid operator for HCURL 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_HCURL_QUAD_I1_FEM::getValues): 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_HCURL_QUAD_I1_FEM::getValues): Invalid operator type");
146  }
147  }
148  }
149  }
150 
151  // -------------------------------------------------------------------------------------
152 
153 
154  template<typename DT, typename OT, typename PT>
157  const ordinal_type spaceDim = 2;
158  this->basisCardinality_ = 4;
159  this->basisDegree_ = 1;
160  this->basisCellTopologyKey_ = shards::Quadrilateral<4>::key;
161  this->basisType_ = BASIS_FEM_DEFAULT;
162  this->basisCoordinates_ = COORDINATES_CARTESIAN;
163  this->functionSpace_ = FUNCTION_SPACE_HCURL;
164 
165  // initialize tags
166  {
167  // Basis-dependent intializations
168  const ordinal_type tagSize = 4; // size of DoF tag
169  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
170  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
171  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
172 
173  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
174  ordinal_type tags[16] = { 1, 0, 0, 1,
175  1, 1, 0, 1,
176  1, 2, 0, 1,
177  1, 3, 0, 1 };
178 
179  // when exec space is device, this wrapping relies on uvm.
180  //Kokkos::View<ordinal_type[16], DT> tagView(tags);
181  OrdinalTypeArray1DHost tagView(&tags[0], 16);
182 
183  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
184  this->setOrdinalTagData(this->tagToOrdinal_,
185  this->ordinalToTag_,
186  tagView,
187  this->basisCardinality_,
188  tagSize,
189  posScDim,
190  posScOrd,
191  posDfOrd);
192  }
193 
194  // dofCoords on host and create its mirror view to device
195  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
196  dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
197 
198  dofCoords(0,0) = 0.0; dofCoords(0,1) = -1.0;
199  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0;
200  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0;
201  dofCoords(3,0) = -1.0; dofCoords(3,1) = 0.0;
202 
203  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
204  Kokkos::deep_copy(this->dofCoords_, dofCoords);
205 
206 
207  // dofCoeffs on host and create its mirror view to device
208  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
209  dofCoeffs("dofCoeffsHost", this->basisCardinality_,spaceDim);
210 
211  // for HCURL_QUAD_I1 dofCoeffs are the tangents on the quadrilateral edges (with tangents magnitude equal to edges' lengths)
212  dofCoeffs(0,0) = 1.0; dofCoeffs(0,1) = 0.0;
213  dofCoeffs(1,0) = 0.0; dofCoeffs(1,1) = 1.0;
214  dofCoeffs(2,0) = -1.0; dofCoeffs(2,1) = 0.0;
215  dofCoeffs(3,0) = 0.0; dofCoeffs(3,1) = -1.0;
216 
217  this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffs);
218  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
219 
220  }
221 
222  template<typename DT, typename OT, typename PT>
223  void
225  ordinal_type& perTeamSpaceSize,
226  ordinal_type& perThreadSpaceSize,
227  const PointViewType inputPoints,
228  const EOperator operatorType) const {
229  perTeamSpaceSize = 0;
230  perThreadSpaceSize = 0;
231  }
232 
233  template<typename DT, typename OT, typename PT>
234  KOKKOS_INLINE_FUNCTION
235  void
236  Basis_HCURL_QUAD_I1_FEM<DT,OT,PT>::getValues(
237  OutputViewType outputValues,
238  const PointViewType inputPoints,
239  const EOperator operatorType,
240  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
241  const typename DT::execution_space::scratch_memory_space & scratchStorage,
242  const ordinal_type subcellDim,
243  const ordinal_type subcellOrdinal) const {
244 
245  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
246  ">>> ERROR: (Intrepid2::Basis_HCURL_QUAD_I1_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
247 
248  (void) scratchStorage; //avoid unused variable warning
249 
250  const int numPoints = inputPoints.extent(0);
251 
252  switch(operatorType) {
253  case OPERATOR_VALUE:
254  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
255  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
256  const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
257  Impl::Basis_HCURL_QUAD_I1_FEM::Serial<OPERATOR_VALUE>::getValues( output, input);
258  });
259  break;
260  case OPERATOR_CURL:
261  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
262  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), pt, Kokkos::ALL() );
263  const auto input = Kokkos::subview( inputPoints, pt, Kokkos::ALL() );
264  Impl::Basis_HCURL_QUAD_I1_FEM::Serial<OPERATOR_CURL>::getValues( output, input);
265  });
266  break;
267  default: {
268  INTREPID2_TEST_FOR_ABORT( true, ">>> ERROR: (Intrepid2::Basis_HCURL_QUAD_I1_FEM::getValues), Operator Type not supported.");
269  }
270  }
271  }
272 
273 }// namespace Intrepid2
274 #endif
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Quadrilateral cell...