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