Intrepid2
Intrepid2_HGRAD_PYR_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_PYR_C1_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_PYR_C1_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_PYR_C1_FEM::Serial<opType>::
64  getValues( outputViewType output,
65  const inputViewType input ) {
66  const auto eps = epsilon();
67 
68  static_assert(std::is_same<
69  typename outputViewType::value_type,
70  typename inputViewType::value_type>::value,"Input/output view has different value types");
71 
72  typedef typename outputViewType::value_type value_type;
73 
74  const value_type x = input(0);
75  const value_type y = input(1);
76  const value_type ztmp = input(2);
77 
78  //be sure that the basis functions are defined when z is very close to 1.
79  const value_type z = ( (value_type(1.0) - ztmp) < value_type(eps) ? value_type(1.0 - eps) : ztmp );
80 
81  switch (opType) {
82 
83  case OPERATOR_VALUE: {
84  const value_type factor = 0.25/(1.0 - z);
85 
86  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
87  output.access(0) = (1.0 - x - z) * (1.0 - y - z) * factor;
88  output.access(1) = (1.0 + x - z) * (1.0 - y - z) * factor;
89  output.access(2) = (1.0 + x - z) * (1.0 + y - z) * factor;
90  output.access(3) = (1.0 - x - z) * (1.0 + y - z) * factor;
91  output.access(4) = z;
92  break;
93  }
94  case OPERATOR_GRAD: {
95  const value_type factor = 0.25/(1.0 - z);
96  const value_type factor2 = 4.0 * factor * factor;
97 
98  // output.accessValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
99  output.access(0, 0) = (y + z - 1.0) * factor;
100  output.access(0, 1) = (x + z - 1.0) * factor;
101  output.access(0, 2) = x * y * factor2 - 0.25;
102 
103  output.access(1, 0) = (1.0 - y - z) * factor;
104  output.access(1, 1) = (z - x - 1.0) * factor;
105  output.access(1, 2) = - x*y * factor2 - 0.25;
106 
107  output.access(2, 0) = (1.0 + y - z) * factor;
108  output.access(2, 1) = (1.0 + x - z) * factor;
109  output.access(2, 2) = x * y * factor2 - 0.25;
110 
111  output.access(3, 0) = (z - y - 1.0) * factor;
112  output.access(3, 1) = (1.0 - x - z) * factor;
113  output.access(3, 2) = - x*y * factor2 - 0.25;
114 
115  output.access(4, 0) = 0.0;
116  output.access(4, 1) = 0.0;
117  output.access(4, 2) = 1;
118  break;
119  }
120  case OPERATOR_D2: {
121  const value_type factor = 0.25/(1.0 - z);
122  const value_type factor2 = 4.0 * factor * factor;
123  const value_type factor3 = 8.0 * factor * factor2;
124 
125  // output.accessValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6)
126  output.access(0, 0) = 0.0; // {2, 0, 0}
127  output.access(0, 1) = factor; // {1, 1, 0}
128  output.access(0, 2) = y*factor2; // {1, 0, 1}
129  output.access(0, 3) = 0.0; // {0, 2, 0}
130  output.access(0, 4) = x*factor2; // {0, 1, 1}
131  output.access(0, 5) = x*y*factor3; // {0, 0, 2}
132 
133  output.access(1, 0) = 0.0; // {2, 0, 0}
134  output.access(1, 1) = -factor; // {1, 1, 0}
135  output.access(1, 2) = -y*factor2; // {1, 0, 1}
136  output.access(1, 3) = 0.0; // {0, 2, 0}
137  output.access(1, 4) = -x*factor2; // {0, 1, 1}
138  output.access(1, 5) = -x*y*factor3; // {0, 0, 2}
139 
140  output.access(2, 0) = 0.0; // {2, 0, 0}
141  output.access(2, 1) = factor; // {1, 1, 0}
142  output.access(2, 2) = y*factor2; // {1, 0, 1}
143  output.access(2, 3) = 0.0; // {0, 2, 0}
144  output.access(2, 4) = x*factor2; // {0, 1, 1}
145  output.access(2, 5) = x*y*factor3; // {0, 0, 2}
146 
147  output.access(3, 0) = 0.0; // {2, 0, 0}
148  output.access(3, 1) = -factor; // {1, 1, 0}
149  output.access(3, 2) = -y*factor2; // {1, 0, 1}
150  output.access(3, 3) = 0.0; // {0, 2, 0}
151  output.access(3, 4) = -x*factor2; // {0, 1, 1}
152  output.access(3, 5) = -x*y*factor3; // {0, 0, 2}
153 
154  output.access(4, 0) = 0.0; // {2, 0, 0}
155  output.access(4, 1) = 0.0; // {1, 1, 0}
156  output.access(4, 2) = 0.0; // {1, 0, 1}
157  output.access(4, 3) = 0.0; // {0, 2, 0}
158  output.access(4, 4) = 0.0; // {0, 1, 1}
159  output.access(4, 5) = 0.0; // {0, 0, 2}
160  break;
161  }
162  case OPERATOR_MAX: {
163  const ordinal_type jend = output.extent(1);
164  const ordinal_type iend = output.extent(0);
165 
166  for (ordinal_type j=0;j<jend;++j)
167  for (ordinal_type i=0;i<iend;++i)
168  output.access(i, j) = 0.0;
169  break;
170  }
171  default: {
172  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
173  opType != OPERATOR_GRAD &&
174  opType != OPERATOR_D2 &&
175  opType != OPERATOR_MAX,
176  ">>> ERROR: (Intrepid2::Basis_HGRAD_PYR_C1_FEM::Serial::getValues) operator is not supported");
177  }
178  }
179  }
180 
181  template<typename SpT,
182  typename outputValueValueType, class ...outputValueProperties,
183  typename inputPointValueType, class ...inputPointProperties>
184  void
185  Basis_HGRAD_PYR_C1_FEM::
186  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
187  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
188  const EOperator operatorType ) {
189  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
190  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
191  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
192 
193  // Number of evaluation points = dim 0 of inputPoints
194  const auto loopSize = inputPoints.extent(0);
195  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
196 
197  switch (operatorType) {
198 
199  case OPERATOR_VALUE: {
200  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
201  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
202  break;
203  }
204  case OPERATOR_GRAD:
205  case OPERATOR_D1: {
206  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
207  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
208  break;
209  }
210  case OPERATOR_CURL: {
211  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
212  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
213  break;
214  }
215  case OPERATOR_DIV: {
216  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
217  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
218  break;
219  }
220  case OPERATOR_D2: {
221  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
222  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
223  break;
224  }
225  case OPERATOR_D3:
226  case OPERATOR_D4:
227  case OPERATOR_D5:
228  case OPERATOR_D6:
229  case OPERATOR_D7:
230  case OPERATOR_D8:
231  case OPERATOR_D9:
232  case OPERATOR_D10: {
233  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
234  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
235  break;
236  }
237  default: {
238  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
239  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): Invalid operator type");
240  }
241  }
242  }
243  }
244 
245  // -------------------------------------------------------------------------------------
246 
247  template<typename SpT, typename OT, typename PT>
250  this->basisCardinality_ = 5;
251  this->basisDegree_ = 1;
252  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() );
253  this->basisType_ = BASIS_FEM_DEFAULT;
254  this->basisCoordinates_ = COORDINATES_CARTESIAN;
255 
256  // initialize tags
257  {
258  // Basis-dependent intializations
259  const ordinal_type tagSize = 4; // size of DoF tag
260  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
261  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
262  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
263 
264  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
265  ordinal_type tags[20] = { 0, 0, 0, 1,
266  0, 1, 0, 1,
267  0, 2, 0, 1,
268  0, 3, 0, 1,
269  0, 4, 0, 1 };
270 
271 
272  // host tags
273  ordinal_type_array_1d_host tagView(&tags[0], 20);
274 
275  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
276  //ordinal_type_array_2d_host ordinalToTag;
277  //ordinal_type_array_3d_host tagToOrdinal;
278  this->setOrdinalTagData(this->tagToOrdinal_,
279  this->ordinalToTag_,
280  tagView,
281  this->basisCardinality_,
282  tagSize,
283  posScDim,
284  posScOrd,
285  posDfOrd);
286 
287  //this->tagToOrdinal_ = Kokkos::create_mirror_view(typename SpT::memory_space(), tagToOrdinal);
288  //Kokkos::deep_copy(this->tagToOrdinal_, tagToOrdinal);
289 
290  //this->ordinalToTag_ = Kokkos::create_mirror_view(typename SpT::memory_space(), ordinalToTag);
291  //Kokkos::deep_copy(this->ordinalToTag_, ordinalToTag);
292  }
293 
294  // dofCoords on host and create its mirror view to device
295  Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
296  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
297 
298  dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = 0.0;
299  dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = 0.0;
300  dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
301  dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = 0.0;
302  dofCoords(4,0) = 0.0; dofCoords(4,1) = 0.0; dofCoords(4,2) = 1.0;
303 
304  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
305  Kokkos::deep_copy(this->dofCoords_, dofCoords);
306  }
307 
308 }
309 
310 #endif
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.