Intrepid2
Intrepid2_HGRAD_TET_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_TET_C2_FEM_DEF_HPP__
50#define __INTREPID2_HGRAD_TET_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_TET_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 const auto z = input(2);
71
72 // output is a rank-2 array with dimensions (basisCardinality_, dim0)
73 output.access(0) = (-1. + x + y + z)*(-1. + 2.*x + 2.*y + 2.*z);
74 output.access(1) = x*(-1. + 2.*x);
75 output.access(2) = y*(-1. + 2.*y);
76 output.access(3) = z*(-1. + 2.*z);
77
78 output.access(4) = -4.*x*(-1. + x + y + z);
79 output.access(5) = 4.*x*y;
80 output.access(6) = -4.*y*(-1. + x + y + z);
81 output.access(7) = -4.*z*(-1. + x + y + z);
82 output.access(8) = 4.*x*z;
83 output.access(9) = 4.*y*z;
84 break;
85 }
86 case OPERATOR_D1:
87 case OPERATOR_GRAD: {
88 const auto x = input(0);
89 const auto y = input(1);
90 const auto z = input(2);
91
92 // output.access is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
93 output.access(0, 0) = -3.+ 4.*x + 4.*y + 4.*z;
94 output.access(0, 1) = -3.+ 4.*x + 4.*y + 4.*z;
95 output.access(0, 2) = -3.+ 4.*x + 4.*y + 4.*z;
96
97 output.access(1, 0) = -1.+ 4.*x;
98 output.access(1, 1) = 0.;
99 output.access(1, 2) = 0.;
100
101 output.access(2, 0) = 0.;
102 output.access(2, 1) = -1.+ 4.*y;
103 output.access(2, 2) = 0.;
104
105 output.access(3, 0) = 0.;
106 output.access(3, 1) = 0.;
107 output.access(3, 2) = -1.+ 4.*z;
108
109
110 output.access(4, 0) = -4.*(-1.+ 2*x + y + z);
111 output.access(4, 1) = -4.*x;
112 output.access(4, 2) = -4.*x;
113
114 output.access(5, 0) = 4.*y;
115 output.access(5, 1) = 4.*x;
116 output.access(5, 2) = 0.;
117
118 output.access(6, 0) = -4.*y;
119 output.access(6, 1) = -4.*(-1.+ x + 2*y + z);
120 output.access(6, 2) = -4.*y;
121
122 output.access(7, 0) = -4.*z;
123 output.access(7, 1) = -4.*z;
124 output.access(7, 2) = -4.*(-1.+ x + y + 2*z);
125
126 output.access(8, 0) = 4.*z;
127 output.access(8, 1) = 0.;
128 output.access(8, 2) = 4.*x;
129
130 output.access(9, 0) = 0.;
131 output.access(9, 1) = 4.*z;
132 output.access(9, 2) = 4.*y;
133 break;
134 }
135 case OPERATOR_D2: {
136 output.access(0, 0) = 4.;
137 output.access(0, 1) = 4.;
138 output.access(0, 2) = 4.;
139 output.access(0, 3) = 4.;
140 output.access(0, 4) = 4.;
141 output.access(0, 5) = 4.;
142
143 output.access(1, 0) = 4.;
144 output.access(1, 1) = 0.;
145 output.access(1, 2) = 0.;
146 output.access(1, 3) = 0.;
147 output.access(1, 4) = 0.;
148 output.access(1, 5) = 0.;
149
150 output.access(2, 0) = 0.;
151 output.access(2, 1) = 0.;
152 output.access(2, 2) = 0.;
153 output.access(2, 3) = 4.;
154 output.access(2, 4) = 0.;
155 output.access(2, 5) = 0.;
156
157 output.access(3, 0) = 0.;
158 output.access(3, 1) = 0.;
159 output.access(3, 2) = 0.;
160 output.access(3, 3) = 0.;
161 output.access(3, 4) = 0.;
162 output.access(3, 5) = 4.;
163
164 output.access(4, 0) = -8.;
165 output.access(4, 1) = -4.;
166 output.access(4, 2) = -4.;
167 output.access(4, 3) = 0.;
168 output.access(4, 4) = 0.;
169 output.access(4, 5) = 0.;
170
171 output.access(5, 0) = 0.;
172 output.access(5, 1) = 4.;
173 output.access(5, 2) = 0.;
174 output.access(5, 3) = 0.;
175 output.access(5, 4) = 0.;
176 output.access(5, 5) = 0.;
177
178 output.access(6, 0) = 0.;
179 output.access(6, 1) = -4.;
180 output.access(6, 2) = 0.;
181 output.access(6, 3) = -8.;
182 output.access(6, 4) = -4.;
183 output.access(6, 5) = 0;
184
185 output.access(7, 0) = 0.;
186 output.access(7, 1) = 0.;
187 output.access(7, 2) = -4.;
188 output.access(7, 3) = 0.;
189 output.access(7, 4) = -4.;
190 output.access(7, 5) = -8.;
191
192 output.access(8, 0) = 0.;
193 output.access(8, 1) = 0.;
194 output.access(8, 2) = 4.;
195 output.access(8, 3) = 0.;
196 output.access(8, 4) = 0.;
197 output.access(8, 5) = 0.;
198
199 output.access(9, 0) = 0.;
200 output.access(9, 1) = 0.;
201 output.access(9, 2) = 0.;
202 output.access(9, 3) = 0.;
203 output.access(9, 4) = 4.;
204 output.access(9, 5) = 0.;
205 break;
206 }
207 case OPERATOR_MAX: {
208 const ordinal_type jend = output.extent(1);
209 const ordinal_type iend = output.extent(0);
210
211 for (ordinal_type j=0;j<jend;++j)
212 for (ordinal_type i=0;i<iend;++i)
213 output.access(i, j) = 0.0;
214 break;
215 }
216 default: {
217 INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
218 opType != OPERATOR_GRAD &&
219 opType != OPERATOR_D1 &&
220 opType != OPERATOR_D2 &&
221 opType != OPERATOR_MAX,
222 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_C2_FEM::Serial::getValues) operator is not supported");
223 }
224 }
225 }
226
227 template<typename DT,
228 typename outputValueValueType, class ...outputValueProperties,
229 typename inputPointValueType, class ...inputPointProperties>
230 void
231 Basis_HGRAD_TET_C2_FEM::
232 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
233 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
234 const EOperator operatorType ) {
235 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
236 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
237 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
238
239 // Number of evaluation points = dim 0 of inputPoints
240 const auto loopSize = inputPoints.extent(0);
241 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
242
243 switch (operatorType) {
244
245 case OPERATOR_VALUE: {
246 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
247 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
248 break;
249 }
250 case OPERATOR_GRAD:
251 case OPERATOR_D1: {
252 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
253 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
254 break;
255 }
256 case OPERATOR_CURL: {
257 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
258 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
259 break;
260 }
261 case OPERATOR_DIV: {
262 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
263 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
264 break;
265 }
266 case OPERATOR_D2: {
267 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
268 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
269 break;
270 }
271 case OPERATOR_D3:
272 case OPERATOR_D4:
273 case OPERATOR_D5:
274 case OPERATOR_D6:
275 case OPERATOR_D7:
276 case OPERATOR_D8:
277 case OPERATOR_D9:
278 case OPERATOR_D10: {
279 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
280 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
281 break;
282 }
283 default: {
284 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
285 ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
286 }
287 }
288 }
289 }
290 // -------------------------------------------------------------------------------------
291
292 template<typename DT, typename OT, typename PT>
295 this->basisCardinality_ = 10;
296 this->basisDegree_ = 2;
297 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
298 this->basisType_ = BASIS_FEM_DEFAULT;
299 this->basisCoordinates_ = COORDINATES_CARTESIAN;
300 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
301
302 // intialize tags
303 {
304 // Basis-dependent intializations
305 const ordinal_type tagSize = 4; // size of DoF tag
306 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
307 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
308 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
309
310 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
311 ordinal_type tags[40] = { 0, 0, 0, 1,
312 0, 1, 0, 1,
313 0, 2, 0, 1,
314 0, 3, 0, 1,
315 1, 0, 0, 1,
316 1, 1, 0, 1,
317 1, 2, 0, 1,
318 1, 3, 0, 1,
319 1, 4, 0, 1,
320 1, 5, 0, 1,
321 };
322
323 // host tags
324 OrdinalTypeArray1DHost tagView(&tags[0], 40);
325
326 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
327 this->setOrdinalTagData(this->tagToOrdinal_,
328 this->ordinalToTag_,
329 tagView,
330 this->basisCardinality_,
331 tagSize,
332 posScDim,
333 posScOrd,
334 posDfOrd);
335 }
336
337 // dofCoords on host and create its mirror view to device
338 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
339 dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
340
341 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0; dofCoords(0,2) = 0.0;
342 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0; dofCoords(1,2) = 0.0;
343 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0; dofCoords(2,2) = 0.0;
344 dofCoords(3,0) = 0.0; dofCoords(3,1) = 0.0; dofCoords(3,2) = 1.0;
345
346 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.0; dofCoords(4,2) = 0.0;
347 dofCoords(5,0) = 0.5; dofCoords(5,1) = 0.5; dofCoords(5,2) = 0.0;
348 dofCoords(6,0) = 0.0; dofCoords(6,1) = 0.5; dofCoords(6,2) = 0.0;
349 dofCoords(7,0) = 0.0; dofCoords(7,1) = 0.0; dofCoords(7,2) = 0.5;
350 dofCoords(8,0) = 0.5; dofCoords(8,1) = 0.0; dofCoords(8,2) = 0.5;
351 dofCoords(9,0) = 0.0; dofCoords(9,1) = 0.5; dofCoords(9,2) = 0.5;
352
353 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
354 Kokkos::deep_copy(this->dofCoords_, dofCoords);
355 }
356
357}// namespace Intrepid2
358#endif