Intrepid2
Intrepid2_HGRAD_TRI_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_TRI_C2_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_TRI_C2_FEM_DEF_HPP__
51 
52 namespace 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_TRI_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 + y - 1.0)*(2.0*x + 2.0*y - 1.0);
73  output.access(1) = x*(2.0*x - 1.0);
74  output.access(2) = y*(2.0*y - 1.0);
75  output.access(3) = -4.0*x*(x + y - 1.0);
76  output.access(4) = 4.0*x*y;
77  output.access(5) = -4.0*y*(x + y - 1.0);
78  break;
79  }
80  case OPERATOR_D1:
81  case OPERATOR_GRAD: {
82  const auto x = input(0);
83  const auto y = input(1);
84  // output is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
85  output.access(0, 0) = 4.0*x + 4.0*y - 3.0;
86  output.access(0, 1) = 4.0*x + 4.0*y - 3.0;
87 
88  output.access(1, 0) = 4.0*x - 1.0;
89  output.access(1, 1) = 0.0;
90 
91  output.access(2, 0) = 0.0;
92  output.access(2, 1) = 4.0*y - 1.0;
93 
94  output.access(3, 0) = -4.0*(2.0*x + y - 1.0);
95  output.access(3, 1) = -4.0*x;
96 
97  output.access(4, 0) = 4.0*y;
98  output.access(4, 1) = 4.0*x;
99 
100  output.access(5, 0) = -4.0*y;
101  output.access(5, 1) = -4.0*(x + 2.0*y - 1.0);
102  break;
103  }
104  case OPERATOR_CURL: {
105  const auto x = input(0);
106  const auto y = input(1);
107  // CURL(u) = (u_y, -u_x), is rotated GRAD
108  output.access(0, 1) =-(4.0*x + 4.0*y - 3.0);
109  output.access(0, 0) = 4.0*x + 4.0*y - 3.0;
110 
111  output.access(1, 1) =-(4.0*x - 1.0);
112  output.access(1, 0) = 0.0;
113 
114  output.access(2, 1) = 0.0;
115  output.access(2, 0) = 4.0*y - 1.0;
116 
117  output.access(3, 1) = 4.0*(2.0*x + y - 1.0);
118  output.access(3, 0) = -4.0*x;
119 
120  output.access(4, 1) = -4.0*y;
121  output.access(4, 0) = 4.0*x;
122 
123  output.access(5, 1) = 4.0*y;
124  output.access(5, 0) = -4.0*(x + 2.0*y - 1.0);
125  break;
126  }
127  case OPERATOR_D2: {
128  // output is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
129  // D2 -> (2,0) -> dx^2.
130  output.access(0, 0) = 4.0;
131  output.access(1, 0) = 4.0;
132  output.access(2, 0) = 0.0;
133  output.access(3, 0) =-8.0;
134  output.access(4, 0) = 0.0;
135  output.access(5, 0) = 0.0;
136 
137  // D2 -> (1,1) -> dx dy
138  output.access(0, 1) = 4.0;
139  output.access(1, 1) = 0.0;
140  output.access(2, 1) = 0.0;
141  output.access(3, 1) =-4.0;
142  output.access(4, 1) = 4.0;
143  output.access(5, 1) =-4.0;
144 
145  // D2 -> (0,2) -> dy^2
146  output.access(0, 2) = 4.0;
147  output.access(1, 2) = 0.0;
148  output.access(2, 2) = 4.0;
149  output.access(3, 2) = 0.0;
150  output.access(4, 2) = 0.0;
151  output.access(5, 2) =-8.0;
152  break;
153  }
154  case OPERATOR_MAX: {
155  const ordinal_type jend = output.extent(1);
156  const ordinal_type iend = output.extent(0);
157 
158  for (ordinal_type j=0;j<jend;++j)
159  for (ordinal_type i=0;i<iend;++i)
160  output.access(i, j) = 0.0;
161  break;
162  }
163  default: {
164  INTREPID2_TEST_FOR_ABORT( opType != OPERATOR_VALUE &&
165  opType != OPERATOR_GRAD &&
166  opType != OPERATOR_CURL &&
167  opType != OPERATOR_D1 &&
168  opType != OPERATOR_D2 &&
169  opType != OPERATOR_MAX,
170  ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_C2_FEM::Serial::getValues) operator is not supported");
171  }
172  }
173  }
174 
175  template<typename DT,
176  typename outputValueValueType, class ...outputValueProperties,
177  typename inputPointValueType, class ...inputPointProperties>
178  void
179  Basis_HGRAD_TRI_C2_FEM::
180  getValues( const typename DT::execution_space& space,
181  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
182  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
183  const EOperator operatorType ) {
184  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
185  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
186  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
187 
188  // Number of evaluation points = dim 0 of inputPoints
189  const auto loopSize = inputPoints.extent(0);
190  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
191 
192  switch (operatorType) {
193 
194  case OPERATOR_VALUE: {
195  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
196  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
197  break;
198  }
199  case OPERATOR_GRAD:
200  case OPERATOR_D1: {
201  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
202  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
203  break;
204  }
205  case OPERATOR_CURL: {
206  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
207  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
208  break;
209  }
210  case OPERATOR_DIV: {
211  INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
212  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D.");
213  break;
214  }
215  case OPERATOR_D2: {
216  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
217  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
218  break;
219  }
220  case OPERATOR_D3:
221  case OPERATOR_D4:
222  case OPERATOR_D5:
223  case OPERATOR_D6:
224  case OPERATOR_D7:
225  case OPERATOR_D8:
226  case OPERATOR_D9:
227  case OPERATOR_D10: {
228  typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
229  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
230  break;
231  }
232  default: {
233  INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
234  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type");
235  }
236  }
237  }
238 
239  }
240  // -------------------------------------------------------------------------------------
241 
242 
243  template< typename DT, typename OT, typename PT>
244  Basis_HGRAD_TRI_C2_FEM<DT,OT,PT>::
245  Basis_HGRAD_TRI_C2_FEM() {
246  this->basisCardinality_ = 6;
247  this->basisDegree_ = 2;
248  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
249  this->basisType_ = BASIS_FEM_DEFAULT;
250  this->basisCoordinates_ = COORDINATES_CARTESIAN;
251  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
252 
253  // initialize tags
254  {
255  // Basis-dependent initializations
256  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
257  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
258  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
259  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
260 
261  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
262  ordinal_type tags[24] = { 0, 0, 0, 1,
263  0, 1, 0, 1,
264  0, 2, 0, 1,
265  1, 0, 0, 1,
266  1, 1, 0, 1,
267  1, 2, 0, 1};
268 
269  //host tags
270  OrdinalTypeArray1DHost tagView(&tags[0], 24);
271 
272  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
273  this->setOrdinalTagData(this->tagToOrdinal_,
274  this->ordinalToTag_,
275  tagView,
276  this->basisCardinality_,
277  tagSize,
278  posScDim,
279  posScOrd,
280  posDfOrd);
281  }
282 
283  // dofCoords on host and create its mirror view to device
284  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
285  dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
286 
287  dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0;
288  dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0;
289  dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0;
290  dofCoords(3,0) = 0.5; dofCoords(3,1) = 0.0;
291  dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.5;
292  dofCoords(5,0) = 0.0; dofCoords(5,1) = 0.5;
293 
294  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
295  Kokkos::deep_copy(this->dofCoords_, dofCoords);
296  }
297 
298 }// namespace Intrepid2
299 #endif
Intrepid2::Basis_HGRAD_TRI_C2_FEM::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:158
Intrepid2::Basis_HGRAD_TRI_C2_FEM::Basis_HGRAD_TRI_C2_FEM
Basis_HGRAD_TRI_C2_FEM()
Constructor.
Definition: Intrepid2_HGRAD_TRI_C2_FEMDef.hpp:245
Generated by   doxygen 1.8.5