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
55namespace Intrepid {
56
57
58// No-arg constructor uses cubic Hermite interpolants based at the cell vertices
59template<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
79template<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
120template<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
159template<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
200template<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
221template<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
251template<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
376template<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
489template<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
538template<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
548template<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
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Implements Hermite interpolant basis of degree n on the reference Line cell. The basis has cardinalit...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
FieldContainer< Scalar > latticePts_
Holds the points defining the Hermite basis.
void setupVandermonde(bool factor=true)
Form the Legendre/Derivative Vandermonde matrix at the given lattice points and have the linear solve...
virtual void getDofCoords(ArrayScalar &DofCoords) const
implements the dofcoords interface
Basis_HGRAD_LINE_Hermite_FEM()
Default Constructor assumes the two interpolation points are the cell vertices. Cubic Hermite Interpo...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void recurrence(ArrayScalar &P, ArrayScalar &Px, const Scalar x) const
Evaluates and at a particular point .
void legendre_d(ArrayScalar &Pm, ArrayScalar &Pm1, const int m, const Scalar pt) const
Evaluates and at a particular point .
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....