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