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 
53 namespace 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 SpT, 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 ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType 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 SpT, 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 
222  // jacobi
223  this->alpha_ = alpha;
224  this->beta_ = beta;
225 
226  // initialize tags
227  {
228  // Basis-dependent intializations
229  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
230  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
231  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
232  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
233 
234  ordinal_type tags[Parameters::MaxOrder+1][4];
235  const ordinal_type card = this->basisCardinality_;
236  for (ordinal_type i=0;i<card;++i) {
237  tags[i][0] = 1; // these are all "internal" i.e. "volume" DoFs
238  tags[i][1] = 0; // there is only one line
239  tags[i][2] = i; // local DoF id
240  tags[i][3] = card; // total number of DoFs
241  }
242 
243  ordinal_type_array_1d_host tagView(&tags[0][0], card*4);
244 
245  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
246  // tags are constructed on host
247  this->setOrdinalTagData(this->tagToOrdinal_,
248  this->ordinalToTag_,
249  tagView,
250  this->basisCardinality_,
251  tagSize,
252  posScDim,
253  posScOrd,
254  posDfOrd);
255  }
256 
257  // dof coords is not applicable to hierarchical functions
258  }
259 
260 
261 }
262 
263 #endif
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, .
Basis_HGRAD_LINE_Cn_FEM_JACOBI(const ordinal_type order, const double alpha=0, const double beta=0)
Constructor.
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.