Intrepid2
Intrepid2_HGRAD_QUAD_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_QUAD_C2_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_QUAD_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_QUAD_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 
71  // output is a rank-2 array with dimensions (basisCardinality_, dim0)
72  output.access(0) = x*(x - 1.0)*y*(y - 1.0)/4.0;
73  output.access(1) = x*(x + 1.0)*y*(y - 1.0)/4.0;
74  output.access(2) = x*(x + 1.0)*y*(y + 1.0)/4.0;
75  output.access(3) = x*(x - 1.0)*y*(y + 1.0)/4.0;
76  // edge midpoints basis functions
77  output.access(4) = (1.0 - x)*(1.0 + x)*y*(y - 1.0)/2.0;
78  output.access(5) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0;
79  output.access(6) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0;
80  output.access(7) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0;
81  // quad bubble basis function
82  output.access(8) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y);
83  break;
84  }
85  case OPERATOR_D1 :
86  case OPERATOR_GRAD : {
87  const auto x = input(0);
88  const auto y = input(1);
89 
90  // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
91  output.access(0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y;
92  output.access(0, 1) = (-1.0 + x)*x*(-0.25 + 0.5*y);
93 
94  output.access(1, 0) = (0.25 + 0.5*x)*(-1. + y)*y;
95  output.access(1, 1) = x*(1. + x)*(-0.25 + 0.5*y);
96 
97  output.access(2, 0) = (0.25 + 0.5*x)*y*(1. + y);
98  output.access(2, 1) = x*(1. + x)*(0.25 + 0.5*y);
99 
100  output.access(3, 0) = (-0.25 + 0.5*x)*y*(1. + y);
101  output.access(3, 1) = (-1. + x)*x*(0.25 + 0.5*y);
102 
103  output.access(4, 0) = x*(1.0 - y)*y;
104  output.access(4, 1) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
105 
106  output.access(5, 0) = 0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
107  output.access(5, 1) =-x*(1.0 + x)*y;
108 
109  output.access(6, 0) =-y*(1.0 + y)*x;
110  output.access(6, 1) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
111 
112  output.access(7, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x);
113  output.access(7, 1) = (1.0 - x)*x*y;
114 
115  output.access(8, 0) =-2.0*(1.0 - y)*(1.0 + y)*x;
116  output.access(8, 1) =-2.0*(1.0 - x)*(1.0 + x)*y;
117  break;
118  }
119  case OPERATOR_CURL : {
120  const auto x = input(0);
121  const auto y = input(1);
122 
123  // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
124  // CURL(u) = (u_y, -u_x), is rotated GRAD
125  output.access(0, 1) =-(-0.25 + 0.5*x)*(-1. + y)*y;
126  output.access(0, 0) = (-1.0 + x)*x*(-0.25 + 0.5*y);
127 
128  output.access(1, 1) =-(0.25 + 0.5*x)*(-1. + y)*y;
129  output.access(1, 0) = x*(1. + x)*(-0.25 + 0.5*y);
130 
131  output.access(2, 1) =-(0.25 + 0.5*x)*y*(1. + y);
132  output.access(2, 0) = x*(1. + x)*(0.25 + 0.5*y);
133 
134  output.access(3, 1) =-(-0.25 + 0.5*x)*y*(1. + y);
135  output.access(3, 0) = (-1. + x)*x*(0.25 + 0.5*y);
136 
137  output.access(4, 1) =-x*(1.0 - y)*y;
138  output.access(4, 0) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
139 
140  output.access(5, 1) =-0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
141  output.access(5, 0) =-x*(1.0 + x)*y;
142 
143  output.access(6, 1) = y*(1.0 + y)*x;
144  output.access(6, 0) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
145 
146  output.access(7, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x);
147  output.access(7, 0) = (1.0 - x)*x*y;
148 
149  output.access(8, 1) = 2.0*(1.0 - y)*(1.0 + y)*x;
150  output.access(8, 0) =-2.0*(1.0 - x)*(1.0 + x)*y;
151  break;
152  }
153  case OPERATOR_D2 : {
154  const auto x = input(0);
155  const auto y = input(1);
156  // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality=3)
157  output.access(0, 0) = 0.5*(-1.0 + y)*y;
158  output.access(0, 1) = 0.25 - 0.5*y + x*(-0.5 + 1.*y);
159  output.access(0, 2) = 0.5*(-1.0 + x)*x;
160 
161  output.access(1, 0) = 0.5*(-1.0 + y)*y;
162  output.access(1, 1) =-0.25 + 0.5*y + x*(-0.5 + 1.*y);
163  output.access(1, 2) = 0.5*x*(1.0 + x);
164 
165  output.access(2, 0) = 0.5*y*(1.0 + y);
166  output.access(2, 1) = 0.25 + 0.5*y + x*(0.5 + 1.*y);
167  output.access(2, 2) = 0.5*x*(1.0 + x);
168 
169  output.access(3, 0) = 0.5*y*(1.0 + y);
170  output.access(3, 1) =-0.25 - 0.5*y + x*(0.5 + 1.*y);
171  output.access(3, 2) = 0.5*(-1.0 + x)*x;
172 
173  output.access(4, 0) = (1.0 - y)*y;
174  output.access(4, 1) = x*(1. - 2.*y);
175  output.access(4, 2) = (1.0 - x)*(1.0 + x);
176 
177  output.access(5, 0) = (1.0 - y)*(1.0 + y);
178  output.access(5, 1) = x*(0. - 2.*y) - 1.*y;
179  output.access(5, 2) =-x*(1.0 + x);
180 
181  output.access(6, 0) =-y*(1.0 + y);
182  output.access(6, 1) = x*(-1. - 2.*y);
183  output.access(6, 2) = (1.0 - x)*(1.0 + x);
184 
185  output.access(7, 0) = (1.0 - y)*(1.0 + y);
186  output.access(7, 1) = x*(0. - 2.*y) + 1.*y;
187  output.access(7, 2) = (1.0 - x)*x;
188 
189  output.access(8, 0) =-2.0 + 2.0*y*y;
190  output.access(8, 1) = 4*x*y;
191  output.access(8, 2) =-2.0 + 2.0*x*x;
192  break;
193  }
194  case OPERATOR_D3 : {
195  const auto x = input(0);
196  const auto y = input(1);
197  output.access(0, 0) = 0.0;
198  output.access(0, 1) =-0.5 + y;
199  output.access(0, 2) =-0.5 + x;
200  output.access(0, 3) = 0.0;
201 
202  output.access(1, 0) = 0.0;
203  output.access(1, 1) =-0.5 + y;
204  output.access(1, 2) = 0.5 + x;
205  output.access(1, 3) = 0.0;
206 
207  output.access(2, 0) = 0.0;
208  output.access(2, 1) = 0.5 + y;
209  output.access(2, 2) = 0.5 + x;
210  output.access(2, 3) = 0.0;
211 
212  output.access(3, 0) = 0.0;
213  output.access(3, 1) = 0.5 + y;
214  output.access(3, 2) =-0.5 + x;
215  output.access(3, 3) = 0.0;
216 
217  output.access(4, 0) = 0.0;
218  output.access(4, 1) = 1.0 - 2.0*y;
219  output.access(4, 2) =-2.0*x;
220  output.access(4, 3) = 0.0;
221 
222  output.access(5, 0) = 0.0;
223  output.access(5, 1) =-2.0*y;
224  output.access(5, 2) =-1.0 - 2.0*x;
225  output.access(5, 3) = 0.0;
226 
227  output.access(6, 0) = 0.0;
228  output.access(6, 1) =-1.0 - 2.0*y;
229  output.access(6, 2) =-2.0*x;
230  output.access(6, 3) = 0.0;
231 
232  output.access(7, 0) = 0.0;
233  output.access(7, 1) =-2.0*y;
234  output.access(7, 2) = 1.0 - 2.0*x;
235  output.access(7, 3) = 0.0;
236 
237  output.access(8, 0) = 0.0;
238  output.access(8, 1) = 4.0*y;
239  output.access(8, 2) = 4.0*x;
240  output.access(8, 3) = 0.0;
241  break;
242  }
243  case OPERATOR_D4 : {
244  output.access(0, 0) = 0.0;
245  output.access(0, 1) = 0.0;
246  output.access(0, 2) = 1.0;
247  output.access(0, 3) = 0.0;
248  output.access(0, 4) = 0.0;
249 
250  output.access(1, 0) = 0.0;
251  output.access(1, 1) = 0.0;
252  output.access(1, 2) = 1.0;
253  output.access(1, 3) = 0.0;
254  output.access(1, 4) = 0.0;
255 
256  output.access(2, 0) = 0.0;
257  output.access(2, 1) = 0.0;
258  output.access(2, 2) = 1.0;
259  output.access(2, 3) = 0.0;
260  output.access(2, 4) = 0.0;
261 
262  output.access(3, 0) = 0.0;
263  output.access(3, 1) = 0.0;
264  output.access(3, 2) = 1.0;
265  output.access(3, 3) = 0.0;
266  output.access(3, 4) = 0.0;
267 
268  output.access(4, 0) = 0.0;
269  output.access(4, 1) = 0.0;
270  output.access(4, 2) =-2.0;
271  output.access(4, 3) = 0.0;
272  output.access(4, 4) = 0.0;
273 
274  output.access(5, 0) = 0.0;
275  output.access(5, 1) = 0.0;
276  output.access(5, 2) =-2.0;
277  output.access(5, 3) = 0.0;
278  output.access(5, 4) = 0.0;
279 
280  output.access(6, 0) = 0.0;
281  output.access(6, 1) = 0.0;
282  output.access(6, 2) =-2.0;
283  output.access(6, 3) = 0.0;
284  output.access(6, 4) = 0.0;
285 
286  output.access(7, 0) = 0.0;
287  output.access(7, 1) = 0.0;
288  output.access(7, 2) =-2.0;
289  output.access(7, 3) = 0.0;
290  output.access(7, 4) = 0.0;
291 
292  output.access(8, 0) = 0.0;
293  output.access(8, 1) = 0.0;
294  output.access(8, 2) = 4.0;
295  output.access(8, 3) = 0.0;
296  output.access(8, 4) = 0.0;
297  break;
298  }
299  case OPERATOR_MAX : {
300  const ordinal_type jend = output.extent(1);
301  const ordinal_type iend = output.extent(0);
302 
303  for (ordinal_type j=0;j<jend;++j)
304  for (ordinal_type i=0;i<iend;++i)
305  output.access(i, j) = 0.0;
306  break;
307  }
308  default: {
309  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
310  opType != OPERATOR_GRAD &&
311  opType != OPERATOR_CURL &&
312  opType != OPERATOR_D1 &&
313  opType != OPERATOR_D2 &&
314  opType != OPERATOR_D3 &&
315  opType != OPERATOR_D4 &&
316  opType != OPERATOR_MAX,
317  ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_C2_FEM::Serial::getValues) operator is not supported");
318 
319  }
320  }
321  }
322 
323  template<typename SpT,
324  typename outputValueValueType, class ...outputValueProperties,
325  typename inputPointValueType, class ...inputPointProperties>
326  void
327  Basis_HGRAD_QUAD_C2_FEM::
328  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
329  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
330  const EOperator operatorType ) {
331  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
332  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
333  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
334 
335  // Number of evaluation points = dim 0 of inputPoints
336  const auto loopSize = inputPoints.extent(0);
337  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
338 
339  switch (operatorType) {
340 
341  case OPERATOR_VALUE: {
342  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
343  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
344  break;
345  }
346  case OPERATOR_GRAD:
347  case OPERATOR_D1: {
348  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
349  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
350  break;
351  }
352  case OPERATOR_CURL: {
353  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
354  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
355  break;
356  }
357  case OPERATOR_DIV: {
358  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
359  ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
360  break;
361  }
362  case OPERATOR_D2: {
363  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
364  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
365  break;
366  }
367  case OPERATOR_D3: {
368  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D3Cardinality=4)
369  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D3> FunctorType;
370  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
371  break;
372  }
373  case OPERATOR_D4: {
374  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D4Cardinality=5)
375  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D4> FunctorType;
376  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
377  break;
378  }
379  case OPERATOR_D5:
380  case OPERATOR_D6:
381  case OPERATOR_D7:
382  case OPERATOR_D8:
383  case OPERATOR_D9:
384  case OPERATOR_D10: {
385  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
386  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
387  break;
388  }
389  default: {
390  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
391  ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
392  }
393  }
394  }
395 
396 
397 
398 
399  }
400  // -------------------------------------------------------------------------------------
401 
402 
403  template<typename SpT, typename OT, typename PT>
404  Basis_HGRAD_QUAD_C2_FEM<SpT,OT,PT>::
405  Basis_HGRAD_QUAD_C2_FEM() {
406  this->basisCardinality_ = 9;
407  this->basisDegree_ = 2;
408  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
409  this->basisType_ = BASIS_FEM_DEFAULT;
410  this->basisCoordinates_ = COORDINATES_CARTESIAN;
411 
412  {
413  // Basis-dependent intializations
414  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
415  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
416  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
417  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
418 
419  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
420  ordinal_type tags[36] = { 0, 0, 0, 1,
421  0, 1, 0, 1,
422  0, 2, 0, 1,
423  0, 3, 0, 1,
424  // edge midpoints
425  1, 0, 0, 1,
426  1, 1, 0, 1,
427  1, 2, 0, 1,
428  1, 3, 0, 1,
429  // quad center
430  2, 0, 0, 1};
431 
432  //host view
433  ordinal_type_array_1d_host tagView(&tags[0],36);
434 
435  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
436  this->setOrdinalTagData(this->tagToOrdinal_,
437  this->ordinalToTag_,
438  tagView,
439  this->basisCardinality_,
440  tagSize,
441  posScDim,
442  posScOrd,
443  posDfOrd);
444  }
445 
446  // dofCoords on host and create its mirror view to device
447  Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
448  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
449 
450  dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0;
451  dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0;
452  dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0;
453  dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0;
454 
455  dofCoords(4,0) = 0.0; dofCoords(4,1) = -1.0;
456  dofCoords(5,0) = 1.0; dofCoords(5,1) = 0.0;
457  dofCoords(6,0) = 0.0; dofCoords(6,1) = 1.0;
458  dofCoords(7,0) = -1.0; dofCoords(7,1) = 0.0;
459 
460  dofCoords(8,0) = 0.0; dofCoords(8,1) = 0.0;
461 
462  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
463  Kokkos::deep_copy(this->dofCoords_, dofCoords);
464  }
465 
466 }// namespace Intrepid2
467 #endif
Intrepid2::Basis_HGRAD_QUAD_C2_FEM::Basis_HGRAD_QUAD_C2_FEM
Basis_HGRAD_QUAD_C2_FEM()
Constructor.
Definition: Intrepid2_HGRAD_QUAD_C2_FEMDef.hpp:405
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