Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_LINE_Cn_FEM_JACOBIDef.hpp
Go to the documentation of this file.
00001 #ifndef INTREPID_HGRAD_LINE_CN_FEM_JACOBIDEF_HPP
00002 #define INTREPID_HGRAD_LINE_CN_FEM_JACOBIDEF_HPP
00003 // @HEADER
00004 // ************************************************************************
00005 //
00006 //                           Intrepid Package
00007 //                 Copyright (2007) Sandia Corporation
00008 //
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 //
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00040 //                    Denis Ridzal  (dridzal@sandia.gov), or
00041 //                    Kara Peterson (kjpeter@sandia.gov)
00042 //
00043 // ************************************************************************
00044 // @HEADER
00045 
00051 namespace Intrepid {
00052 
00053 
00054 template<class Scalar, class ArrayScalar>
00055 Basis_HGRAD_LINE_Cn_FEM_JACOBI<Scalar,ArrayScalar>::Basis_HGRAD_LINE_Cn_FEM_JACOBI(int order, Scalar alpha, Scalar beta) {
00056     this -> basisCardinality_  = order+1;
00057     this -> basisDegree_       = order;    
00058     this -> jacobiAlpha_       = alpha;    
00059     this -> jacobiBeta_        = beta;    
00060     this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
00061     this -> basisType_         = BASIS_FEM_HIERARCHICAL;
00062     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00063     this -> basisTagsAreSet_   = false;
00064 }
00065 
00066 
00067 
00068 template<class Scalar, class ArrayScalar> 
00069 void Basis_HGRAD_LINE_Cn_FEM_JACOBI<Scalar, ArrayScalar>::getValues(ArrayScalar &        outputValues,
00070                                                                     const ArrayScalar &  inputPoints,
00071                                                                     const EOperator      operatorType) const {
00072   
00073   // Verify arguments
00074 #ifdef HAVE_INTREPID_DEBUG
00075   Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
00076                                                       inputPoints,
00077                                                       operatorType,
00078                                                       this -> getBaseCellTopology(),
00079                                                       this -> getCardinality() );
00080 #endif
00081   
00082   // Number of evaluation points = dimension 0 of inputPoints
00083   int numPoints = inputPoints.dimension(0);  
00084   
00085   Teuchos::Array<Scalar> tmpPoints(numPoints);
00086   Teuchos::Array<Scalar> jacobiPolyAtPoints(numPoints);
00087 
00088   // Copy inputPoints into tmpPoints, to prepare for call to jacobfd
00089   for (int i=0; i<numPoints; i++) {
00090     tmpPoints[i] = inputPoints(i, 0);
00091   }
00092 
00093   try {
00094     switch (operatorType) {
00095     case OPERATOR_VALUE: {
00096       for (int ord = 0; ord < this -> basisCardinality_; ord++) {
00097         IntrepidPolylib::jacobfd(numPoints, &tmpPoints[0], &jacobiPolyAtPoints[0], (Scalar*)0, ord, jacobiAlpha_, jacobiBeta_);
00098         for (int pt = 0; pt < numPoints; pt++) {
00099           // outputValues is a rank-2 array with dimensions (basisCardinality_, numPoints)
00100           outputValues(ord, pt) = jacobiPolyAtPoints[pt];
00101         }
00102       }
00103     }
00104     break;
00105       
00106     case OPERATOR_GRAD:
00107     case OPERATOR_D1: {
00108       for (int ord = 0; ord < this -> basisCardinality_; ord++) {
00109         IntrepidPolylib::jacobd(numPoints, &tmpPoints[0], &jacobiPolyAtPoints[0], ord, jacobiAlpha_, jacobiBeta_);
00110         for (int pt = 0; pt < numPoints; pt++) {
00111           // outputValues is a rank-2 array with dimensions (basisCardinality_, numPoints)
00112           outputValues(ord, pt, 0) = jacobiPolyAtPoints[pt];
00113         }
00114       }
00115     }
00116     break;
00117 
00118     case OPERATOR_D2:
00119     case OPERATOR_D3:
00120     case OPERATOR_D4:
00121     case OPERATOR_D5:
00122     case OPERATOR_D6:
00123     case OPERATOR_D7:
00124     case OPERATOR_D8:
00125     case OPERATOR_D9:
00126     case OPERATOR_D10: {
00127       int d_order = getOperatorOrder( operatorType );
00128       // fill in derivatives of polynomials of degree 0 through d_order - 1  with 0
00129       // e.g. D2 annhialates linears.
00130       int stop_order;
00131       if (d_order > this->getDegree()) {
00132         stop_order = this->getDegree();
00133       }
00134       else {
00135         stop_order = d_order;
00136       }
00137       for (int p_order=0;p_order<stop_order;p_order++) {
00138         for (int pt=0;pt<numPoints;pt++) {
00139           outputValues(p_order,pt,0) = 0.0;
00140         }
00141       }
00142       // fill in rest of derivatives with the differentiation rule for Jacobi polynomials
00143       for (int p_order=d_order;p_order<=this->getDegree();p_order++) {
00144         // calculate the scaling factor with a little loop.
00145         Scalar scalefactor = 1.0;
00146         for (int d=1;d<=d_order;d++) {
00147           scalefactor *= 0.5 * ( p_order + jacobiAlpha_ + jacobiBeta_ + d );
00148         }
00149 
00150         // put in the right call to IntrepidPolyLib
00151         IntrepidPolylib::jacobfd(numPoints, &tmpPoints[0], 
00152                                  &jacobiPolyAtPoints[0], 
00153                                  (Scalar*)0, p_order-d_order, 
00154                                  jacobiAlpha_ + d_order, 
00155                                  jacobiBeta_ + d_order);
00156         for (int pt = 0; pt < numPoints; pt++) {
00157           // outputValues is a rank-3 array with dimensions (basisCardinality_, numPoints)
00158           outputValues(p_order, pt,0) = scalefactor *jacobiPolyAtPoints[pt];
00159         }
00160         
00161       }
00162       
00163     }
00164     break;
00165     case OPERATOR_DIV:
00166     case OPERATOR_CURL:
00167         TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00168                             ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Invalid operator type.");
00169       break;
00170     default:
00171       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
00172                           ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Invalid operator type.");
00173       break;
00174     }
00175   }
00176   catch (std::invalid_argument &exception){
00177     TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
00178                         ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): Operator failed");    
00179   }
00180 }
00181 
00182 
00183 
00184 template<class Scalar, class ArrayScalar>
00185 void Basis_HGRAD_LINE_Cn_FEM_JACOBI<Scalar, ArrayScalar>::getValues(ArrayScalar&           outputValues,
00186                                                                     const ArrayScalar &    inputPoints,
00187                                                                     const ArrayScalar &    cellVertices,
00188                                                                     const EOperator        operatorType) const {
00189   TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
00190                       ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM_JACOBI): FEM Basis calling an FVD member function");
00191 }
00192 
00193 
00194 
00195 template<class Scalar, class ArrayScalar>
00196 void Basis_HGRAD_LINE_Cn_FEM_JACOBI<Scalar, ArrayScalar>::setBasisParameters(int n, Scalar alpha, Scalar beta) {
00197   this -> basisCardinality_  = n+1;
00198   this -> basisDegree_       = n;
00199   this -> jacobiAlpha_       = alpha;
00200   this -> jacobiBeta_        = beta;
00201   this -> initializeTags();
00202 }
00203 
00204 
00205 
00206 template<class Scalar, class ArrayScalar>
00207 void Basis_HGRAD_LINE_Cn_FEM_JACOBI<Scalar, ArrayScalar>::initializeTags() {
00208 
00209   // Basis-dependent initializations
00210 
00211   int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00212   int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim
00213   int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00214   int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00215 
00216   FieldContainer<int> tags(this->basisCardinality_, 4);
00217 
00218   for (int i=0; i < this->basisCardinality_; i++) {
00219     tags(i, 0) = 1;                        // these are all "internal" i.e. "volume" DoFs
00220     tags(i, 1) = 0;                        // there is only one line
00221     tags(i, 2) = i;                        // local DoF id 
00222     tags(i, 3) = this->basisCardinality_;  // total number of DoFs 
00223   }
00224 
00225   // Basis-independent function, sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00226   Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00227                               this -> ordinalToTag_,
00228                               &tags[0],
00229                               this -> basisCardinality_,
00230                               tagSize,
00231                               posScDim,
00232                               posScOrd,
00233                               posDfOrd);
00234 }
00235 
00236 }// namespace Intrepid
00237 #endif