|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00035 namespace Intrepid { 00036 00037 template<class Scalar, class ArrayScalar> 00038 Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_LINE_Cn_FEM( const int n , 00039 const ArrayScalar &pts ): 00040 latticePts_( n+1 , 1 ), 00041 Phis_( n ), 00042 V_(n+1,n+1), 00043 Vinv_(n+1,n+1) 00044 { 00045 const int N = n+1; 00046 this -> basisCardinality_ = N; 00047 this -> basisDegree_ = n; 00048 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() ); 00049 this -> basisType_ = BASIS_FEM_FIAT; 00050 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00051 this -> basisTagsAreSet_ = false; 00052 00053 00054 // check validity of points 00055 for (int i=0;i<n;i++) { 00056 TEST_FOR_EXCEPTION( pts(i,0) >= pts(i+1,0) , 00057 std::runtime_error , 00058 "Intrepid::Basis_HGRAD_LINE_Cn_FEM Illegal points given to constructor" ); 00059 } 00060 00061 // copy points int latticePts, correcting endpoints if needed 00062 if (std::abs(pts(0,0)+1.0) < INTREPID_TOL) { 00063 latticePts_(0,0) = -1.0; 00064 } 00065 else { 00066 latticePts_(0,0) = pts(0,0); 00067 } 00068 for (int i=1;i<n;i++) { 00069 latticePts_(i,0) = pts(i,0); 00070 } 00071 if (std::abs(pts(n,0)-1.0) < INTREPID_TOL) { 00072 latticePts_(n,0) = 1.0; 00073 } 00074 else { 00075 latticePts_(n,0) = pts(n,0); 00076 } 00077 00078 // form Vandermonde matrix. Actually, this is the transpose of the VDM, 00079 // so we transpose on copy below. 00080 00081 Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE ); 00082 00083 // now I need to copy V into a Teuchos array to do the inversion 00084 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N); 00085 for (int i=0;i<N;i++) { 00086 for (int j=0;j<N;j++) { 00087 Vsdm(i,j) = V_(i,j); 00088 } 00089 } 00090 00091 // invert the matrix 00092 Teuchos::SerialDenseSolver<int,Scalar> solver; 00093 solver.setMatrix( rcp( &Vsdm , false ) ); 00094 solver.invert( ); 00095 00096 // now I need to copy the inverse into Vinv 00097 for (int i=0;i<N;i++) { 00098 for (int j=0;j<N;j++) { 00099 Vinv_(i,j) = Vsdm(j,i); 00100 } 00101 } 00102 00103 } 00104 00105 template<class Scalar, class ArrayScalar> 00106 Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_LINE_Cn_FEM( const int n , 00107 const EPointType &pointType ): 00108 latticePts_( n+1 , 1 ), 00109 Phis_( n ), 00110 V_(n+1,n+1), 00111 Vinv_(n+1,n+1) 00112 { 00113 const int N = n+1; 00114 this -> basisCardinality_ = N; 00115 this -> basisDegree_ = n; 00116 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() ); 00117 this -> basisType_ = BASIS_FEM_FIAT; 00118 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00119 this -> basisTagsAreSet_ = false; 00120 00121 switch(pointType) { 00122 case POINTTYPE_EQUISPACED: 00123 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ , this->basisCellTopology_ , n , 0 , POINTTYPE_EQUISPACED ); 00124 break; 00125 case POINTTYPE_SPECTRAL: 00126 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ , this->basisCellTopology_ , n , 0 , POINTTYPE_WARPBLEND ); 00127 break; 00128 case POINTTYPE_SPECTRAL_OPEN: 00129 PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( latticePts_ , n ); 00130 break; 00131 default: 00132 TEST_FOR_EXCEPTION( true , std::invalid_argument , "Basis_HGRAD_LINE_Cn_FEM:: invalid point type" ); 00133 break; 00134 } 00135 00136 // form Vandermonde matrix. Actually, this is the transpose of the VDM, 00137 // so we transpose on copy below. 00138 00139 Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE ); 00140 00141 // now I need to copy V into a Teuchos array to do the inversion 00142 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N); 00143 for (int i=0;i<N;i++) { 00144 for (int j=0;j<N;j++) { 00145 Vsdm(i,j) = V_(i,j); 00146 } 00147 } 00148 00149 // invert the matrix 00150 Teuchos::SerialDenseSolver<int,Scalar> solver; 00151 solver.setMatrix( rcp( &Vsdm , false ) ); 00152 solver.invert( ); 00153 00154 // now I need to copy the inverse into Vinv 00155 for (int i=0;i<N;i++) { 00156 for (int j=0;j<N;j++) { 00157 Vinv_(i,j) = Vsdm(j,i); 00158 } 00159 } 00160 } 00161 00162 00163 template<class Scalar, class ArrayScalar> 00164 void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::initializeTags() { 00165 00166 // Basis-dependent initializations 00167 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00168 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00169 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00170 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00171 00172 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00173 00174 int *tags = new int[ tagSize * this->getCardinality() ]; 00175 00176 int internal_dof; 00177 int edge_dof; 00178 00179 const int n = this->getDegree(); 00180 00181 // now we check the points for association 00182 if (latticePts_(0,0) == -1.0) { 00183 tags[0] = 0; 00184 tags[1] = 0; 00185 tags[2] = 0; 00186 tags[3] = 1; 00187 edge_dof = 1; 00188 internal_dof = n-1; 00189 } 00190 else { 00191 tags[0] = 1; 00192 tags[1] = 0; 00193 tags[2] = 0; 00194 tags[3] = n+1; 00195 edge_dof = 0; 00196 internal_dof = n+1; 00197 } 00198 for (int i=1;i<n;i++) { 00199 tags[4*i] = 1; 00200 tags[4*i+1] = 0; 00201 tags[4*i+2] = -edge_dof + i; 00202 tags[4*i+3] = internal_dof; 00203 } 00204 if (latticePts_(n,0) == 1.0) { 00205 tags[4*n] = 0; 00206 tags[4*n+1] = 1; 00207 tags[4*n+2] = 0; 00208 tags[4*n+3] = 1; 00209 } 00210 else { 00211 tags[4*n] = 1; 00212 tags[4*n+1] = 0; 00213 tags[4*n+2] = n; 00214 tags[4*n+3] = n; 00215 } 00216 00217 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00218 this -> ordinalToTag_, 00219 tags, 00220 this -> basisCardinality_, 00221 tagSize, 00222 posScDim, 00223 posScOrd, 00224 posDfOrd); 00225 00226 delete []tags; 00227 00228 } 00229 00230 00231 00232 template<class Scalar, class ArrayScalar> 00233 void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00234 const ArrayScalar & inputPoints, 00235 const EOperator operatorType) const { 00236 00237 // Verify arguments 00238 #ifdef HAVE_INTREPID_DEBUG 00239 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00240 inputPoints, 00241 operatorType, 00242 this -> getBaseCellTopology(), 00243 this -> getCardinality() ); 00244 #endif 00245 const int numPts = inputPoints.dimension(0); 00246 const int numBf = this->getCardinality(); 00247 00248 try { 00249 switch (operatorType) { 00250 case OPERATOR_VALUE: 00251 { 00252 FieldContainer<Scalar> phisCur( numBf , numPts ); 00253 Phis_.getValues( phisCur , inputPoints , operatorType ); 00254 for (int i=0;i<outputValues.dimension(0);i++) { 00255 for (int j=0;j<outputValues.dimension(1);j++) { 00256 outputValues(i,j) = 0.0; 00257 for (int k=0;k<this->getCardinality();k++) { 00258 outputValues(i,j) += this->Vinv_(k,i) * phisCur(k,j); 00259 } 00260 } 00261 } 00262 } 00263 break; 00264 case OPERATOR_GRAD: 00265 case OPERATOR_D1: 00266 case OPERATOR_D2: 00267 case OPERATOR_D3: 00268 case OPERATOR_D4: 00269 case OPERATOR_D5: 00270 case OPERATOR_D6: 00271 case OPERATOR_D7: 00272 case OPERATOR_D8: 00273 case OPERATOR_D9: 00274 case OPERATOR_D10: 00275 { 00276 const int dkcard = 00277 (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,1): getDkCardinality(operatorType,1); 00278 00279 FieldContainer<Scalar> phisCur( numBf , numPts , dkcard ); 00280 Phis_.getValues( phisCur , inputPoints , operatorType ); 00281 00282 for (int i=0;i<outputValues.dimension(0);i++) { 00283 for (int j=0;j<outputValues.dimension(1);j++) { 00284 for (int k=0;k<outputValues.dimension(2);k++) { 00285 outputValues(i,j,k) = 0.0; 00286 for (int l=0;l<this->getCardinality();l++) { 00287 outputValues(i,j,k) += this->Vinv_(l,i) * phisCur(l,j,k); 00288 } 00289 } 00290 } 00291 } 00292 } 00293 break; 00294 default: 00295 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00296 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator type not implemented" ); 00297 break; 00298 } 00299 } 00300 catch (std::invalid_argument &exception){ 00301 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00302 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator failed"); 00303 } 00304 00305 } 00306 00307 00308 00309 template<class Scalar, class ArrayScalar> 00310 void Basis_HGRAD_LINE_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00311 const ArrayScalar & inputPoints, 00312 const ArrayScalar & cellVertices, 00313 const EOperator operatorType) const { 00314 TEST_FOR_EXCEPTION( (true), std::logic_error, 00315 ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): FEM Basis calling an FVD member function"); 00316 } 00317 00318 00319 }// namespace Intrepid
1.7.4