Intrepid
Intrepid_HGRAD_LINE_Hermite_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_LINE_HERMITE_FEMDEF_HPP
2 #define INTREPID_HGRAD_LINE_HERMITE_FEMDEF_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 #include<array>
52 #include<iostream>
53 #include<iomanip>
54 
55 namespace Intrepid {
56 
57 
58 // No-arg constructor uses cubic Hermite interpolants based at the cell vertices
59 template<class Scalar, class ArrayScalar>
61  latticePts_(2,1) {
62  this->basisCardinality_ = 4; // Four basis functions
63  this->basisDegree_ = 3; // Cubic polynomials
64  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
65  this->basisType_ = BASIS_FEM_DEFAULT;
66  this->basisCoordinates_ = COORDINATES_CARTESIAN;
67  this->basisTagsAreSet_ = false;
68 
69  latticePts_(0,0) = -1.0;
70  latticePts_(1,0) = 1.0;
71 
73 
74 } // no-arg constructor
75 
76 
77 
78 // Constructor with points as argument
79 template<class Scalar, class ArrayScalar>
81  latticePts_( pts.dimension(0), 1 ) {
82 
83  int n = pts.dimension(0);
84 
85  this->basisCardinality_ = 2*n;
86  this->basisDegree_ = 2*n-1;
87  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
88  this->basisType_ = BASIS_FEM_DEFAULT;
89  this->basisCoordinates_ = COORDINATES_CARTESIAN;
90  this->basisTagsAreSet_ = false;
91 
92  for( int i=0; i<n-1; ++i ) {
93  TEUCHOS_TEST_FOR_EXCEPTION( pts(i,0) >= pts(i+1,0), std::runtime_error ,
94  "Intrepid::Basis_HGRAD_LINE_Hermite_FEM Illegal points given to constructor" );
95  }
96 
97  // copy points int latticePts, correcting endpoints if needed
98  if (std::abs(pts(0,0)+1.0) < INTREPID_TOL) {
99  latticePts_(0,0) = -1.0;
100  }
101  else {
102  latticePts_(0,0) = pts(0,0);
103  }
104  for (int i=1;i<n-1;i++) {
105  latticePts_(i,0) = pts(i,0);
106  }
107  if (std::abs(pts(n-1,0)-1.0) < INTREPID_TOL) {
108  latticePts_(n-1,0) = 1.0;
109  }
110  else {
111  latticePts_(n-1,0) = pts(n-1,0);
112  }
113 
115 
116 } // Constructor with points given
117 
118 
119 // Constructor with point type as argument
120 template<class Scalar, class ArrayScalar>
122  const EPointType &pointType ) :
123  latticePts_(n,1) {
124 
125  TEUCHOS_TEST_FOR_EXCEPTION(n<2,std::invalid_argument,"Intrepid::Basis_HGRAD_LINE_Hermite_FEM requires the "
126  "number of interpolation points to be at least 2");
127 
128  this->basisCardinality_ = 2*n;
129  this->basisDegree_ = 2*n-1;
130  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
131  this->basisType_ = BASIS_FEM_DEFAULT;
132  this->basisCoordinates_ = COORDINATES_CARTESIAN;
133  this->basisTagsAreSet_ = false;
134 
135  switch(pointType) {
136  case POINTTYPE_EQUISPACED:
137  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ ,
138  this->basisCellTopology_ , n-1 , 0 , POINTTYPE_EQUISPACED );
139  break;
140  case POINTTYPE_SPECTRAL:
141  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ ,
142  this->basisCellTopology_ , n-1 , 0 , POINTTYPE_WARPBLEND );
143  break;
144  case POINTTYPE_SPECTRAL_OPEN:
145  PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( latticePts_ , n-1 );
146  break;
147  default:
148  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
149  "Basis_HGRAD_LINE_Hermite_FEM:: invalid point type" );
150  break;
151  }
152 
154 
155 } // Constructor with point type given
156 
157 
158 
159 template<class Scalar, class ArrayScalar>
161 
162  initializeTags();
163 
164  int nBf = this->getCardinality();
165  int n = nBf/2;
166 
167  V_.shape(nBf,nBf);
168 
169  // Make containers to store the Legendre polynomials and their derivatives
170  // at a given point
171  ArrayScalar P ( nBf );
172  ArrayScalar Px( nBf );
173 
174  // Loop over grid points
175  for( int i=0; i<n; ++i ) {
176 
177  recurrence(P,Px,latticePts_(i,0));
178 
179  // Loop over basis functions
180  for(int j=0; j<nBf; ++j ) {
181  V_(j, i ) = P (j);
182  V_(j, i+n) = Px(j);
183  }
184  }
185 
186  solver_.setMatrix(Teuchos::rcpFromRef(V_));
187 
188  if(factor) {
189  solver_.factorWithEquilibration(true);
190  solver_.factor();
191  isFactored_ = true;
192  }
193  else {
194  isFactored_ = false;
195  }
196 
197 }
198 
199 
200 template<class Scalar, class ArrayScalar>
202  ArrayScalar &Px,
203  const Scalar x ) const {
204 
205  int n = P.dimension(0);
206  Scalar q = x*x-1.0;
207 
208  P (0) = 1.0;
209  Px(0) = 0.0;
210 
211  // Loop over basis indices
212  for( int j=0; j<n-1; ++j ) {
213  P (j+1) = x*P(j) + q*Px(j)/(j+1); // Compute \f$P_{j+1}(x_i)\f$
214  Px(j+1) = (j+1)*P(j) + x*Px(j); // Compute \f$P'_{j+1}(x_i)\f$
215  }
216 
217 } // recurrence()
218 
219 
220 // Computes the derivatives of Legendre polynomials
221 template<class Scalar, class ArrayScalar>
223  ArrayScalar &Px,
224  const int m,
225  const Scalar x ) const {
226  // Compute P,P'
227  recurrence(P,Px,x);
228 
229  int C = this->getCardinality();
230 
231  // Loop over derivative orders
232  for( int k=1;k<m;++k) {
233  P = Px;
234 
235  // Loop over polynomial indices
236  for( int j=0; j<C; ++j ) {
237 
238  if( j<k ) {
239  Px(j) = 0;
240  }
241 
242  else {
243  Px(j) = (j+k)*P(j-1) + x*Px(j-1);
244  }
245  }
246  }
247 
248 } // legendre_d()
249 
250 
251 template<class Scalar, class ArrayScalar>
253  const ArrayScalar & inputPoints,
254  const EOperator operatorType) const {
255 
256  // Verify arguments
257 #ifdef HAVE_INTREPID_DEBUG
258  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
259  inputPoints,
260  operatorType,
261  this -> getBaseCellTopology(),
262  this -> getCardinality() );
263 #endif
264  // Number of evaluation points = dim 0 of inputPoints
265  int nPts = inputPoints.dimension(0);
266  int nBf = this->getCardinality();
267 
268  int n = nBf/2;
269 
270  // Legendre polynomials and their derivatives evaluated on inputPoints
271  SerialDenseMatrix legendre(nBf,nPts);
272 
273  // Hermite interpolants evaluated on inputPoints
274  SerialDenseMatrix hermite(nBf,nPts);
275 
276  ArrayScalar P (nBf);
277  ArrayScalar Px(nBf);
278 
279  int derivative_order;
280  int derivative_case = static_cast<int>(operatorType);
281 
282  if( derivative_case == 0 ) {
283  derivative_order = 0;
284  }
285  else if( derivative_case > 0 && derivative_case < 5 ) {
286  derivative_order = 1;
287  }
288  else {
289  derivative_order = derivative_case - 3;
290  }
291 
292  try {
293  // GRAD,CURL,DIV, and D1 are all the first derivative
294  switch (operatorType) {
295  case OPERATOR_VALUE:
296  {
297  for( int i=0; i<nPts; ++i ) {
298  recurrence( P, Px, inputPoints(i,0) );
299  for( int j=0; j<nBf; ++j ) {
300  legendre(j,i) = P(j);
301  }
302  }
303  break;
304  }
305  case OPERATOR_GRAD:
306  case OPERATOR_DIV:
307  case OPERATOR_CURL:
308  case OPERATOR_D1:
309  {
310  for( int i=0; i<nPts; ++i ) {
311  recurrence( P, Px, inputPoints(i,0) );
312  for( int j=0; j<nBf; ++j ) {
313  legendre(j,i) = Px(j);
314  }
315  }
316  break;
317  }
318  case OPERATOR_D2:
319  case OPERATOR_D3:
320  case OPERATOR_D4:
321  case OPERATOR_D5:
322  case OPERATOR_D6:
323  case OPERATOR_D7:
324  case OPERATOR_D8:
325  case OPERATOR_D9:
326  case OPERATOR_D10:
327  {
328  for( int i=0; i<nPts; ++i ) {
329  legendre_d( P, Px, derivative_order, inputPoints(i,0));
330  for( int j=0; j<nBf; ++j ) {
331  legendre(j,i) = Px(j);
332  }
333  }
334  break;
335  }
336  default:
337  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
338  ">>> ERROR (Basis_HGRAD_LINE_Hermite_FEM): Invalid operator type");
339 
340  } // switch(operatorType)
341  }
342  catch (std::invalid_argument &exception){
343  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
344  ">>> ERROR (Basis_HGRAD_LINE_Hermite_FEM): Operator failed");
345  }
346 
347  if( !isFactored_ ) {
348  solver_.factorWithEquilibration(true);
349  solver_.factor();
350  }
351 
352  solver_.setVectors(Teuchos::rcpFromRef(hermite),Teuchos::rcpFromRef(legendre));
353  solver_.solve();
354 
355  if(derivative_order > 0)
356  {
357  for( int i=0; i<n; ++i ) {
358  for( int j=0; j<nPts; ++j ) {
359  outputValues(2*i, j,0) = hermite(i, j);
360  outputValues(2*i+1,j,0) = hermite(i+n,j);
361  }
362  }
363  }
364  else {
365  for( int i=0; i<n; ++i ) {
366  for( int j=0; j<nPts; ++j ) {
367  outputValues(2*i ,j) = hermite(i, j);
368  outputValues(2*i+1,j) = hermite(i+n,j);
369  }
370  }
371  }
372 
373 } // getValues()
374 
375 
376 template<class Scalar, class ArrayScalar>
378 
379  // Basis-dependent intializations
380  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
381  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
382  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
383  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
384 
385  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
386 
387  int C = this->getCardinality();
388  tags_.reserve( tagSize * C );
389 
390  int n = C/2;
391 
392  int hasLeftVertex = static_cast<int>( latticePts_(0 , 0) == -1.0 );
393  int hasRightVertex = static_cast<int>( latticePts_(n-1, 0) == 1.0 );
394 
395  int internal_dof = C - 2*(hasLeftVertex+hasRightVertex);
396 
397  if( hasLeftVertex ) {
398 
399  // Value interpolant
400  tags_[0] = 0; // this is a vertex (interval end point)
401  tags_[1] = 0; // this is the first subcell
402  tags_[2] = 0; // this is the first DoF for this vertex
403  tags_[3] = 2; // this vertex has 2 DoF
404 
405  // Derivative interpolant
406  tags_[4] = 0; // this is a vertex (interval end point)
407  tags_[5] = 0; // this is the first subcell
408  tags_[6] = 1; // this is the second DoF for this vertex
409  tags_[7] = 2; // this vertex has 2 DoF
410 
411  }
412  else { // no left vertex
413 
414  // Value interpolant
415  tags_[0] = 1; // this point is on a line
416  tags_[1] = 0; // no subcells
417  tags_[2] = 0; // this is the first DoF for this line
418  tags_[3] = C-2*hasRightVertex; // this cell has 2n DoF
419 
420  // Derivative interpolant
421  tags_[4] = 1; // this point is on a line
422  tags_[5] = 0; // no subcells
423  tags_[6] = 1; // this is the second DoF for this line
424  tags_[7] = C-2*hasRightVertex; // this cell has 2n DoF
425  }
426 
427  if( hasRightVertex ) {
428  int i0 = C-2;
429  int i1 = C-1;
430 
431  // Value interpolant
432  tags_[4*i0 ] = 0;
433  tags_[4*i0+1] = hasLeftVertex;
434  tags_[4*i0+2] = 0;
435  tags_[4*i0+3] = 2;
436 
437  // Derivative interpolant
438  tags_[4*i1 ] = 0;
439  tags_[4*i1+1] = hasLeftVertex;
440  tags_[4*i1+2] = 1;
441  tags_[4*i1+3] = 2;
442  }
443  else { // no right vertex
444  int i0 = C-2;
445  int i1 = C-1;
446 
447  // Value interpolant
448  tags_[4*i0 ] = 1;
449  tags_[4*i0+1] = 0;
450  tags_[4*i0+2] = internal_dof-2;
451  tags_[4*i0+3] = internal_dof;
452 
453  // Derivative interpolant
454  tags_[4*i1 ] = 1;
455  tags_[4*i1+1] = 0;
456  tags_[4*i1+2] = internal_dof-1;
457  tags_[4*i1+3] = internal_dof;
458  }
459 
460  for( int i=1; i<n-1; ++i ) {
461  int i0 = 2*i; int i1 = 2*i+1;
462 
463  // Value interpolant
464  tags_[4*i0 ] = 1; // Points on a line (1 dimensional)
465  tags_[4*i0+1] = 0;
466  tags_[4*i0+2] = i0 - 2*hasLeftVertex;
467  tags_[4*i0+3] = internal_dof;
468 
469  // Derivative interpolant
470  tags_[4*i1 ] = 1;
471  tags_[4*i1+1] = 0;
472  tags_[4*i1+2] = i1 - 2*hasLeftVertex;
473  tags_[4*i1+3] = internal_dof;
474  }
475 
476  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
477  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
478  this -> ordinalToTag_,
479  tags_.data(),
480  this -> basisCardinality_,
481  tagSize,
482  posScDim,
483  posScOrd,
484  posDfOrd);
485 
486 }
487 
488 
489 template<class Scalar, class ArrayScalar>
491 
492  int nBf = this->getCardinality();
493 
494  os << "Tags:" << std::endl;
495  os << "-----" << std::endl;
496 
497 
498  os << "Index: ";
499  for( int i=0; i<nBf; ++i ) {
500  os << std::setw(4) << i;
501  }
502  os << std::endl;
503 
504 
505  os << "Subcell dim: ";
506  for( int i=0; i<nBf; ++i ) {
507  os << std::setw(4) << tags_[4*i];
508  }
509  os << std::endl;
510 
511 
512  os << "Subcell ord: ";
513  for( int i=0; i<nBf; ++i ) {
514  os << std::setw(4) << tags_[4*i+1];
515  }
516  os << std::endl;
517 
518 
519  os << "Subcell DoF: ";
520  for( int i=0; i<nBf; ++i ) {
521  os << std::setw(4) << tags_[4*i+2];
522  }
523  os << std::endl;
524 
525 
526  os << "Total Sc DoF: ";
527  for( int i=0; i<nBf; ++i ) {
528  os << std::setw(4) << tags_[4*i+3];
529  }
530  os << std::endl;
531  os << std::endl;
532 
533 }
534 
535 
536 
537 
538 template<class Scalar, class ArrayScalar>
540  const ArrayScalar & inputPoints,
541  const ArrayScalar & cellVertices,
542  const EOperator operatorType) const {
543  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
544  ">>> ERROR (Basis_HGRAD_LINE_Hermite_FEM): FEM Basis calling an FVD member function");
545 }
546 
547 
548 template<class Scalar, class ArrayScalar>
550 {
551  for (int i=0;i<latticePts_.dimension(0);i++)
552  {
553  for (int j=0;j<latticePts_.dimension(1);j++)
554  {
555  dofCoords(i,j) = latticePts_(i,j);
556  }
557  }
558  return;
559 }
560 
561 
562 
563 
564 
565 }// namespace Intrepid
566 #endif
567 
568 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
569 #ifdef __GNUC__
570 #warning "The Intrepid package is deprecated"
571 #endif
572 #endif
573 
virtual void getDofCoords(ArrayScalar &DofCoords) const
implements the dofcoords interface
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
FieldContainer< Scalar > latticePts_
Holds the points defining the Hermite basis.
EBasis basisType_
Type of the basis.
void legendre_d(ArrayScalar &Pm, ArrayScalar &Pm1, const int m, const Scalar pt) const
Evaluates and at a particular point .
void recurrence(ArrayScalar &P, ArrayScalar &Px, const Scalar x) const
Evaluates and at a particular point .
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
void setupVandermonde(bool factor=true)
Form the Legendre/Derivative Vandermonde matrix at the given lattice points and have the linear solve...
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
Basis_HGRAD_LINE_Hermite_FEM()
Default Constructor assumes the two interpolation points are the cell vertices. Cubic Hermite Interpo...