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