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 
53 namespace 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)/8.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)/8.0;
80  output.access(1, 2) = 0.0;
81 
82  output.access(2, 0) = -(1.0 + y)*(1.0 - z)/8.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)/8.0;
88  output.access(3, 2) = 0.0;
89 
90  output.access(4, 0) = (1.0 - y)*(1.0 + z)/8.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)/8.0;
96  output.access(5, 2) = 0.0;
97 
98  output.access(6, 0) = -(1.0 + y)*(1.0 + z)/8.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)/8.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)/8.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)/8.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)/8.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)/8.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)/8.0;
131  output.access(0, 2) = (1.0 - z)/8.0;
132 
133  output.access(1, 0) = (1.0 + x)/8.0;
134  output.access(1, 1) = 0.0;
135  output.access(1, 2) = (1.0 - z)/8.0;
136 
137  output.access(2, 0) = 0.0;
138  output.access(2, 1) = (1.0 + y)/8.0;
139  output.access(2, 2) = (1.0 - z)/8.0;
140 
141  output.access(3, 0) = -(1.0 - x)/8.0;
142  output.access(3, 1) = 0.0;
143  output.access(3, 2) = (1.0 - z)/8.0;
144 
145  output.access(4, 0) = 0.0;
146  output.access(4, 1) = (1.0 - y)/8.0;
147  output.access(4, 2) = (1.0 + z)/8.0;
148 
149  output.access(5, 0) = -(1.0 + x)/8.0;
150  output.access(5, 1) = 0.0;
151  output.access(5, 2) = (1.0 + z)/8.0;
152 
153  output.access(6, 0) = 0.0;
154  output.access(6, 1) = -(1.0 + y)/8.0;
155  output.access(6, 2) = (1.0 + z)/8.0;
156 
157  output.access(7, 0) = (1.0 - x)/8.0;
158  output.access(7, 1) = 0.0;
159  output.access(7, 2) = (1.0 + z)/8.0;
160 
161  output.access(8, 0) = -(1.0 - x)/8.0;
162  output.access(8, 1) = (1.0 - y)/8.0;
163  output.access(8, 2) = 0.0;
164 
165  output.access(9, 0) = -(1.0 + x)/8.0;
166  output.access(9, 1) = -(1.0 - y)/8.0;
167  output.access(9, 2) = 0.0;
168 
169  output.access(10, 0) = (1.0 + x)/8.0;
170  output.access(10, 1) = -(1.0 + y)/8.0;
171  output.access(10, 2) = 0.0;
172 
173  output.access(11, 0) = (1.0 - x)/8.0;
174  output.access(11, 1) = (1.0 + y)/8.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 SpT,
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,SpT>::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 SpT, 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 SpT::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 SpT::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 SpT::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 (with tangents magnitude equal to edges' lengths)
347  dofCoeffs(0,0) = 2.0; dofCoeffs(0,1) = 0.0; dofCoeffs(0,2) = 0.0;
348  dofCoeffs(1,0) = 0.0; dofCoeffs(1,1) = 2.0; dofCoeffs(1,2) = 0.0;
349  dofCoeffs(2,0) = -2.0; dofCoeffs(2,1) = 0.0; dofCoeffs(2,2) = 0.0;
350  dofCoeffs(3,0) = 0.0; dofCoeffs(3,1) = -2.0; dofCoeffs(3,2) = 0.0;
351  dofCoeffs(4,0) = 2.0; dofCoeffs(4,1) = 0.0; dofCoeffs(4,2) = 0.0;
352  dofCoeffs(5,0) = 0.0; dofCoeffs(5,1) = 2.0; dofCoeffs(5,2) = 0.0;
353  dofCoeffs(6,0) = -2.0; dofCoeffs(6,1) = 0.0; dofCoeffs(6,2) = 0.0;
354  dofCoeffs(7,0) = 0.0; dofCoeffs(7,1) = -2.0; dofCoeffs(7,2) = 0.0;
355  dofCoeffs(8,0) = 0.0; dofCoeffs(8,1) = 0.0; dofCoeffs(8,2) = 2.0;
356  dofCoeffs(9,0) = 0.0; dofCoeffs(9,1) = 0.0; dofCoeffs(9,2) = 2.0;
357  dofCoeffs(10,0) = 0.0; dofCoeffs(10,1) = 0.0; dofCoeffs(10,2) = 2.0;
358  dofCoeffs(11,0) = 0.0; dofCoeffs(11,1) = 0.0; dofCoeffs(11,2) = 2.0;
359 
360  this->dofCoeffs_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoeffs);
361  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
362 
363  }
364 
365 }// namespace Intrepid2
366 
367 #endif
Kokkos::View< ordinal_type *, typename ExecSpaceType::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.