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<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 DT,
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,typename DT::execution_space>::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 DT, typename OT, typename PT>
404 Basis_HGRAD_QUAD_C2_FEM<DT,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 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
412
413 {
414 // Basis-dependent intializations
415 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
416 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
417 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
418 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
419
420 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
421 ordinal_type tags[36] = { 0, 0, 0, 1,
422 0, 1, 0, 1,
423 0, 2, 0, 1,
424 0, 3, 0, 1,
425 // edge midpoints
426 1, 0, 0, 1,
427 1, 1, 0, 1,
428 1, 2, 0, 1,
429 1, 3, 0, 1,
430 // quad center
431 2, 0, 0, 1};
432
433 //host view
434 OrdinalTypeArray1DHost tagView(&tags[0],36);
435
436 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
437 this->setOrdinalTagData(this->tagToOrdinal_,
438 this->ordinalToTag_,
439 tagView,
440 this->basisCardinality_,
441 tagSize,
442 posScDim,
443 posScOrd,
444 posDfOrd);
445 }
446
447 // dofCoords on host and create its mirror view to device
448 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
449 dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
450
451 dofCoords(0,0) = -1.0; dofCoords(0,1) = -1.0;
452 dofCoords(1,0) = 1.0; dofCoords(1,1) = -1.0;
453 dofCoords(2,0) = 1.0; dofCoords(2,1) = 1.0;
454 dofCoords(3,0) = -1.0; dofCoords(3,1) = 1.0;
455
456 dofCoords(4,0) = 0.0; dofCoords(4,1) = -1.0;
457 dofCoords(5,0) = 1.0; dofCoords(5,1) = 0.0;
458 dofCoords(6,0) = 0.0; dofCoords(6,1) = 1.0;
459 dofCoords(7,0) = -1.0; dofCoords(7,1) = 0.0;
460
461 dofCoords(8,0) = 0.0; dofCoords(8,1) = 0.0;
462
463 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
464 Kokkos::deep_copy(this->dofCoords_, dofCoords);
465 }
466
467}// namespace Intrepid2
468#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< void, double, double >::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< void, double, double >::dofCoords_
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
Definition Intrepid2_Basis.hpp:318
Intrepid2::Impl::Basis_HGRAD_QUAD_C2_FEM::Functor
See Intrepid2::Basis_HGRAD_QUAD_C2_FEM.
Definition Intrepid2_HGRAD_QUAD_C2_FEM.hpp:129
Intrepid2::Impl::Basis_HGRAD_QUAD_C2_FEM::Serial
See Intrepid2::Basis_HGRAD_QUAD_C2_FEM.
Definition Intrepid2_HGRAD_QUAD_C2_FEM.hpp:105
Generated by doxygen 1.13.2