Intrepid2
Intrepid2_HGRAD_WEDGE_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_WEDGE_C1_FEM_DEF_HPP__
17 #define __INTREPID2_HGRAD_WEDGE_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_WEDGE_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  const auto z = input(2);
38 
39  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
40  output.access(0) = (1.0 - x - y)*(1.0 - z)/2.0;
41  output.access(1) = x*(1.0 - z)/2.0;
42  output.access(2) = y*(1.0 - z)/2.0;
43  output.access(3) = (1.0 - x - y)*(1.0 + z)/2.0;
44  output.access(4) = x*(1.0 + z)/2.0;
45  output.access(5) = y*(1.0 + z)/2.0;
46  break;
47  }
48  case OPERATOR_GRAD: {
49  const auto x = input(0);
50  const auto y = input(1);
51  const auto z = input(2);
52 
53  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
54  output.access(0, 0) = -(1.0 - z)/2.0;
55  output.access(0, 1) = -(1.0 - z)/2.0;
56  output.access(0, 2) = -(1.0 - x - y)/2.0;
57 
58  output.access(1, 0) = (1.0 - z)/2.0;
59  output.access(1, 1) = 0.0;
60  output.access(1, 2) = -x/2.0;
61 
62  output.access(2, 0) = 0.0;
63  output.access(2, 1) = (1.0 - z)/2.0;
64  output.access(2, 2) = -y/2.0;
65 
66 
67  output.access(3, 0) = -(1.0 + z)/2.0;
68  output.access(3, 1) = -(1.0 + z)/2.0;
69  output.access(3, 2) = (1.0 - x - y)/2.0;
70 
71  output.access(4, 0) = (1.0 + z)/2.0;
72  output.access(4, 1) = 0.0;
73  output.access(4, 2) = x/2.0;
74 
75  output.access(5, 0) = 0.0;
76  output.access(5, 1) = (1.0 + z)/2.0;
77  output.access(5, 2) = y/2.0;
78  break;
79  }
80  case OPERATOR_D2: {
81  output.access(0, 0) = 0.0; output.access(3, 0) = 0.0;
82  output.access(0, 1) = 0.0; output.access(3, 1) = 0.0;
83  output.access(0, 2) = 0.5; output.access(3, 2) =-0.5;
84  output.access(0, 3) = 0.0; output.access(3, 3) = 0.0;
85  output.access(0, 4) = 0.5; output.access(3, 4) =-0.5;
86  output.access(0, 5) = 0.0; output.access(3, 5) = 0.0;
87 
88  output.access(1, 0) = 0.0; output.access(4, 0) = 0.0;
89  output.access(1, 1) = 0.0; output.access(4, 1) = 0.0;
90  output.access(1, 2) =-0.5; output.access(4, 2) = 0.5;
91  output.access(1, 3) = 0.0; output.access(4, 3) = 0.0;
92  output.access(1, 4) = 0.0; output.access(4, 4) = 0.0;
93  output.access(1, 5) = 0.0; output.access(4, 5) = 0.0;
94 
95  output.access(2, 0) = 0.0; output.access(5, 0) = 0.0;
96  output.access(2, 1) = 0.0; output.access(5, 1) = 0.0;
97  output.access(2, 2) = 0.0; output.access(5, 2) = 0.0;
98  output.access(2, 3) = 0.0; output.access(5, 3) = 0.0;
99  output.access(2, 4) =-0.5; output.access(5, 4) = 0.5;
100  output.access(2, 5) = 0.0; output.access(5, 5) = 0.0;
101  break;
102  }
103  case OPERATOR_MAX : {
104  const ordinal_type jend = output.extent(1);
105  const ordinal_type iend = output.extent(0);
106 
107  for (ordinal_type j=0;j<jend;++j)
108  for (ordinal_type i=0;i<iend;++i)
109  output.access(i, j) = 0.0;
110  break;
111  }
112  default: {
113  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
114  opType != OPERATOR_GRAD &&
115  opType != OPERATOR_D2 &&
116  opType != OPERATOR_MAX,
117  ">>> ERROR: (Intrepid2::Basis_HGRAD_WEDGE_C1_FEM::Serial::getValues) operator is not supported");
118 
119  }
120  }
121  }
122 
123  template<typename DT,
124  typename outputValueValueType, class ...outputValueProperties,
125  typename inputPointValueType, class ...inputPointProperties>
126  void
127  Basis_HGRAD_WEDGE_C1_FEM::
128  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
129  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
130  const EOperator operatorType ) {
131  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
132  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
133  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
134 
135  // Number of evaluation points = dim 0 of inputPoints
136  const auto loopSize = inputPoints.extent(0);
137  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
138 
139  switch (operatorType) {
140 
141  case OPERATOR_VALUE: {
142  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
143  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
144  break;
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  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
154  ">>> ERROR (Basis_HGRAD_WEDGE_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
155  break;
156  }
157 
158  case OPERATOR_DIV: {
159  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
160  ">>> ERROR (Basis_HGRAD_WEDGE_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
161  break;
162  }
163 
164  case OPERATOR_D2: {
165  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
166  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
167  break;
168  }
169  case OPERATOR_D3:
170  case OPERATOR_D4:
171  case OPERATOR_D5:
172  case OPERATOR_D6:
173  case OPERATOR_D7:
174  case OPERATOR_D8:
175  case OPERATOR_D9:
176  case OPERATOR_D10: {
177  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
178  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
179  break;
180  }
181  default: {
182  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
183  ">>> ERROR (Basis_HGRAD_WEDGE_C1_FEM): Invalid operator type");
184  }
185  }
186  }
187  }
188  // -------------------------------------------------------------------------------------
189 
190  template<typename DT, typename OT, typename PT>
191  Basis_HGRAD_WEDGE_C1_FEM<DT,OT,PT>::
192  Basis_HGRAD_WEDGE_C1_FEM() {
193  const ordinal_type spaceDim = 3;
194  this->basisCardinality_ = 6;
195  this->basisDegree_ = 1;
196  this->basisCellTopologyKey_ = shards::Wedge<6>::key;
197  this->basisType_ = BASIS_FEM_DEFAULT;
198  this->basisCoordinates_ = COORDINATES_CARTESIAN;
199  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
200 
201  // initialize tags
202  {
203  // Basis-dependent intializations
204  const ordinal_type tagSize = 4; // size of DoF tag
205  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
206  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
207  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
208 
209  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
210  ordinal_type tags[24] = { 0, 0, 0, 1,
211  0, 1, 0, 1,
212  0, 2, 0, 1,
213  0, 3, 0, 1,
214  0, 4, 0, 1,
215  0, 5, 0, 1 };
216 
217  // host tags
218  OrdinalTypeArray1DHost tagView(&tags[0], 24);
219 
220  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
221  //OrdinalTypeArray2DHost ordinalToTag;
222  //OrdinalTypeArray3DHost tagToOrdinal;
223  this->setOrdinalTagData(this->tagToOrdinal_,
224  this->ordinalToTag_,
225  tagView,
226  this->basisCardinality_,
227  tagSize,
228  posScDim,
229  posScOrd,
230  posDfOrd);
231  }
232 
233  // dofCoords on host and create its mirror view to device
234  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
235  dofCoords("dofCoordsHost", this->basisCardinality_,spaceDim);
236 
237  dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = -1.0;
238  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = -1.0;
239  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = -1.0;
240  dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
241  dofCoords(4,0) = 1.0; dofCoords(4,1) = 0.0; dofCoords(4,2) = 1.0;
242  dofCoords(5,0) = 0.0; dofCoords(5,1) = 1.0; dofCoords(5,2) = 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
250  Basis_HGRAD_WEDGE_C1_FEM<DT,OT,PT>::getScratchSpaceSize(
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
262  Basis_HGRAD_WEDGE_C1_FEM<DT,OT,PT>::getValues(
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_WEDGE_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() );
283  Impl::Basis_HGRAD_WEDGE_C1_FEM::template Serial<OPERATOR_VALUE>::getValues( output, input);
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_WEDGE_C1_FEM::template Serial<OPERATOR_GRAD>::getValues( output, input);
291  });
292  break;
293  default: {}
294  }
295  }
296 
297 }// namespace Intrepid2
298 #endif
Intrepid2::Basis_HGRAD_WEDGE_C1_FEM
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Wedge cell.
Definition: Intrepid2_HGRAD_WEDGE_C1_FEM.hpp:134
Intrepid2::Basis_HGRAD_WEDGE_C1_FEM::Basis_HGRAD_WEDGE_C1_FEM
Basis_HGRAD_WEDGE_C1_FEM()
Constructor.
Definition: Intrepid2_HGRAD_WEDGE_C1_FEMDef.hpp:192
Generated by   doxygen 1.8.5