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
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_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( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
181 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
182 const EOperator operatorType ) {
183 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
184 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
185 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
186
187 // Number of evaluation points = dim 0 of inputPoints
188 const auto loopSize = inputPoints.extent(0);
189 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
190
191 switch (operatorType) {
192
193 case OPERATOR_VALUE: {
194 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE> FunctorType;
195 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
196 break;
197 }
198 case OPERATOR_GRAD:
199 case OPERATOR_D1: {
200 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD> FunctorType;
201 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
202 break;
203 }
204 case OPERATOR_CURL: {
205 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_CURL> FunctorType;
206 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
207 break;
208 }
209 case OPERATOR_DIV: {
210 INTREPID2_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
211 ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D.");
212 break;
213 }
214 case OPERATOR_D2: {
215 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_D2> FunctorType;
216 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
217 break;
218 }
219 case OPERATOR_D3:
220 case OPERATOR_D4:
221 case OPERATOR_D5:
222 case OPERATOR_D6:
223 case OPERATOR_D7:
224 case OPERATOR_D8:
225 case OPERATOR_D9:
226 case OPERATOR_D10: {
227 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_MAX> FunctorType;
228 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints) );
229 break;
230 }
231 default: {
232 INTREPID2_TEST_FOR_EXCEPTION( !( Intrepid2::isValidOperator(operatorType) ), std::invalid_argument,
233 ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type");
234 }
235 }
236 }
237
238 }
239 // -------------------------------------------------------------------------------------
240
241
242 template< typename DT, typename OT, typename PT>
245 this->basisCardinality_ = 6;
246 this->basisDegree_ = 2;
247 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
248 this->basisType_ = BASIS_FEM_DEFAULT;
249 this->basisCoordinates_ = COORDINATES_CARTESIAN;
250 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
251
252 // initialize tags
253 {
254 // Basis-dependent initializations
255 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
256 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
257 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
258 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
259
260 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
261 ordinal_type tags[24] = { 0, 0, 0, 1,
262 0, 1, 0, 1,
263 0, 2, 0, 1,
264 1, 0, 0, 1,
265 1, 1, 0, 1,
266 1, 2, 0, 1};
267
268 //host tags
269 OrdinalTypeArray1DHost tagView(&tags[0], 24);
270
271 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
272 this->setOrdinalTagData(this->tagToOrdinal_,
273 this->ordinalToTag_,
274 tagView,
275 this->basisCardinality_,
276 tagSize,
277 posScDim,
278 posScOrd,
279 posDfOrd);
280 }
281
282 // dofCoords on host and create its mirror view to device
283 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
284 dofCoords("dofCoordsHost", this->basisCardinality_,this->basisCellTopology_.getDimension());
285
286 dofCoords(0,0) = 0.0; dofCoords(0,1) = 0.0;
287 dofCoords(1,0) = 1.0; dofCoords(1,1) = 0.0;
288 dofCoords(2,0) = 0.0; dofCoords(2,1) = 1.0;
289 dofCoords(3,0) = 0.5; dofCoords(3,1) = 0.0;
290 dofCoords(4,0) = 0.5; dofCoords(4,1) = 0.5;
291 dofCoords(5,0) = 0.0; dofCoords(5,1) = 0.5;
292
293 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
294 Kokkos::deep_copy(this->dofCoords_, dofCoords);
295 }
296
297}// namespace Intrepid2
298#endif