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