Intrepid2
Intrepid2_HGRAD_HEX_C1_FEMDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
49 #ifndef __INTREPID2_HGRAD_HEX_C1_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_HEX_C1_FEM_DEF_HPP__
51 
52 namespace Intrepid2 {
53 
54  // -------------------------------------------------------------------------------------
55  namespace Impl {
56 
57  template<EOperator opType>
58  template<typename OutputViewType,
59  typename inputViewType>
60  KOKKOS_INLINE_FUNCTION
61  void
62  Basis_HGRAD_HEX_C1_FEM::Serial<opType>::
63  getValues( OutputViewType output,
64  const inputViewType input ) {
65  switch (opType) {
66  case OPERATOR_VALUE : {
67  const auto x = input(0);
68  const auto y = input(1);
69  const auto z = input(2);
70 
71  // output is a rank-2 array with dimensions (basisCardinality_, dim0)
72  output.access(0) = (1.0 - x)*(1.0 - y)*(1.0 - z)/8.0;
73  output.access(1) = (1.0 + x)*(1.0 - y)*(1.0 - z)/8.0;
74  output.access(2) = (1.0 + x)*(1.0 + y)*(1.0 - z)/8.0;
75  output.access(3) = (1.0 - x)*(1.0 + y)*(1.0 - z)/8.0;
76 
77  output.access(4) = (1.0 - x)*(1.0 - y)*(1.0 + z)/8.0;
78  output.access(5) = (1.0 + x)*(1.0 - y)*(1.0 + z)/8.0;
79  output.access(6) = (1.0 + x)*(1.0 + y)*(1.0 + z)/8.0;
80  output.access(7) = (1.0 - x)*(1.0 + y)*(1.0 + z)/8.0;
81  break;
82  }
83  case OPERATOR_GRAD : {
84  const auto x = input(0);
85  const auto y = input(1);
86  const auto z = input(2);
87 
88  // output is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
89  output.access(0, 0) = -(1.0 - y)*(1.0 - z)/8.0;
90  output.access(0, 1) = -(1.0 - x)*(1.0 - z)/8.0;
91  output.access(0, 2) = -(1.0 - x)*(1.0 - y)/8.0;
92 
93  output.access(1, 0) = (1.0 - y)*(1.0 - z)/8.0;
94  output.access(1, 1) = -(1.0 + x)*(1.0 - z)/8.0;
95  output.access(1, 2) = -(1.0 + x)*(1.0 - y)/8.0;
96 
97  output.access(2, 0) = (1.0 + y)*(1.0 - z)/8.0;
98  output.access(2, 1) = (1.0 + x)*(1.0 - z)/8.0;
99  output.access(2, 2) = -(1.0 + x)*(1.0 + y)/8.0;
100 
101  output.access(3, 0) = -(1.0 + y)*(1.0 - z)/8.0;
102  output.access(3, 1) = (1.0 - x)*(1.0 - z)/8.0;
103  output.access(3, 2) = -(1.0 - x)*(1.0 + y)/8.0;
104 
105  output.access(4, 0) = -(1.0 - y)*(1.0 + z)/8.0;
106  output.access(4, 1) = -(1.0 - x)*(1.0 + z)/8.0;
107  output.access(4, 2) = (1.0 - x)*(1.0 - y)/8.0;
108 
109  output.access(5, 0) = (1.0 - y)*(1.0 + z)/8.0;
110  output.access(5, 1) = -(1.0 + x)*(1.0 + z)/8.0;
111  output.access(5, 2) = (1.0 + x)*(1.0 - y)/8.0;
112 
113  output.access(6, 0) = (1.0 + y)*(1.0 + z)/8.0;
114  output.access(6, 1) = (1.0 + x)*(1.0 + z)/8.0;
115  output.access(6, 2) = (1.0 + x)*(1.0 + y)/8.0;
116 
117  output.access(7, 0) = -(1.0 + y)*(1.0 + z)/8.0;
118  output.access(7, 1) = (1.0 - x)*(1.0 + z)/8.0;
119  output.access(7, 2) = (1.0 - x)*(1.0 + y)/8.0;
120  break;
121  }
122  case OPERATOR_D2 : {
123  const auto x = input(0);
124  const auto y = input(1);
125  const auto z = input(2);
126 
127  // output is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6)
128  output.access(0, 0) = 0.0; // {2, 0, 0}
129  output.access(0, 1) = (1.0 - z)/8.0; // {1, 1, 0}
130  output.access(0, 2) = (1.0 - y)/8.0; // {1, 0, 1}
131  output.access(0, 3) = 0.0; // {0, 2, 0}
132  output.access(0, 4) = (1.0 - x)/8.0; // {0, 1, 1}
133  output.access(0, 5) = 0.0; // {0, 0, 2}
134 
135  output.access(1, 0) = 0.0; // {2, 0, 0}
136  output.access(1, 1) = -(1.0 - z)/8.0; // {1, 1, 0}
137  output.access(1, 2) = -(1.0 - y)/8.0; // {1, 0, 1}
138  output.access(1, 3) = 0.0; // {0, 2, 0}
139  output.access(1, 4) = (1.0 + x)/8.0; // {0, 1, 1}
140  output.access(1, 5) = 0.0; // {0, 0, 2}
141 
142  output.access(2, 0) = 0.0; // {2, 0, 0}
143  output.access(2, 1) = (1.0 - z)/8.0; // {1, 1, 0}
144  output.access(2, 2) = -(1.0 + y)/8.0; // {1, 0, 1}
145  output.access(2, 3) = 0.0; // {0, 2, 0}
146  output.access(2, 4) = -(1.0 + x)/8.0; // {0, 1, 1}
147  output.access(2, 5) = 0.0; // {0, 0, 2}
148 
149  output.access(3, 0) = 0.0; // {2, 0, 0}
150  output.access(3, 1) = -(1.0 - z)/8.0; // {1, 1, 0}
151  output.access(3, 2) = (1.0 + y)/8.0; // {1, 0, 1}
152  output.access(3, 3) = 0.0; // {0, 2, 0}
153  output.access(3, 4) = -(1.0 - x)/8.0; // {0, 1, 1}
154  output.access(3, 5) = 0.0; // {0, 0, 2}
155 
156 
157  output.access(4, 0) = 0.0; // {2, 0, 0}
158  output.access(4, 1) = (1.0 + z)/8.0; // {1, 1, 0}
159  output.access(4, 2) = -(1.0 - y)/8.0; // {1, 0, 1}
160  output.access(4, 3) = 0.0; // {0, 2, 0}
161  output.access(4, 4) = -(1.0 - x)/8.0; // {0, 1, 1}
162  output.access(4, 5) = 0.0; // {0, 0, 2}
163 
164  output.access(5, 0) = 0.0; // {2, 0, 0}
165  output.access(5, 1) = -(1.0 + z)/8.0; // {1, 1, 0}
166  output.access(5, 2) = (1.0 - y)/8.0; // {1, 0, 1}
167  output.access(5, 3) = 0.0; // {0, 2, 0}
168  output.access(5, 4) = -(1.0 + x)/8.0; // {0, 1, 1}
169  output.access(5, 5) = 0.0; // {0, 0, 2}
170 
171  output.access(6, 0) = 0.0; // {2, 0, 0}
172  output.access(6, 1) = (1.0 + z)/8.0; // {1, 1, 0}
173  output.access(6, 2) = (1.0 + y)/8.0; // {1, 0, 1}
174  output.access(6, 3) = 0.0; // {0, 2, 0}
175  output.access(6, 4) = (1.0 + x)/8.0; // {0, 1, 1}
176  output.access(6, 5) = 0.0; // {0, 0, 2}
177 
178  output.access(7, 0) = 0.0; // {2, 0, 0}
179  output.access(7, 1) = -(1.0 + z)/8.0; // {1, 1, 0}
180  output.access(7, 2) = -(1.0 + y)/8.0; // {1, 0, 1}
181  output.access(7, 3) = 0.0; // {0, 2, 0}
182  output.access(7, 4) = (1.0 - x)/8.0; // {0, 1, 1}
183  output.access(7, 5) = 0.0; // {0, 0, 2}
184  break;
185  }
186  case OPERATOR_MAX : {
187  const ordinal_type jend = output.extent(1);
188  const ordinal_type iend = output.extent(0);
189 
190  for (ordinal_type j=0;j<jend;++j)
191  for (ordinal_type i=0;i<iend;++i)
192  output.access(i, j) = 0.0;
193  break;
194  }
195  default: {
196  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
197  opType != OPERATOR_GRAD &&
198  opType != OPERATOR_CURL &&
199  opType != OPERATOR_D2 &&
200  opType != OPERATOR_MAX,
201  ">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C1_FEM::Serial::getValues) operator is not supported");
202 
203  }
204  }
205  }
206 
207  template<typename SpT,
208  typename outputValueValueType, class ...outputValueProperties,
209  typename inputPointValueType, class ...inputPointProperties>
210  void
211  Basis_HGRAD_HEX_C1_FEM::
212  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
213  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
214  const EOperator operatorType ) {
215  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
216  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
217  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
218 
219  // Number of evaluation points = dim 0 of inputPoints
220  const auto loopSize = inputPoints.extent(0);
221  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
222 
223  switch (operatorType) {
224 
225  case OPERATOR_VALUE: {
226  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
227  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
228  break;
229  }
230  case OPERATOR_GRAD:
231  case OPERATOR_D1: {
232  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
233  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
234  break;
235  }
236  case OPERATOR_CURL: {
237  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
238  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
239  break;
240  }
241 
242  case OPERATOR_DIV: {
243  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
244  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
245  break;
246  }
247 
248  case OPERATOR_D2: {
249  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
250  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
251  break;
252  }
253  case OPERATOR_D3:
254  case OPERATOR_D4:
255  case OPERATOR_D5:
256  case OPERATOR_D6:
257  case OPERATOR_D7:
258  case OPERATOR_D8:
259  case OPERATOR_D9:
260  case OPERATOR_D10: {
261  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
262  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
263  break;
264  }
265  default: {
266  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
267  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): Invalid operator type");
268  }
269  }
270  }
271  }
272 
273  // -------------------------------------------------------------------------------------
274 
275  template<typename SpT, typename OT, typename PT>
276  Basis_HGRAD_HEX_C1_FEM<SpT,OT,PT>::
277  Basis_HGRAD_HEX_C1_FEM() {
278  this->basisCardinality_ = 8;
279  this->basisDegree_ = 1;
280  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
281  this->basisType_ = BASIS_FEM_DEFAULT;
282  this->basisCoordinates_ = COORDINATES_CARTESIAN;
283  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
284 
285  // initialize tags
286  {
287  // Basis-dependent intializations
288  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
289  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
290  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
291  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
292 
293  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
294  ordinal_type tags[32] = { 0, 0, 0, 1,
295  0, 1, 0, 1,
296  0, 2, 0, 1,
297  0, 3, 0, 1,
298  0, 4, 0, 1,
299  0, 5, 0, 1,
300  0, 6, 0, 1,
301  0, 7, 0, 1 };
302 
303  // host tags
304  OrdinalTypeArray1DHost tagView(&tags[0], 32);
305 
306  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
307  //OrdinalTypeArray2DHost ordinalToTag;
308  //OrdinalTypeArray3DHost tagToOrdinal;
309  this->setOrdinalTagData(this->tagToOrdinal_,
310  this->ordinalToTag_,
311  tagView,
312  this->basisCardinality_,
313  tagSize,
314  posScDim,
315  posScOrd,
316  posDfOrd);
317 
318  //this->tagToOrdinal_ = Kokkos::create_mirror_view(typename SpT::memory_space(), tagToOrdinal);
319  //Kokkos::deep_copy(this->tagToOrdinal_, tagToOrdinal);
320 
321  //this->ordinalToTag_ = Kokkos::create_mirror_view(typename SpT::memory_space(), ordinalToTag);
322  //Kokkos::deep_copy(this->ordinalToTag_, ordinalToTag);
323  }
324 
325  // dofCoords on host and create its mirror view to device
326  Kokkos::DynRankView<typename ScalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
327  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
328 
329  dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = -1.0;
330  dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = -1.0;
331  dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = -1.0;
332  dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = -1.0;
333  dofCoords(4,0) = -1.0; dofCoords(4,1) = -1.0; dofCoords(4,2) = 1.0;
334  dofCoords(5,0) = 1.0; dofCoords(5,1) = -1.0; dofCoords(5,2) = 1.0;
335  dofCoords(6,0) = 1.0; dofCoords(6,1) = 1.0; dofCoords(6,2) = 1.0;
336  dofCoords(7,0) = -1.0; dofCoords(7,1) = 1.0; dofCoords(7,2) = 1.0;
337 
338  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
339  Kokkos::deep_copy(this->dofCoords_, dofCoords);
340  }
341 
342 }// namespace Intrepid2
343 
344 #endif
345 
Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecSpaceType::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:129
Intrepid2::Basis_HGRAD_HEX_C1_FEM::Basis_HGRAD_HEX_C1_FEM
Basis_HGRAD_HEX_C1_FEM()
Constructor.
Definition: Intrepid2_HGRAD_HEX_C1_FEMDef.hpp:277
Generated by   doxygen 1.8.5