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