Intrepid2
Intrepid2_HGRAD_LINE_Cn_FEM_JACOBIDef.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
50#ifndef __INTREPID2_HGRAD_LINE_CN_FEM_JACOBI_DEF_HPP__
51#define __INTREPID2_HGRAD_LINE_CN_FEM_JACOBI_DEF_HPP__
52
53namespace Intrepid2 {
54 // -------------------------------------------------------------------------------------
55
56 namespace Impl {
57
58 // output (N,P,D)
59 // input (P,D) - assumes that it has a set of points to amortize the function call cost for jacobi polynomial.
60 template<EOperator opType>
61 template<typename OutputViewType,
62 typename inputViewType>
63 KOKKOS_INLINE_FUNCTION
64 void
65 Basis_HGRAD_LINE_Cn_FEM_JACOBI::Serial<opType>::
66 getValues( OutputViewType output,
67 const inputViewType input,
68 const ordinal_type order,
69 const double alpha,
70 const double beta,
71 const ordinal_type operatorDn ) {
72 // cardinality of the evaluation order
73 const ordinal_type card = order + 1;
74 ordinal_type opDn = operatorDn;
75
76 const auto pts = Kokkos::subview( input, Kokkos::ALL(), 0 );
77 const ordinal_type np = input.extent(0);
78
79 switch (opType) {
80 case OPERATOR_VALUE: {
81 const Kokkos::View<typename inputViewType::value_type*,
82 typename inputViewType::memory_space,Kokkos::MemoryUnmanaged> null;
83 for (ordinal_type p=0;p<card;++p) {
84 auto poly = Kokkos::subview( output, p, Kokkos::ALL() );
85 Polylib::Serial::JacobiPolynomial(np, pts, poly, null, p, alpha, beta);
86 }
87 break;
88 }
89 case OPERATOR_GRAD:
90 case OPERATOR_D1: {
91 for (ordinal_type p=0;p<card;++p) {
92 auto polyd = Kokkos::subview( output, p, Kokkos::ALL(), 0 );
93 Polylib::Serial::JacobiPolynomialDerivative(np, pts, polyd, p, alpha, beta);
94 }
95 break;
96 }
97 case OPERATOR_D2:
98 case OPERATOR_D3:
99 case OPERATOR_D4:
100 case OPERATOR_D5:
101 case OPERATOR_D6:
102 case OPERATOR_D7:
103 case OPERATOR_D8:
104 case OPERATOR_D9:
105 case OPERATOR_D10:
106 opDn = getOperatorOrder(opType);
107 case OPERATOR_Dn: {
108 {
109 const ordinal_type pend = output.extent(0);
110 const ordinal_type iend = output.extent(1);
111 const ordinal_type jend = output.extent(2);
112
113 for (ordinal_type p=0;p<pend;++p)
114 for (ordinal_type i=0;i<iend;++i)
115 for (ordinal_type j=0;j<jend;++j)
116 output.access(p, i, j) = 0.0;
117 }
118 {
119 const Kokkos::View<typename inputViewType::value_type*,
120 typename inputViewType::memory_space,Kokkos::MemoryUnmanaged> null;
121
122 for (ordinal_type p=opDn;p<card;++p) {
123 double scaleFactor = 1.0;
124 for (ordinal_type i=1;i<=opDn;++i)
125 scaleFactor *= 0.5*(p + alpha + beta + i);
126
127 const auto poly = Kokkos::subview( output, p, Kokkos::ALL(), 0 );
128 Polylib::Serial::JacobiPolynomial(np, pts, poly, null, p-opDn, alpha+opDn, beta+opDn);
129 for (ordinal_type i=0;i<np;++i)
130 poly(i) = scaleFactor*poly(i);
131 }
132 }
133 break;
134 }
135 default: {
136 INTREPID2_TEST_FOR_ABORT( true,
137 ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM_JACOBI::Serial::getValues) operator is not supported");
138 }
139 }
140 }
141
142 // -------------------------------------------------------------------------------------
143
144 template<typename DT, ordinal_type numPtsPerEval,
145 typename outputValueValueType, class ...outputValueProperties,
146 typename inputPointValueType, class ...inputPointProperties>
147 void
148 Basis_HGRAD_LINE_Cn_FEM_JACOBI::
149 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
150 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
151 const ordinal_type order,
152 const double alpha,
153 const double beta,
154 const EOperator operatorType ) {
155 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
156 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
157 typedef typename DT::execution_space ExecSpaceType;
158
159 // loopSize corresponds to the # of points
160 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
161 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
162 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
163 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
164
165 switch (operatorType) {
166 case OPERATOR_VALUE: {
167 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
168 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
169 order, alpha, beta) );
170 break;
171 }
172 case OPERATOR_GRAD:
173 case OPERATOR_D1: {
174 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD,numPtsPerEval> FunctorType;
175 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
176 order, alpha, beta) );
177 break;
178 }
179 case OPERATOR_D2:
180 case OPERATOR_D3:
181 case OPERATOR_D4:
182 case OPERATOR_D5:
183 case OPERATOR_D6:
184 case OPERATOR_D7:
185 case OPERATOR_D8:
186 case OPERATOR_D9:
187 case OPERATOR_D10: {
188 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_Dn,numPtsPerEval> FunctorType;
189 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
190 order, alpha, beta,
191 getOperatorOrder(operatorType)) );
192 break;
193 }
194 case OPERATOR_DIV:
195 case OPERATOR_CURL: {
196 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV ||
197 operatorType == OPERATOR_CURL, std::invalid_argument,
198 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): invalid operator type (div and curl).");
199 break;
200 }
201 default: {
202 INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
203 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): invalid operator type");
204 }
205 }
206 }
207 }
208
209 // -------------------------------------------------------------------------------------
210
211 template<typename DT, typename OT, typename PT>
213 Basis_HGRAD_LINE_Cn_FEM_JACOBI( const ordinal_type order,
214 const double alpha,
215 const double beta ) {
216 this->basisCardinality_ = order+1;
217 this->basisDegree_ = order;
218 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
219 this->basisType_ = BASIS_FEM_HIERARCHICAL;
220 this->basisCoordinates_ = COORDINATES_CARTESIAN;
221 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
222
223 // jacobi
224 this->alpha_ = alpha;
225 this->beta_ = beta;
226
227 // initialize tags
228 {
229 // Basis-dependent intializations
230 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
231 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
232 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
233 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
234
235 ordinal_type tags[Parameters::MaxOrder+1][4];
236 const ordinal_type card = this->basisCardinality_;
237 for (ordinal_type i=0;i<card;++i) {
238 tags[i][0] = 1; // these are all "internal" i.e. "volume" DoFs
239 tags[i][1] = 0; // there is only one line
240 tags[i][2] = i; // local DoF id
241 tags[i][3] = card; // total number of DoFs
242 }
243
244 OrdinalTypeArray1DHost tagView(&tags[0][0], card*4);
245
246 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
247 // tags are constructed on host
248 this->setOrdinalTagData(this->tagToOrdinal_,
249 this->ordinalToTag_,
250 tagView,
251 this->basisCardinality_,
252 tagSize,
253 posScDim,
254 posScOrd,
255 posDfOrd);
256 }
257
258 // dof coords is not applicable to hierarchical functions
259 }
260
261
262}
263
264#endif
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Basis_HGRAD_LINE_Cn_FEM_JACOBI(const ordinal_type order, const double alpha=0, const double beta=0)
Constructor.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
static KOKKOS_INLINE_FUNCTION void JacobiPolynomialDerivative(const ordinal_type np, const zViewType z, polydViewType polyd, const ordinal_type n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
static KOKKOS_INLINE_FUNCTION void JacobiPolynomial(const ordinal_type np, const zViewType z, polyiViewType poly_in, polydViewType polyd, const ordinal_type n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .