Intrepid2
Intrepid2_HGRAD_TET_C2_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_TET_C2_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_TET_C2_FEM_DEF_HPP__
51 
52 namespace Intrepid2 {
53 
54  // -------------------------------------------------------------------------------------
55 
56  namespace Impl {
57 
58  template<EOperator opType>
59  template<typename outputViewType,
60  typename inputViewType>
61  KOKKOS_INLINE_FUNCTION
62  void
63  Basis_HGRAD_TET_C2_FEM::Serial<opType>::
64  getValues( outputViewType output,
65  const inputViewType input ) {
66  switch (opType) {
67  case OPERATOR_VALUE: {
68  const auto x = input(0);
69  const auto y = input(1);
70  const auto z = input(2);
71 
72  // output is a rank-2 array with dimensions (basisCardinality_, dim0)
73  output.access(0) = (-1. + x + y + z)*(-1. + 2.*x + 2.*y + 2.*z);
74  output.access(1) = x*(-1. + 2.*x);
75  output.access(2) = y*(-1. + 2.*y);
76  output.access(3) = z*(-1. + 2.*z);
77 
78  output.access(4) = -4.*x*(-1. + x + y + z);
79  output.access(5) = 4.*x*y;
80  output.access(6) = -4.*y*(-1. + x + y + z);
81  output.access(7) = -4.*z*(-1. + x + y + z);
82  output.access(8) = 4.*x*z;
83  output.access(9) = 4.*y*z;
84  break;
85  }
86  case OPERATOR_D1:
87  case OPERATOR_GRAD: {
88  const auto x = input(0);
89  const auto y = input(1);
90  const auto z = input(2);
91 
92  // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
93  output.access(0, 0) = -3.+ 4.*x + 4.*y + 4.*z;
94  output.access(0, 1) = -3.+ 4.*x + 4.*y + 4.*z;
95  output.access(0, 2) = -3.+ 4.*x + 4.*y + 4.*z;
96 
97  output.access(1, 0) = -1.+ 4.*x;
98  output.access(1, 1) = 0.;
99  output.access(1, 2) = 0.;
100 
101  output.access(2, 0) = 0.;
102  output.access(2, 1) = -1.+ 4.*y;
103  output.access(2, 2) = 0.;
104 
105  output.access(3, 0) = 0.;
106  output.access(3, 1) = 0.;
107  output.access(3, 2) = -1.+ 4.*z;
108 
109 
110  output.access(4, 0) = -4.*(-1.+ 2*x + y + z);
111  output.access(4, 1) = -4.*x;
112  output.access(4, 2) = -4.*x;
113 
114  output.access(5, 0) = 4.*y;
115  output.access(5, 1) = 4.*x;
116  output.access(5, 2) = 0.;
117 
118  output.access(6, 0) = -4.*y;
119  output.access(6, 1) = -4.*(-1.+ x + 2*y + z);
120  output.access(6, 2) = -4.*y;
121 
122  output.access(7, 0) = -4.*z;
123  output.access(7, 1) = -4.*z;
124  output.access(7, 2) = -4.*(-1.+ x + y + 2*z);
125 
126  output.access(8, 0) = 4.*z;
127  output.access(8, 1) = 0.;
128  output.access(8, 2) = 4.*x;
129 
130  output.access(9, 0) = 0.;
131  output.access(9, 1) = 4.*z;
132  output.access(9, 2) = 4.*y;
133  break;
134  }
135  case OPERATOR_D2: {
136  output.access(0, 0) = 4.;
137  output.access(0, 1) = 4.;
138  output.access(0, 2) = 4.;
139  output.access(0, 3) = 4.;
140  output.access(0, 4) = 4.;
141  output.access(0, 5) = 4.;
142 
143  output.access(1, 0) = 4.;
144  output.access(1, 1) = 0.;
145  output.access(1, 2) = 0.;
146  output.access(1, 3) = 0.;
147  output.access(1, 4) = 0.;
148  output.access(1, 5) = 0.;
149 
150  output.access(2, 0) = 0.;
151  output.access(2, 1) = 0.;
152  output.access(2, 2) = 0.;
153  output.access(2, 3) = 4.;
154  output.access(2, 4) = 0.;
155  output.access(2, 5) = 0.;
156 
157  output.access(3, 0) = 0.;
158  output.access(3, 1) = 0.;
159  output.access(3, 2) = 0.;
160  output.access(3, 3) = 0.;
161  output.access(3, 4) = 0.;
162  output.access(3, 5) = 4.;
163 
164  output.access(4, 0) = -8.;
165  output.access(4, 1) = -4.;
166  output.access(4, 2) = -4.;
167  output.access(4, 3) = 0.;
168  output.access(4, 4) = 0.;
169  output.access(4, 5) = 0.;
170 
171  output.access(5, 0) = 0.;
172  output.access(5, 1) = 4.;
173  output.access(5, 2) = 0.;
174  output.access(5, 3) = 0.;
175  output.access(5, 4) = 0.;
176  output.access(5, 5) = 0.;
177 
178  output.access(6, 0) = 0.;
179  output.access(6, 1) = -4.;
180  output.access(6, 2) = 0.;
181  output.access(6, 3) = -8.;
182  output.access(6, 4) = -4.;
183  output.access(6, 5) = 0;
184 
185  output.access(7, 0) = 0.;
186  output.access(7, 1) = 0.;
187  output.access(7, 2) = -4.;
188  output.access(7, 3) = 0.;
189  output.access(7, 4) = -4.;
190  output.access(7, 5) = -8.;
191 
192  output.access(8, 0) = 0.;
193  output.access(8, 1) = 0.;
194  output.access(8, 2) = 4.;
195  output.access(8, 3) = 0.;
196  output.access(8, 4) = 0.;
197  output.access(8, 5) = 0.;
198 
199  output.access(9, 0) = 0.;
200  output.access(9, 1) = 0.;
201  output.access(9, 2) = 0.;
202  output.access(9, 3) = 0.;
203  output.access(9, 4) = 4.;
204  output.access(9, 5) = 0.;
205  break;
206  }
207  case OPERATOR_MAX: {
208  const ordinal_type jend = output.extent(1);
209  const ordinal_type iend = output.extent(0);
210 
211  for (ordinal_type j=0;j<jend;++j)
212  for (ordinal_type i=0;i<iend;++i)
213  output.access(i, j) = 0.0;
214  break;
215  }
216  default: {
217  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
218  opType != OPERATOR_GRAD &&
219  opType != OPERATOR_D1 &&
220  opType != OPERATOR_D2 &&
221  opType != OPERATOR_MAX,
222  ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_C2_FEM::Serial::getValues) operator is not supported");
223  }
224  }
225  }
226 
227  template<typename SpT,
228  typename outputValueValueType, class ...outputValueProperties,
229  typename inputPointValueType, class ...inputPointProperties>
230  void
231  Basis_HGRAD_TET_C2_FEM::
232  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
233  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
234  const EOperator operatorType ) {
235  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
236  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
237  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
238 
239  // Number of evaluation points = dim 0 of inputPoints
240  const auto loopSize = inputPoints.extent(0);
241  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
242 
243  switch (operatorType) {
244 
245  case OPERATOR_VALUE: {
246  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
247  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
248  break;
249  }
250  case OPERATOR_GRAD:
251  case OPERATOR_D1: {
252  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
253  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
254  break;
255  }
256  case OPERATOR_CURL: {
257  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
258  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
259  break;
260  }
261  case OPERATOR_DIV: {
262  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
263  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
264  break;
265  }
266  case OPERATOR_D2: {
267  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
268  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
269  break;
270  }
271  case OPERATOR_D3:
272  case OPERATOR_D4:
273  case OPERATOR_D5:
274  case OPERATOR_D6:
275  case OPERATOR_D7:
276  case OPERATOR_D8:
277  case OPERATOR_D9:
278  case OPERATOR_D10: {
279  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
280  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
281  break;
282  }
283  default: {
284  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
285  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
286  }
287  }
288  }
289  }
290  // -------------------------------------------------------------------------------------
291 
292  template<typename SpT, typename OT, typename PT>
293  Basis_HGRAD_TET_C2_FEM<SpT,OT,PT>::
294  Basis_HGRAD_TET_C2_FEM() {
295  this->basisCardinality_ = 10;
296  this->basisDegree_ = 2;
297  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
298  this->basisType_ = BASIS_FEM_DEFAULT;
299  this->basisCoordinates_ = COORDINATES_CARTESIAN;
300 
301  // intialize tags
302  {
303  // Basis-dependent intializations
304  const ordinal_type tagSize = 4; // size of DoF tag
305  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
306  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
307  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
308 
309  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
310  ordinal_type tags[40] = { 0, 0, 0, 1,
311  0, 1, 0, 1,
312  0, 2, 0, 1,
313  0, 3, 0, 1,
314  1, 0, 0, 1,
315  1, 1, 0, 1,
316  1, 2, 0, 1,
317  1, 3, 0, 1,
318  1, 4, 0, 1,
319  1, 5, 0, 1,
320  };
321 
322  // host tags
323  ordinal_type_array_1d_host tagView(&tags[0], 40);
324 
325  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
326  this->setOrdinalTagData(this->tagToOrdinal_,
327  this->ordinalToTag_,
328  tagView,
329  this->basisCardinality_,
330  tagSize,
331  posScDim,
332  posScOrd,
333  posDfOrd);
334  }
335 
336  // dofCoords on host and create its mirror view to device
337  Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
338  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
339 
340  dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
341  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
342  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
343  dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
344 
345  dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
346  dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
347  dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
348  dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
349  dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
350  dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
351 
352  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
353  Kokkos::deep_copy(this->dofCoords_, dofCoords);
354  }
355 
356 }// namespace Intrepid2
357 #endif
Intrepid2::Basis_HGRAD_TET_C2_FEM::Basis_HGRAD_TET_C2_FEM
Basis_HGRAD_TET_C2_FEM()
Constructor.
Definition: Intrepid2_HGRAD_TET_C2_FEMDef.hpp:294
Intrepid2::Basis::ordinal_type_array_1d_host
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:113
Generated by   doxygen 1.8.5