Intrepid: /usr/src/RPM/BUILD/trilinos-11.12.1/packages/intrepid/src/Discretization/Basis/Intrepid_HGRAD_POLY_C1

00001 #include "Intrepid_FieldContainer.hpp"
00002 #include "Sacado_Fad_DFad.hpp"
00003 
00004 
00005 namespace Intrepid{
00006 
00007   template<class Scalar, class ArrayScalar>
00008   Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::Basis_HGRAD_POLY_C1_FEM(const shards::CellTopology& cellTopology){
00009     this -> basisCardinality_  = cellTopology.getNodeCount();
00010     this -> basisDegree_       = 1;
00011     this -> basisCellTopology_ = cellTopology;
00012     this -> basisType_         = BASIS_FEM_DEFAULT;
00013     this -> basisCoordinates_  = COORDINATES_CARTESIAN;
00014     this -> basisTagsAreSet_   = false;
00015   }
00016 
00017   
00018   template<class Scalar, class ArrayScalar>
00019   void Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::initializeTags(){
00020     // Basis-dependent initializations
00021     int tagSize  = 4;        // size of DoF tag, i.e., number of fields in the tag
00022     int posScDim = 0;        // position in the tag, counting from 0, of the subcell dim 
00023     int posScOrd = 1;        // position in the tag, counting from 0, of the subcell ordinal
00024     int posDfOrd = 2;        // position in the tag, counting from 0, of DoF ordinal relative to the subcell
00025   
00026     // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 
00027     
00028     int *tags = new int[tagSize * this->getCardinality()];
00029     for (int i=0;i<this->getCardinality();i++){
00030       tags[4*i] = 0;
00031       tags[4*i+1] = i;
00032       tags[4*i+2] = 0;
00033       tags[4*i+3] = 1;
00034     }
00035     
00036     
00037     // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
00038     Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
00039                                 this -> ordinalToTag_,
00040                                 tags,
00041                                 this -> basisCardinality_,
00042                                 tagSize,
00043                                 posScDim,
00044                                 posScOrd,
00045                                 posDfOrd);
00046 
00047     delete [] tags;
00048   }  
00049   
00050   // this is the FEM reference element version, this should not be called for polygonal basis
00051   // since polygonal basis is defined on physical element.
00052   template<class Scalar, class ArrayScalar>
00053   void Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
00054                                                                const ArrayScalar& inputPoints,
00055                                                                const EOperator operatorType) const{
00056     TEUCHOS_TEST_FOR_EXCEPTION ( true, std::logic_error,
00057                          ">>>ERROR (Basis_HGRAD_POLY_C1_FEM): Polygonal basis calling FEM member function");
00058   }
00059 
00060 
00061   template<class Scalar, class ArrayScalar>
00062   void Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
00063                                                                const ArrayScalar& inputPoints,
00064                                                                const ArrayScalar& cellVertices,
00065                                                                const EOperator operatorType) const{
00066     
00067     
00068     // implement wachspress basis functions
00069     switch (operatorType) {
00070     case OPERATOR_VALUE:
00071       {
00072         shapeFunctions<Scalar,ArrayScalar>(outputValues,inputPoints,cellVertices);
00073       }
00074       break;
00075     case OPERATOR_GRAD:
00076       {
00077         FieldContainer<Sacado::Fad::DFad<Scalar> > dInput(inputPoints.dimension(0),inputPoints.dimension(1));
00078         for (int i=0;i<inputPoints.dimension(0);i++){
00079           for (int j=0;j<2;j++){
00080             dInput(i,j) = Sacado::Fad::DFad<Scalar>( inputPoints(i,j));
00081             dInput(i,j).diff(j,2);
00082           }
00083         }
00084         FieldContainer<Sacado::Fad::DFad<Scalar> > cellVerticesFad(cellVertices.dimension(0),cellVertices.dimension(1));
00085         for (int i=0;i<cellVertices.dimension(0);i++){
00086           for (int j=0;j<cellVertices.dimension(1);j++){
00087             cellVerticesFad(i,j) = Sacado::Fad::DFad<Scalar>( cellVertices(i,j) );
00088           }
00089         }
00090         
00091         FieldContainer<Sacado::Fad::DFad<Scalar> > dOutput(outputValues.dimension(0),outputValues.dimension(1));
00092         
00093         shapeFunctions<Sacado::Fad::DFad<Scalar>, FieldContainer<Sacado::Fad::DFad<Scalar> > >(dOutput,dInput,cellVerticesFad);
00094         
00095         for (int i=0;i<outputValues.dimension(0);i++){
00096           for (int j=0;j<outputValues.dimension(1);j++){
00097             for (int k=0;k<outputValues.dimension(2);k++){
00098               outputValues(i,j,k) = dOutput(i,j).dx(k);
00099             }
00100           }
00101         }
00102       }
00103       break;
00104     case OPERATOR_D1:
00105     case OPERATOR_D2:
00106     case OPERATOR_D3:
00107     case OPERATOR_D4:
00108     case OPERATOR_D5:
00109     case OPERATOR_D6:
00110     case OPERATOR_D7:
00111     case OPERATOR_D8:
00112     case OPERATOR_D9:
00113     case OPERATOR_D10:
00114       {
00115         TEUCHOS_TEST_FOR_EXCEPTION ( true, std::invalid_argument, 
00116                              ">>> ERROR (Basis_HGRAD_POLY_C1_FEM): operator not implemented yet");
00117       }
00118       break;
00119     default:
00120       TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType)), std::invalid_argument,
00121                           ">>> ERROR (Basis_HGRAD_POLY_C1_FEM): Invalid operator type");
00122       break;
00123     }
00124 
00125   }
00126 
00127   
00128   template<class Scalar, class ArrayScalar>
00129   template<class Scalar1, class ArrayScalar1>
00130   void Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::shapeFunctions(ArrayScalar1& outputValues,
00131                                                                     const ArrayScalar1& inputPoints,
00132                                                                     const ArrayScalar1& cellVertices)const{
00133     int numPoints = inputPoints.dimension(0);
00134     FieldContainer<Scalar1> weightFunctions( this->basisCardinality_,numPoints );
00135     evaluateWeightFunctions<Scalar1, ArrayScalar1>(weightFunctions,inputPoints,cellVertices);
00136     for (int pointIndex = 0;pointIndex<outputValues.dimension(1);pointIndex++){
00137       Scalar1 denominator=0;
00138       
00139       for (int k=0;k<weightFunctions.dimension(0);k++){
00140         denominator += weightFunctions(k,pointIndex);
00141       }
00142       for (int k=0;k<outputValues.dimension(0);k++){
00143         outputValues(k,pointIndex) = weightFunctions(k,pointIndex)/denominator;
00144       }
00145     }
00146   }
00147                                                             
00148 
00149 
00150   template<class Scalar, class ArrayScalar>
00151   template<class Scalar1, class ArrayScalar1>
00152   Scalar1 Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::computeArea(const ArrayScalar1& p1,
00153                                                                     const ArrayScalar1& p2,
00154                                                                     const ArrayScalar1& p3) const{
00155     Scalar1 area = 0.5*(p1(0)*(p2(1)-p3(1))
00156                         -p2(0)*(p1(1)-p3(1))
00157                         +p3(0)*(p1(1)-p2(1)));
00158     return area;
00159   }
00160   
00161 
00162   template<class Scalar, class ArrayScalar>
00163   template<class Scalar1, class ArrayScalar1>
00164   void Basis_HGRAD_POLY_C1_FEM<Scalar, ArrayScalar>::evaluateWeightFunctions(ArrayScalar1& outputValues,
00165                                                                              const ArrayScalar1& inputValues,
00166                                                                              const ArrayScalar1& cellVertices)const{
00167     
00168     
00169     int spaceDim = this->basisCellTopology_.getDimension();
00170     for (int k=0;k<outputValues.dimension(0);k++){
00171       
00172       // compute a_k for each weight function
00173       // need a_k = A(p_i,p_j,p_k) where nodes i and j are adjacent to node k
00174       int adjIndex1 = -1, adjIndex2 = -1;
00175       for (int i = 0;i < this->basisCellTopology_.getEdgeCount();i++){
00176         if ( this->basisCellTopology_.getNodeMap(1,i,0) == k )
00177           adjIndex1 = this->basisCellTopology_.getNodeMap(1,i,1);
00178         else if ( this->basisCellTopology_.getNodeMap(1,i,1) == k )
00179           adjIndex2 = this->basisCellTopology_.getNodeMap(1,i,0);
00180       }
00181       TEUCHOS_TEST_FOR_EXCEPTION( (adjIndex1 == -1 || adjIndex2 == -1), std::invalid_argument,
00182                           ">>> ERROR (Intrepid_HGRAD_POLY_C1_FEM): cannot find adjacent nodes when evaluating Wachspress weight function.");
00183       FieldContainer<Scalar1> p1(spaceDim);
00184       FieldContainer<Scalar1> p2(spaceDim);
00185       FieldContainer<Scalar1> p3(spaceDim);
00186       for (int i=0;i<spaceDim;i++){
00187         p1(i) = cellVertices(adjIndex1,i);
00188         p2(i) = cellVertices(k,i);
00189         p3(i) = cellVertices(adjIndex2,i);
00190       }
00191       Scalar1 a_k = computeArea<Scalar1, ArrayScalar1>(p1,p2,p3);
00192       // now compute prod_{ij!=k} r_ij
00193       for (int pointIndex = 0;pointIndex < inputValues.dimension(0);pointIndex++){
00194         Scalar1 product = a_k;
00195         for (int edgeIndex = 0;edgeIndex < this->basisCellTopology_.getEdgeCount();edgeIndex++){
00196           int edgeNode1 = this->basisCellTopology_.getNodeMap(1,edgeIndex,0);
00197           int edgeNode2 = this->basisCellTopology_.getNodeMap(1,edgeIndex,1);
00198           if ( edgeNode1 != k && edgeNode2 != k ){
00199             for (int i=0;i<spaceDim;i++){
00200               p1(i) = inputValues(pointIndex,i);
00201               p2(i) = cellVertices(edgeNode1,i);
00202               p3(i) = cellVertices(edgeNode2,i);
00203             }
00204             product *= computeArea<Scalar1, ArrayScalar1>(p1,p2,p3);
00205           }
00206         }
00207         outputValues(k,pointIndex) = product;
00208       }
00209     }
00210   }
00211 } // namespace Intrepid
00212            
00213