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 DT,
228  typename outputValueValueType, class ...outputValueProperties,
229  typename inputPointValueType, class ...inputPointProperties>
230  void
231  Basis_HGRAD_TET_C2_FEM::
232  getValues( const typename DT::execution_space& space,
233  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
234  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
235  const EOperator operatorType ) {
236  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
237  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
238  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
239 
240  // Number of evaluation points = dim 0 of inputPoints
241  const auto loopSize = inputPoints.extent(0);
242  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
243 
244  switch (operatorType) {
245 
246  case OPERATOR_VALUE: {
247  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
248  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
249  break;
250  }
251  case OPERATOR_GRAD:
252  case OPERATOR_D1: {
253  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
254  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
255  break;
256  }
257  case OPERATOR_CURL: {
258  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
259  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
260  break;
261  }
262  case OPERATOR_DIV: {
263  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
264  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
265  break;
266  }
267  case OPERATOR_D2: {
268  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
269  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
270  break;
271  }
272  case OPERATOR_D3:
273  case OPERATOR_D4:
274  case OPERATOR_D5:
275  case OPERATOR_D6:
276  case OPERATOR_D7:
277  case OPERATOR_D8:
278  case OPERATOR_D9:
279  case OPERATOR_D10: {
280  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
281  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
282  break;
283  }
284  default: {
285  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
286  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
287  }
288  }
289  }
290  }
291  // -------------------------------------------------------------------------------------
292 
293  template<typename DT, typename OT, typename PT>
294  Basis_HGRAD_TET_C2_FEM<DT,OT,PT>::
295  Basis_HGRAD_TET_C2_FEM() {
296  this->basisCardinality_ = 10;
297  this->basisDegree_ = 2;
298  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
299  this->basisType_ = BASIS_FEM_DEFAULT;
300  this->basisCoordinates_ = COORDINATES_CARTESIAN;
301  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
302 
303  // intialize tags
304  {
305  // Basis-dependent intializations
306  const ordinal_type tagSize = 4; // size of DoF tag
307  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
308  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
309  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
310 
311  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
312  ordinal_type tags[40] = { 0, 0, 0, 1,
313  0, 1, 0, 1,
314  0, 2, 0, 1,
315  0, 3, 0, 1,
316  1, 0, 0, 1,
317  1, 1, 0, 1,
318  1, 2, 0, 1,
319  1, 3, 0, 1,
320  1, 4, 0, 1,
321  1, 5, 0, 1,
322  };
323 
324  // host tags
325  OrdinalTypeArray1DHost tagView(&tags[0], 40);
326 
327  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
328  this->setOrdinalTagData(this->tagToOrdinal_,
329  this->ordinalToTag_,
330  tagView,
331  this->basisCardinality_,
332  tagSize,
333  posScDim,
334  posScOrd,
335  posDfOrd);
336  }
337 
338  // dofCoords on host and create its mirror view to device
339  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
340  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
341 
342  dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
343  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
344  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
345  dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
346 
347  dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
348  dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
349  dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
350  dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
351  dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
352  dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
353 
354  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
355  Kokkos::deep_copy(this->dofCoords_, dofCoords);
356  }
357 
358 }// namespace Intrepid2
359 #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:295
Intrepid2::Basis_HGRAD_TET_C2_FEM::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:158
Generated by   doxygen 1.8.5