Intrepid2
Intrepid2_HGRAD_TET_COMP12_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 
50 #ifndef __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
51 #define __INTREPID2_HGRAD_TET_COMP12_FEMDEF_HPP__
52 
53 namespace Intrepid2 {
54 
55  // -------------------------------------------------------------------------------------
56  namespace Impl {
57 
58  // assume that subtets are disjoint and a single point belong to one subtet
59  // at the interface, it returns the first one that satisfy the condition
60  template<typename pointValueType>
61  KOKKOS_INLINE_FUNCTION
62  ordinal_type
63  Basis_HGRAD_TET_COMP12_FEM::getLocalSubTet( const pointValueType x,
64  const pointValueType y,
65  const pointValueType z ) {
66 
67  const pointValueType
68  xyz = x + y + z,
69  xy = x + y,
70  xz = x + z,
71  yz = y + z;
72 
73  // cycle through each subdomain and push back if the point lies within
74 
75  // subtet #0 E0 := 0.0 <= r + s + t <= 0.5 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
76  if ( (0.0 <= xyz && xyz <= 0.5) &&
77  (0.0 <= x && x <= 0.5) &&
78  (0.0 <= y && y <= 0.5) &&
79  (0.0 <= z && z <= 0.5) )
80  return 0;
81 
82  // subtet #1 E1 := 0.5 <= r + s + t <= 1.0 && 0.5 <= r <= 1.0 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
83  if ( (0.5 <= xyz && xyz <= 1.0) &&
84  (0.5 <= x && x <= 1.0) &&
85  (0.0 <= y && y <= 0.5) &&
86  (0.0 <= z && z <= 0.5) )
87  return 1;
88 
89  // subtet #2 E2 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.5 <= s <= 1.0 && 0.0 <= t <= 0.5
90  if ( (0.5 <= xyz && xyz <= 1.0) &&
91  (0.0 <= x && x <= 0.5) &&
92  (0.5 <= y && y <= 1.0) &&
93  (0.0 <= z && z <= 0.5) )
94  return 2;
95 
96  // subtet #3 E3 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.5 <= t <= 1.0
97  if ( (0.5 <= xyz && xyz <= 1.0) &&
98  (0.0 <= x && x <= 0.5) &&
99  (0.0 <= y && y <= 0.5) &&
100  (0.5 <= z && z <= 1.0) )
101  return 3;
102 
103  // subtet #4 E4 := 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= r + t <= 1.0 && 0.0 <= r <= 0.5
104  if ( (0.0 <= yz && yz <= 0.5) &&
105  (0.5 <= xy && xy <= 1.0) &&
106  (0.5 <= xz && xz <= 1.0) &&
107  (0.0 <= x && x <= 0.5) )
108  return 4;
109 
110  // subtet #5 E5 := 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.5 <= r + t <= 1.0 && 0.75 <= r + s + t <= 1.0
111  if ( (0.5 <= xy && xy <= 1.0) &&
112  (0.5 <= yz && yz <= 1.0) &&
113  (0.5 <= xz && xz <= 1.0) &&
114  (0.75 <= xyz && xyz <= 1.0) )
115  return 5;
116 
117  // subtet #6 E6 := 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= t <= 0.5
118  if ( (0.5 <= yz && yz <= 1.0) &&
119  (0.0 <= xy && xy <= 0.5) &&
120  (0.5 <= xz && xz <= 1.0) &&
121  (0.0 <= z && z <= 0.5) )
122  return 6;
123 
124  // subtet #7 E7 := 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= s <= 0.25
125  if ( (0.0 <= yz && yz <= 0.5) &&
126  (0.0 <= xy && xy <= 0.5) &&
127  (0.5 <= xz && xz <= 1.0) &&
128  (0.0 <= y && y <= 0.25) )
129  return 7;
130 
131  // subtet #8 E8 := 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.0 <= t <= 0.25
132  if ( (0.0 <= xz && xz <= 0.5) &&
133  (0.0 <= yz && yz <= 0.5) &&
134  (0.5 <= xy && xy <= 1.0) &&
135  (0.0 <= z && z <= 0.25) )
136  return 8;
137 
138  // subtet #9 E9 := 0.0 <= r + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.0 <= s <= 0.5
139  if ( (0.0 <= xz && xz <= 0.5) &&
140  (0.5 <= xy && xy <= 1.0) &&
141  (0.5 <= yz && yz <= 1.0) &&
142  (0.0 <= y && y <= 0.5) )
143  return 9;
144 
145  // subtet #10 E10 := 0.0 <= r + t <= 0.5 && 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.0 <= r <= 0.25
146  if ( (0.0 <= xz && xz <= 0.5) &&
147  (0.5 <= yz && yz <= 1.0) &&
148  (0.0 <= xy && xy <= 0.5) &&
149  (0.0 <= x && x <= 0.25) )
150  return 10;
151 
152  // subtet #11 E11 := 0.5 <= r + s + t <= 0.75 && 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5
153  if ( (0.5 <= xyz && xyz <= 0.75) &&
154  (0.0 <= xz && xz <= 0.5) &&
155  (0.0 <= yz && yz <= 0.5) &&
156  (0.0 <= xy && xy <= 0.5) )
157  return 11;
158 
159  return -1;
160  }
161 
162  template<EOperator opType>
163  template<typename outputValueViewType,
164  typename inputPointViewType>
165  KOKKOS_INLINE_FUNCTION
166  void
167  Basis_HGRAD_TET_COMP12_FEM::Serial<opType>::
168  getValues( outputValueViewType output,
169  const inputPointViewType input ) {
170  switch (opType) {
171  case OPERATOR_VALUE: {
172  const typename inputPointViewType::value_type r = input(0);
173  const typename inputPointViewType::value_type s = input(1);
174  const typename inputPointViewType::value_type t = input(2);
175 
176  // initialize output
177  for (auto i=0;i<10;++i)
178  output.access(i) = 0.0;
179 
180  const auto subtet = getLocalSubTet( r, s, t );
181 
182  // idependent verification shows that a given point will produce
183  // the same shape functions for each tet that contains it, so
184  // we only need to use the first one returned.
185  if (subtet != -1) {
186  typename inputPointViewType::value_type aux = 0.0;
187  switch (subtet) {
188  case 0:
189  output.access(0) = 1. - 2. * (r + s + t);
190  output.access(4) = 2. * r;
191  output.access(6) = 2. * s;
192  output.access(7) = 2. * t;
193  break;
194  case 1:
195  output.access(1) = 2. * r - 1.;
196  output.access(4) = 2. - 2. * (r + s + t);
197  output.access(5) = 2. * s;
198  output.access(8) = 2. * t;
199  break;
200  case 2:
201  output.access(2) = 2. * s - 1.;
202  output.access(5) = 2. * r;
203  output.access(6) = 2. - 2. * (r + s + t);
204  output.access(9) = 2. * t;
205  break;
206  case 3:
207  output.access(3) = 2. * t - 1.;
208  output.access(7) = 2. - 2. * (r + s + t);
209  output.access(8) = 2. * r;
210  output.access(9) = 2. * s;
211  break;
212  case 4:
213  output.access(4) = 1. - 2. * (s + t);
214  output.access(5) = 2. * (r + s) - 1.;
215  output.access(8) = 2. * (r + t) - 1.;
216  aux = 2. - 4. * r;
217  break;
218  case 5:
219  output.access(5) = 2. * (r + s) - 1.;
220  output.access(8) = 2. * (r + t) - 1.;
221  output.access(9) = 2. * (s + t) - 1.;
222  aux = 4. - 4. * (r + s + t);
223  break;
224  case 6:
225  output.access(7) = 1. - 2. * (r + s);
226  output.access(8) = 2. * (r + t) - 1.;
227  output.access(9) = 2. * (s + t) - 1.;
228  aux = 2. - 4. * t;
229  break;
230  case 7:
231  output.access(4) = 1. - 2. * (s + t);
232  output.access(7) = 1. - 2. * (r + s);
233  output.access(8) = 2. * (r + t) - 1.;
234  aux = 4. * s;
235  break;
236  case 8:
237  output.access(4) = 1. - 2. * (s + t);
238  output.access(5) = 2. * (r + s) - 1.;
239  output.access(6) = 1. - 2. * (r + t);
240  aux = 4. * t;
241  break;
242  case 9:
243  output.access(5) = 2. * (r + s) - 1.;
244  output.access(6) = 1. - 2. * (r + t);
245  output.access(9) = 2. * (s + t) - 1.;
246  aux = 2. - 4. * s;
247  break;
248  case 10:
249  output.access(6) = 1. - 2. * (r + t);
250  output.access(7) = 1. - 2. * (r + s);
251  output.access(9) = 2. * (s + t) - 1.;
252  aux = 4. * r;
253  break;
254  case 11:
255  output.access(4) = 1. - 2. * (s + t);
256  output.access(6) = 1. - 2. * (r + t);
257  output.access(7) = 1. - 2. * (r + s);
258  aux = 4. * (r + s + t) - 2.;
259  break;
260  }
261  for (auto i=4;i<10;++i)
262  output.access(i) += aux/6.0;
263  }
264  break;
265  }
266  case OPERATOR_GRAD: {
267  const typename inputPointViewType::value_type r = input(0);
268  const typename inputPointViewType::value_type s = input(1);
269  const typename inputPointViewType::value_type t = input(2);
270 
271  output.access(0,0) = (-17 + 20*r + 20*s + 20*t)/8.;
272  output.access(0,1) = (-17 + 20*r + 20*s + 20*t)/8.;
273  output.access(0,2) = (-17 + 20*r + 20*s + 20*t)/8.;
274  output.access(1,0) = -0.375 + (5*r)/2.;
275  output.access(1,1) = 0.;
276  output.access(1,2) = 0.;
277  output.access(2,0) = 0.;
278  output.access(2,1) = -0.375 + (5*s)/2.;
279  output.access(2,2) = 0.;
280  output.access(3,0) = 0.;
281  output.access(3,1) = 0.;
282  output.access(3,2) = -0.375 + (5*t)/2.;
283  output.access(4,0) = (-35*(-1 + 2*r + s + t))/12.;
284  output.access(4,1) = (-4 - 35*r + 5*s + 10*t)/12.;
285  output.access(4,2) = (-4 - 35*r + 10*s + 5*t)/12.;
286  output.access(5,0) = (-1 + 5*r + 40*s - 5*t)/12.;
287  output.access(5,1) = (-1 + 40*r + 5*s - 5*t)/12.;
288  output.access(5,2) = (-5*(-1 + r + s + 2*t))/12.;
289  output.access(6,0) = (-4 + 5*r - 35*s + 10*t)/12.;
290  output.access(6,1) = (-35*(-1 + r + 2*s + t))/12.;
291  output.access(6,2) = (-4 + 10*r - 35*s + 5*t)/12.;
292  output.access(7,0) = (-4 + 5*r + 10*s - 35*t)/12.;
293  output.access(7,1) = (-4 + 10*r + 5*s - 35*t)/12.;
294  output.access(7,2) = (-35*(-1 + r + s + 2*t))/12.;
295  output.access(8,0) = (-1 + 5*r - 5*s + 40*t)/12.;
296  output.access(8,1) = (-5*(-1 + r + 2*s + t))/12.;
297  output.access(8,2) = (-1 + 40*r - 5*s + 5*t)/12.;
298  output.access(9,0) = (-5*(-1 + 2*r + s + t))/12.;
299  output.access(9,1) = (-1 - 5*r + 5*s + 40*t)/12.;
300  output.access(9,2) = (-1 - 5*r + 40*s + 5*t)/12.;
301  break;
302  }
303  case OPERATOR_MAX: {
304  const ordinal_type jend = output.extent(1);
305  const ordinal_type iend = output.extent(0);
306 
307  for (ordinal_type j=0;j<jend;++j)
308  for (auto i=0;i<iend;++i)
309  output.access(i, j) = 0.0;
310  break;
311  }
312  default: {
313  INTREPID2_TEST_FOR_ABORT( true ,
314  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
315  }
316  }
317  }
318 
319  template<typename SpT,
320  typename outputValueValueType, class ...outputValueProperties,
321  typename inputPointValueType, class ...inputPointProperties>
322  void
323  Basis_HGRAD_TET_COMP12_FEM::
324  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
325  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
326  const EOperator operatorType ) {
327  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
328  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
329  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
330 
331  // Number of evaluation points = dim 0 of inputPoints
332  const auto loopSize = inputPoints.extent(0);
333  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
334 
335  switch (operatorType) {
336 
337  case OPERATOR_VALUE: {
338  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
339  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
340  break;
341  }
342  case OPERATOR_GRAD:
343  case OPERATOR_D1: {
344  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
345  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
346  break;
347  }
348  case OPERATOR_CURL: {
349  INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_CURL, std::invalid_argument,
350  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
351  break;
352  }
353  case OPERATOR_DIV: {
354  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
355  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
356  break;
357  }
358  case OPERATOR_D2:
359  case OPERATOR_D3:
360  case OPERATOR_D4:
361  case OPERATOR_D5:
362  case OPERATOR_D6:
363  case OPERATOR_D7:
364  case OPERATOR_D8:
365  case OPERATOR_D9:
366  case OPERATOR_D10: {
367  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
368  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
369  break;
370  }
371  default: {
372  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
373  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Operator type not implemented" );
374  }
375  }
376  }
377 
378  }
379 
380  // -------------------------------------------------------------------------------------
381  template<typename SpT, typename OT, typename PT>
382  Basis_HGRAD_TET_COMP12_FEM<SpT,OT,PT>::
383  Basis_HGRAD_TET_COMP12_FEM() {
384  this->basisCardinality_ = 10;
385  this->basisDegree_ = 1;
386  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<11> >() );
387  this->basisType_ = BASIS_FEM_DEFAULT;
388  this->basisCoordinates_ = COORDINATES_CARTESIAN;
389 
390  {
391  // Basis-dependent intializations
392  const ordinal_type tagSize = 4; // size of DoF tag
393  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
394  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
395  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
396 
397  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
398  ordinal_type tags[] = { 0, 0, 0, 1,
399  0, 1, 0, 1,
400  0, 2, 0, 1,
401  0, 3, 0, 1,
402  1, 0, 0, 1,
403  1, 1, 0, 1,
404  1, 2, 0, 1,
405  1, 3, 0, 1,
406  1, 4, 0, 1,
407  1, 5, 0, 1 };
408 
409  // host view
410  ordinal_type_array_1d_host tagView(&tags[0],40);
411 
412  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
413  this->setOrdinalTagData(this->tagToOrdinal_,
414  this->ordinalToTag_,
415  tagView,
416  this->basisCardinality_,
417  tagSize,
418  posScDim,
419  posScOrd,
420  posDfOrd);
421  }
422 
423  // dofCoords on host and create its mirror view to device
424  Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
425  dofCoords("dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
426 
427  dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
428  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
429  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
430  dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
431  dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
432  dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
433  dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
434  dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
435  dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
436  dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
437 
438  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
439  Kokkos::deep_copy(this->dofCoords_, dofCoords);
440  }
441 }
442 
443 #endif
Intrepid2::Basis_HGRAD_TET_COMP12_FEM::Basis_HGRAD_TET_COMP12_FEM
Basis_HGRAD_TET_COMP12_FEM()
Constructor.
Definition: Intrepid2_HGRAD_TET_COMP12_FEMDef.hpp:383
Intrepid2::Impl::Basis_HGRAD_TET_COMP12_FEM::Serial
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.
Definition: Intrepid2_HGRAD_TET_COMP12_FEM.hpp:123
Intrepid2::Impl::Basis_HGRAD_TET_COMP12_FEM::getLocalSubTet
static KOKKOS_INLINE_FUNCTION ordinal_type getLocalSubTet(const pointValueType x, const pointValueType y, const pointValueType z)
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.
Definition: Intrepid2_HGRAD_TET_COMP12_FEMDef.hpp:63
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