50 #ifndef __INTREPID2_HGRAD_LINE_CN_FEM_JACOBI_DEF_HPP__
51 #define __INTREPID2_HGRAD_LINE_CN_FEM_JACOBI_DEF_HPP__
60 template<EOperator opType>
61 template<
typename OutputViewType,
62 typename inputViewType>
63 KOKKOS_INLINE_FUNCTION
65 Basis_HGRAD_LINE_Cn_FEM_JACOBI::Serial<opType>::
66 getValues( OutputViewType output,
67 const inputViewType input,
68 const ordinal_type order,
71 const ordinal_type operatorDn ) {
73 const ordinal_type card = order + 1;
74 ordinal_type opDn = operatorDn;
76 const auto pts = Kokkos::subview( input, Kokkos::ALL(), 0 );
77 const ordinal_type np = input.extent(0);
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() );
91 for (ordinal_type p=0;p<card;++p) {
92 auto polyd = Kokkos::subview( output, p, Kokkos::ALL(), 0 );
106 opDn = getOperatorOrder(opType);
109 const ordinal_type pend = output.extent(0);
110 const ordinal_type iend = output.extent(1);
111 const ordinal_type jend = output.extent(2);
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;
119 const Kokkos::View<
typename inputViewType::value_type*,
120 typename inputViewType::memory_space,Kokkos::MemoryUnmanaged> null;
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);
127 const auto poly = Kokkos::subview( output, p, Kokkos::ALL(), 0 );
129 for (ordinal_type i=0;i<np;++i)
130 poly(i) = scaleFactor*poly(i);
136 INTREPID2_TEST_FOR_ABORT(
true,
137 ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM_JACOBI::Serial::getValues) operator is not supported");
144 template<
typename DT, ordinal_type numPtsPerEval,
145 typename outputValueValueType,
class ...outputValueProperties,
146 typename inputPointValueType,
class ...inputPointProperties>
148 Basis_HGRAD_LINE_Cn_FEM_JACOBI::
149 getValues(
const typename DT::execution_space& space,
150 Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
151 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
152 const ordinal_type order,
155 const EOperator operatorType ) {
156 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
157 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
158 typedef typename DT::execution_space ExecSpaceType;
161 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
162 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
163 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
164 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
166 switch (operatorType) {
167 case OPERATOR_VALUE: {
168 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
169 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
170 order, alpha, beta) );
175 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_GRAD,numPtsPerEval> FunctorType;
176 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
177 order, alpha, beta) );
189 typedef Functor<outputValueViewType,inputPointViewType,OPERATOR_Dn,numPtsPerEval> FunctorType;
190 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints,
192 getOperatorOrder(operatorType)) );
196 case OPERATOR_CURL: {
197 INTREPID2_TEST_FOR_EXCEPTION( operatorType == OPERATOR_DIV ||
198 operatorType == OPERATOR_CURL, std::invalid_argument,
199 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): invalid operator type (div and curl).");
203 INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
204 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): invalid operator type");
212 template<
typename DT,
typename OT,
typename PT>
216 const double beta ) {
217 this->basisCardinality_ = order+1;
218 this->basisDegree_ = order;
219 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
220 this->basisType_ = BASIS_FEM_HIERARCHICAL;
221 this->basisCoordinates_ = COORDINATES_CARTESIAN;
222 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
225 this->alpha_ = alpha;
231 const ordinal_type tagSize = 4;
232 const ordinal_type posScDim = 0;
233 const ordinal_type posScOrd = 1;
234 const ordinal_type posDfOrd = 2;
237 const ordinal_type card = this->basisCardinality_;
238 for (ordinal_type i=0;i<card;++i) {
249 this->setOrdinalTagData(this->tagToOrdinal_,
252 this->basisCardinality_,
Basis_HGRAD_LINE_Cn_FEM_JACOBI(const ordinal_type order, const double alpha=0, const double beta=0)
Constructor.
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, .
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.