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
53namespace 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
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>
166 void
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 DT,
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,typename DT::execution_space>::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 DT, typename OT, typename PT>
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 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
390
391 {
392 // Basis-dependent intializations
393 const ordinal_type tagSize = 4; // size of DoF tag
394 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
395 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
396 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
397
398 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
399 ordinal_type tags[] = { 0, 0, 0, 1,
400 0, 1, 0, 1,
401 0, 2, 0, 1,
402 0, 3, 0, 1,
403 1, 0, 0, 1,
404 1, 1, 0, 1,
405 1, 2, 0, 1,
406 1, 3, 0, 1,
407 1, 4, 0, 1,
408 1, 5, 0, 1 };
409
410 // host view
411 OrdinalTypeArray1DHost tagView(&tags[0],40);
412
413 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
414 this->setOrdinalTagData(this->tagToOrdinal_,
415 this->ordinalToTag_,
416 tagView,
417 this->basisCardinality_,
418 tagSize,
419 posScDim,
420 posScOrd,
421 posDfOrd);
422 }
423
424 // dofCoords on host and create its mirror view to device
425 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
426 dofCoords("dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
427
428 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
429 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
430 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
431 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
432 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
433 dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
434 dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
435 dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
436 dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
437 dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
438
439 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
440 Kokkos::deep_copy(this->dofCoords_, dofCoords);
441 }
442}
443
444#endif
static KOKKOS_INLINE_FUNCTION ordinal_type getLocalSubTet(const pointValueType x, const pointValueType y, const pointValueType z)
See Intrepid2::Basis_HGRAD_TET_COMP12_FEM.