Intrepid2
Intrepid2_HGRAD_PYR_I2_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_I2_SERENDIPITY_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_PYR_I2_SERENDIPITY_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_I2_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 w = 1.0/(1.0 - z);
85 
86  output.access(0) = 0.25 * (-x - y - 1.0)*((1.0-x)*(1.0-y) - z + x*y*z*w);
87  output.access(1) = 0.25 * ( x - y - 1.0)*((1.0+x)*(1.0-y) - z - x*y*z*w);
88  output.access(2) = 0.25 * ( x + y - 1.0)*((1.0+x)*(1.0+y) - z + x*y*z*w);
89  output.access(3) = 0.25 * (-x + y - 1.0)*((1.0-x)*(1.0+y) - z - x*y*z*w);
90 
91  output.access(4) = z * (2.0*z - 1.0);
92 
93  output.access(5) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 - y - z)*w;
94  output.access(6) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 + x - z)*w;
95  output.access(7) = 0.5 * (1.0 + x - z)*(1.0 - x - z)*(1.0 + y - z)*w;
96  output.access(8) = 0.5 * (1.0 + y - z)*(1.0 - y - z)*(1.0 - x - z)*w;
97 
98  output.access(9) = z*(1.0 - x - z)*(1.0 - y - z)*w;
99  output.access(10) = z*(1.0 + x - z)*(1.0 - y - z)*w;
100  output.access(11) = z*(1.0 + x - z)*(1.0 + y - z)*w;
101  output.access(12) = z*(1.0 - x - z)*(1.0 + y - z)*w;
102 
103  break;
104  }
105  case OPERATOR_GRAD:
106  case OPERATOR_D1: {
107  const value_type w = 1.0/(1.0 - z);
108 
109  output.access(0, 0) = 0.25*(-1.0-x-y)*(-1.0+y + y*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
110  output.access(0, 1) = 0.25*(-1.0-x-y)*(-1.0+x + x*z*w) - 0.25*((1.0-x)*(1.0-y)-z + x*y*z*w);
111  output.access(0, 2) = 0.25*(-1.0-x-y)*(-1.0 + x*y*w + x*y*z*w*w);
112 
113  output.access(1, 0) = 0.25*(-1.0+x-y)*( 1.0-y - y*z*w) + 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
114  output.access(1, 1) = 0.25*(-1.0+x-y)*(-1.0-x - x*z*w) - 0.25*((1.0+x)*(1.0-y)-z - x*y*z*w);
115  output.access(1, 2) = 0.25*(-1.0+x-y)*(-1.0 - x*y*w - x*y*z*w*w);
116 
117  output.access(2, 0) = 0.25*(-1.0+x+y)*(1.0+y + y*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
118  output.access(2, 1) = 0.25*(-1.0+x+y)*(1.0+x + x*z*w) + 0.25*((1.0+x)*(1.0+y)-z + x*y*z*w);
119  output.access(2, 2) = 0.25*(-1.0+x+y)*(-1.0 + x*y*w + x*y*z*w*w);
120 
121  output.access(3, 0) = 0.25*(-1.0-x+y)*(-1.0-y - y*z*w) - 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
122  output.access(3, 1) = 0.25*(-1.0-x+y)*( 1.0-x - x*z*w) + 0.25*((1.0-x)*(1.0+y)-z - x*y*z*w);
123  output.access(3, 2) = 0.25*(-1.0-x+y)*(-1.0 - x*y*w - x*y*z*w*w);
124 
125  output.access(4, 0) = 0.0;
126  output.access(4, 1) = 0.0;
127  output.access(4, 2) = -1.0 + 4.0*z;
128 
129  output.access(5, 0) = -x*w*(1.0-y-z);
130  output.access(5, 1) = -0.5*(1.0-x-z)*(1.0+x-z)*w;
131  output.access(5, 2) = 0.5*y*x*x*w*w + 0.5*y - 1.0+z;
132 
133  output.access(6, 0) = 0.5*(1.0-y-z)*(1.0+y-z)*w;
134  output.access(6, 1) = -y*w*(1.0+x-z);
135  output.access(6, 2) = -0.5*x*y*y*w*w - 0.5*x - 1.0+z;
136 
137  output.access(7, 0) = -x*w*(1.0+y-z);
138  output.access(7, 1) = 0.5*(1.0-x-z)*(1.0+x-z)*w;
139  output.access(7, 2) = -0.5*y*x*x*w*w - 0.5*y - 1.0+z;
140 
141  output.access(8, 0) = -0.5*(1.0-y-z)*(1.0+y-z)*w;
142  output.access(8, 1) = -y*w*(1.0-x-z);
143  output.access(8, 2) = 0.5*x*y*y*w*w + 0.5*x - 1.0+z;
144 
145  output.access(9, 0) = -(1.0-y-z)*z*w;
146  output.access(9, 1) = -(1.0-x-z)*z*w;
147  output.access(9, 2) = x*y*w*w + 1.0 - x - y - 2.0*z;
148 
149  output.access(10,0) = (1.0-y-z)*z*w;
150  output.access(10,1) = -(1.0+x-z)*z*w;
151  output.access(10,2) = -x*y*w*w + 1.0 + x - y - 2.0*z;
152 
153  output.access(11,0) = (1.0+y-z)*z*w;
154  output.access(11,1) = (1.0+x-z)*z*w;
155  output.access(11,2) = x*y*w*w + 1.0 + x + y - 2.0*z;
156 
157  output.access(12,0) = -(1.0+y-z)*z*w;
158  output.access(12,1) = (1.0-x-z)*z*w;
159  output.access(12,2) = -x*y*w*w + 1.0 - x + y - 2.0*z;
160 
161  break;
162  }
163  case OPERATOR_D2: {
164  const value_type w = 1.0/(1.0 - z);
165 
166  output.access(0, 0) = -0.5*(-1.0+y+y*z*w);
167  output.access(0, 1) = -(-0.25 + 0.5*x + 0.5*y + 0.5*z)*w;
168  output.access(0, 2) = 0.25 + (-0.25*y-0.5*x*y-0.25*y*y)*w*w;
169  output.access(0, 3) = -0.5*(-1.0+x+x*z*w);
170  output.access(0, 4) = 0.25 + (-0.25*x*x-0.25*x-0.5*x*y)*w*w;
171  output.access(0, 5) = 0.5*x*y*(-1.0-x-y)*w*w*w;
172 
173  output.access(1, 0) = 0.5*(1.0-y-y*z*w);
174  output.access(1, 1) =-(0.25 + 0.5*x - 0.5*y - 0.5*z)*w;
175  output.access(1, 2) = -0.25 + (0.25*y-0.5*x*y+0.25*y*y)*w*w;
176  output.access(1, 3) = -0.5*(-1.0-x-x*z*w);
177  output.access(1, 4) = 0.25 + (-0.25*x*x + 0.25*x + 0.5*x*y)*w*w;
178  output.access(1, 5) = -0.5*x*y*(-1.0+x-y)*w*w*w;
179 
180  output.access(2, 0) = 0.5*(1.0+y+y*z*w);
181  output.access(2, 1) =-(-0.25 - 0.5*x - 0.5*y + 0.5*z)*w;
182  output.access(2, 2) = -0.25 + (-0.25*y+0.5*x*y+0.25*y*y)*w*w;
183  output.access(2, 3) = 0.5*(1.0+x+x*z*w);
184  output.access(2, 4) = -0.25 + (0.25*x*x -0.25*x + 0.5*x*y)*w*w;
185  output.access(2, 5) = 0.5*x*y*(-1.0+x+y)*w*w*w;
186 
187  output.access(3, 0) = -0.5*(-1.0-y-y*z*w);
188  output.access(3, 1) =-(0.25 - 0.5*x + 0.5*y - 0.5*z)*w;
189  output.access(3, 2) = 0.25 + (0.25*y+0.5*x*y-0.25*y*y)*w*w;
190  output.access(3, 3) = 0.5*(1.0-x-x*z*w);
191  output.access(3, 4) = -0.25 + (0.25*x + 0.25*x*x - 0.5*x*y)*w*w;
192  output.access(3, 5) = -0.5*x*y*(-1.0-x+y)*w*w*w;
193 
194  output.access(4, 0) = 0.0;
195  output.access(4, 1) = 0.0;
196  output.access(4, 2) = 0.0;
197  output.access(4, 3) = 0.0;
198  output.access(4, 4) = 0.0;
199  output.access(4, 5) = 4.0;
200 
201  output.access(5, 0) = -(1.0-y-z)*w;
202  output.access(5, 1) = x*w;
203  output.access(5, 2) = x*y*w*w;
204  output.access(5, 3) = 0.0;
205  output.access(5, 4) = 0.5*x*x*w*w + 0.5;
206  output.access(5, 5) = x*x*y*w*w*w + 1.0;
207 
208  output.access(6, 0) = 0.0;
209  output.access(6, 1) = -y*w;
210  output.access(6, 2) = -0.5*y*y*w*w - 0.5;
211  output.access(6, 3) =-(1.0+x-z)*w;
212  output.access(6, 4) = -x*y*w*w;
213  output.access(6, 5) = -x*y*y*w*w*w + 1.0;
214 
215  output.access(7, 0) = -(1.0+y-z)*w;
216  output.access(7, 1) = -x*w;
217  output.access(7, 2) = -x*y*w*w;
218  output.access(7, 3) = 0.0;
219  output.access(7, 4) = -0.5*x*x*w*w - 0.5;
220  output.access(7, 5) = -x*x*y*w*w*w + 1.0;
221 
222  output.access(8, 0) = 0.0;
223  output.access(8, 1) = y*w;
224  output.access(8, 2) = 0.5*y*y*w*w + 0.5;
225  output.access(8, 3) = -(1.0-x-z)*w;
226  output.access(8, 4) = x*y*w*w;
227  output.access(8, 5) = x*y*y*w*w*w + 1.0;
228 
229  output.access(9, 0) = 0.0;
230  output.access(9, 1) = z*w;
231  output.access(9, 2) = y*w*w - 1.0;
232  output.access(9, 3) = 0.0;
233  output.access(9, 4) = x*w*w - 1.0;
234  output.access(9, 5) = 2.0*x*y*w*w*w - 2.0;
235 
236  output.access(10,0) = 0.0;
237  output.access(10,1) = -z*w;
238  output.access(10,2) = -y*w*w + 1.0;
239  output.access(10,3) = 0.0;
240  output.access(10,4) = -x*w*w - 1.0;
241  output.access(10,5) = -2.0*x*y*w*w*w - 2.0;
242 
243  output.access(11,0) = 0.0;
244  output.access(11,1) = z*w;
245  output.access(11,2) = y*w*w + 1.0;
246  output.access(11,3) = 0.0;
247  output.access(11,4) = x*w*w + 1.0;
248  output.access(11,5) = 2.0*x*y*w*w*w - 2.0;
249 
250  output.access(12,0) = 0.0;
251  output.access(12,1) = -z*w;
252  output.access(12,2) = -y*w*w - 1.0;
253  output.access(12,3) = 0.0;
254  output.access(12,4) = -x*w*w + 1.0;
255  output.access(12,5) = -2.0*x*y*w*w*w - 2.0;
256 
257  break;
258  }
259  default: {
260  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
261  opType != OPERATOR_GRAD &&
262  opType != OPERATOR_D1 &&
263  opType != OPERATOR_D2,
264  ">>> ERROR: (Intrepid2::Basis_HGRAD_PYR_I2_FEM::Serial::getValues) operator is not supported");
265  }
266  }
267  }
268 
269  template<typename DT,
270  typename outputValueValueType, class ...outputValueProperties,
271  typename inputPointValueType, class ...inputPointProperties>
272  void
273  Basis_HGRAD_PYR_I2_FEM::
274  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
275  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
276  const EOperator operatorType ) {
277  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
278  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
279  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
280 
281  // Number of evaluation points = dim 0 of inputPoints
282  const auto loopSize = inputPoints.extent(0);
283  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
284 
285  switch (operatorType) {
286 
287  case OPERATOR_VALUE: {
288  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
289  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
290  break;
291  }
292  case OPERATOR_GRAD:
293  case OPERATOR_D1: {
294  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
295  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
296  break;
297  }
298  case OPERATOR_CURL: {
299  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
300  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
301  break;
302  }
303  case OPERATOR_DIV: {
304  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
305  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
306  break;
307  }
308  case OPERATOR_D2: {
309  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
310  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
311  break;
312  }
313  case OPERATOR_D3:
314  case OPERATOR_D4:
315  case OPERATOR_D5:
316  case OPERATOR_D6:
317  case OPERATOR_D7:
318  case OPERATOR_D8:
319  case OPERATOR_D9:
320  case OPERATOR_D10: {
321  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
322  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Operator not implemented yet");
323  break;
324  }
325  default: {
326  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
327  ">>> ERROR (Basis_HGRAD_PYR_I2_FEM): Invalid operator type");
328  }
329  }
330  }
331  }
332 
333  // -------------------------------------------------------------------------------------
334 
335  template<typename DT, typename OT, typename PT>
336  Basis_HGRAD_PYR_I2_FEM<DT,OT,PT>::
337  Basis_HGRAD_PYR_I2_FEM() {
338  this->basisCardinality_ = 13;
339  this->basisDegree_ = 2;
340  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() );
341  this->basisType_ = BASIS_FEM_DEFAULT;
342  this->basisCoordinates_ = COORDINATES_CARTESIAN;
343  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
344 
345  // initialize tags
346  {
347  // Basis-dependent intializations
348  const ordinal_type tagSize = 4; // size of DoF tag
349  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
350  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
351  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
352 
353  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
354  ordinal_type tags[52] = { 0, 0, 0, 1,
355  0, 1, 0, 1,
356  0, 2, 0, 1,
357  0, 3, 0, 1,
358  0, 4, 0, 1,
359  1, 0, 0, 1,
360  1, 1, 0, 1,
361  1, 2, 0, 1,
362  1, 3, 0, 1,
363  1, 4, 0, 1,
364  1, 5, 0, 1,
365  1, 6, 0, 1,
366  1, 7, 0, 1 };
367 
368 
369  // host tags
370  OrdinalTypeArray1DHost tagView(&tags[0], 52);
371 
372  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
373  this->setOrdinalTagData(this->tagToOrdinal_,
374  this->ordinalToTag_,
375  tagView,
376  this->basisCardinality_,
377  tagSize,
378  posScDim,
379  posScOrd,
380  posDfOrd);
381  }
382 
383  // dofCoords on host and create its mirror view to device
384  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
385  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
386 
387  dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0; dofCoords(0,2) = 0.0;
388  dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0; dofCoords(1,2) = 0.0;
389  dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
390  dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0; dofCoords(3,2) = 0.0;
391  dofCoords(4,0) = 0.0; dofCoords(4,1) = 0.0; dofCoords(4,2) = 1.0;
392 
393  dofCoords(5,0) = 0.0; dofCoords(5,1) = -1.0; dofCoords(5,2) = 0.0;
394  dofCoords(6,0) = 1.0; dofCoords(6,1) = 0.0; dofCoords(6,2) = 0.0;
395  dofCoords(7,0) = 0.0; dofCoords(7,1) = 1.0; dofCoords(7,2) = 0.0;
396  dofCoords(8,0) = -1.0; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.0;
397  dofCoords(9,0) = -0.5; dofCoords(9,1) = -0.5; dofCoords(9,2) = 0.5;
398  dofCoords(10,0) = 0.5; dofCoords(10,1) = -0.5; dofCoords(10,2) = 0.5;
399  dofCoords(11,0) = 0.5; dofCoords(11,1) = 0.5; dofCoords(11,2) = 0.5;
400  dofCoords(12,0) = -0.5; dofCoords(12,1) = 0.5; dofCoords(12,2) = 0.5;
401 
402  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
403  Kokkos::deep_copy(this->dofCoords_, dofCoords);
404  }
405 
406 }
407 
408 #endif
Intrepid2::Basis_HGRAD_PYR_I2_FEM::Basis_HGRAD_PYR_I2_FEM
Basis_HGRAD_PYR_I2_FEM()
Constructor.
Definition: Intrepid2_HGRAD_PYR_I2_FEMDef.hpp:337
Generated by   doxygen 1.8.5