Intrepid
Intrepid_HGRAD_LINE_Cn_FEM_JACOBIDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_LINE_CN_FEM_JACOBIDEF_HPP
2 #define INTREPID_HGRAD_LINE_CN_FEM_JACOBIDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53 
54 template<class Scalar, class ArrayScalar>
56  this -> basisCardinality_ = order+1;
57  this -> basisDegree_ = order;
58  this -> jacobiAlpha_ = alpha;
59  this -> jacobiBeta_ = beta;
60  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
61  this -> basisType_ = BASIS_FEM_HIERARCHICAL;
62  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
63  this -> basisTagsAreSet_ = false;
64 }
65 
66 
67 
68 template<class Scalar, class ArrayScalar>
70  const ArrayScalar & inputPoints,
71  const EOperator operatorType) const {
72 
73  // Verify arguments
74 #ifdef HAVE_INTREPID_DEBUG
75  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
76  inputPoints,
77  operatorType,
78  this -> getBaseCellTopology(),
79  this -> getCardinality() );
80 #endif
81 
82  // Number of evaluation points = dimension 0 of inputPoints
83  int numPoints = inputPoints.dimension(0);
84 
85  Teuchos::Array<Scalar> tmpPoints(numPoints);
86  Teuchos::Array<Scalar> jacobiPolyAtPoints(numPoints);
87 
88  // Copy inputPoints into tmpPoints, to prepare for call to jacobfd
89  for (int i=0; i<numPoints; i++) {
90  tmpPoints[i] = inputPoints(i, 0);
91  }
92 
93  try {
94  switch (operatorType) {
95  case OPERATOR_VALUE: {
96  for (int ord = 0; ord < this -> basisCardinality_; ord++) {
97  IntrepidPolylib::jacobfd(numPoints, &tmpPoints[0], &jacobiPolyAtPoints[0], (Scalar*)0, ord, jacobiAlpha_, jacobiBeta_);
98  for (int pt = 0; pt < numPoints; pt++) {
99  // outputValues is a rank-2 array with dimensions (basisCardinality_, numPoints)
100  outputValues(ord, pt) = jacobiPolyAtPoints[pt];
101  }
102  }
103  }
104  break;
105 
106  case OPERATOR_GRAD:
107  case OPERATOR_D1: {
108  for (int ord = 0; ord < this -> basisCardinality_; ord++) {
109  IntrepidPolylib::jacobd(numPoints, &tmpPoints[0], &jacobiPolyAtPoints[0], ord, jacobiAlpha_, jacobiBeta_);
110  for (int pt = 0; pt < numPoints; pt++) {
111  // outputValues is a rank-2 array with dimensions (basisCardinality_, numPoints)
112  outputValues(ord, pt, 0) = jacobiPolyAtPoints[pt];
113  }
114  }
115  }
116  break;
117 
118  case OPERATOR_D2:
119  case OPERATOR_D3:
120  case OPERATOR_D4:
121  case OPERATOR_D5:
122  case OPERATOR_D6:
123  case OPERATOR_D7:
124  case OPERATOR_D8:
125  case OPERATOR_D9:
126  case OPERATOR_D10: {
127  int d_order = getOperatorOrder( operatorType );
128  // fill in derivatives of polynomials of degree 0 through d_order - 1 with 0
129  // e.g. D2 annhialates linears.
130  int stop_order;
131  if (d_order > this->getDegree()) {
132  stop_order = this->getDegree();
133  }
134  else {
135  stop_order = d_order;
136  }
137  for (int p_order=0;p_order<stop_order;p_order++) {
138  for (int pt=0;pt<numPoints;pt++) {
139  outputValues(p_order,pt,0) = 0.0;
140  }
141  }
142  // fill in rest of derivatives with the differentiation rule for Jacobi polynomials
143  for (int p_order=d_order;p_order<=this->getDegree();p_order++) {
144  // calculate the scaling factor with a little loop.
145  Scalar scalefactor = 1.0;
146  for (int d=1;d<=d_order;d++) {
147  scalefactor *= 0.5 * ( p_order + jacobiAlpha_ + jacobiBeta_ + d );
148  }
149 
150  // put in the right call to IntrepidPolyLib
151  IntrepidPolylib::jacobfd(numPoints, &tmpPoints[0],
152  &jacobiPolyAtPoints[0],
153  (Scalar*)0, p_order-d_order,
154  jacobiAlpha_ + d_order,
155  jacobiBeta_ + d_order);
156  for (int pt = 0; pt < numPoints; pt++) {
157  // outputValues is a rank-3 array with dimensions (basisCardinality_, numPoints)
158  outputValues(p_order, pt,0) = scalefactor *jacobiPolyAtPoints[pt];
159  }
160 
161  }
162 
163  }
164  break;
165  case OPERATOR_DIV:
166  case OPERATOR_CURL:
167  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
168  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Invalid operator type.");
169  break;
170  default:
171  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
172  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Invalid operator type.");
173  break;
174  }
175  }
176  catch (std::invalid_argument &exception){
177  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
178  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Operator failed");
179  }
180 }
181 
182 
183 
184 template<class Scalar, class ArrayScalar>
186  const ArrayScalar & inputPoints,
187  const ArrayScalar & cellVertices,
188  const EOperator operatorType) const {
189  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
190  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): FEM Basis calling an FVD member function");
191 }
192 
193 
194 
195 template<class Scalar, class ArrayScalar>
197  this -> basisCardinality_ = n+1;
198  this -> basisDegree_ = n;
199  this -> jacobiAlpha_ = alpha;
200  this -> jacobiBeta_ = beta;
201  this -> initializeTags();
202 }
203 
204 
205 
206 template<class Scalar, class ArrayScalar>
208 
209  // Basis-dependent initializations
210 
211  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
212  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
213  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
214  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
215 
216  FieldContainer<int> tags(this->basisCardinality_, 4);
217 
218  for (int i=0; i < this->basisCardinality_; i++) {
219  tags(i, 0) = 1; // these are all "internal" i.e. "volume" DoFs
220  tags(i, 1) = 0; // there is only one line
221  tags(i, 2) = i; // local DoF id
222  tags(i, 3) = this->basisCardinality_; // total number of DoFs
223  }
224 
225  // Basis-independent function, sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
226  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
227  this -> ordinalToTag_,
228  &tags[0],
229  this -> basisCardinality_,
230  tagSize,
231  posScDim,
232  posScOrd,
233  posDfOrd);
234 }
235 
236 }// namespace Intrepid
237 #endif
238 
239 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
240 #ifdef __GNUC__
241 #warning "The Intrepid package is deprecated"
242 #endif
243 #endif
244 
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
static void jacobd(const int np, const Scalar *z, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Calculate the derivative of Jacobi polynomials.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
static void jacobfd(const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
void setBasisParameters(int n, Scalar alpha=0, Scalar beta=0)
Sets private data basisDegree_, basisCardinality_, jacobiAlpha_, and jacobiBeta_, to n...
Basis_HGRAD_LINE_Cn_FEM_JACOBI(int order, Scalar alpha=0, Scalar beta=0)
Constructor.