|
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) or 00026 // Robert Kirby (robert.c.kirby@ttu.edu) 00027 // 00028 // ************************************************************************ 00029 // @HEADER 00030 00037 namespace Intrepid { 00038 00044 static int idx(int p, int q); 00045 00064 template<typename Scalar> 00065 static void jrc(const Scalar &alpha , const Scalar &beta , 00066 const int &n , 00067 Scalar &an , Scalar &bn, Scalar &cn ); 00068 00069 00070 template<class Scalar, class ArrayScalar> 00071 Basis_HGRAD_TRI_Cn_FEM_ORTH<Scalar,ArrayScalar>::Basis_HGRAD_TRI_Cn_FEM_ORTH( int degree ) 00072 { 00073 this -> basisCardinality_ = (degree+1)*(degree+2)/2; 00074 this -> basisDegree_ = degree; 00075 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() ); 00076 this -> basisType_ = BASIS_FEM_HIERARCHICAL; 00077 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00078 this -> basisTagsAreSet_ = false; 00079 } 00080 00081 00082 00083 template<class Scalar, class ArrayScalar> 00084 void Basis_HGRAD_TRI_Cn_FEM_ORTH<Scalar, ArrayScalar>::initializeTags() { 00085 00086 // Basis-dependent initializations 00087 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00088 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00089 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00090 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00091 00092 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00093 int *tags = new int[tagSize * this->getCardinality()]; 00094 for (int i=0;i<this->getCardinality();i++) { 00095 tags[4*i] = 2; 00096 tags[4*i+1] = 0; 00097 tags[4*i+2] = i; 00098 tags[4*i+3] = this->getCardinality(); 00099 } 00100 00101 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00102 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00103 this -> ordinalToTag_, 00104 tags, 00105 this -> basisCardinality_, 00106 tagSize, 00107 posScDim, 00108 posScOrd, 00109 posDfOrd); 00110 } 00111 00112 00113 00114 template<class Scalar, class ArrayScalar> 00115 void Basis_HGRAD_TRI_Cn_FEM_ORTH<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00116 const ArrayScalar & inputPoints, 00117 const EOperator operatorType) const { 00118 00119 // Verify arguments 00120 #ifdef HAVE_INTREPID_DEBUG 00121 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00122 inputPoints, 00123 operatorType, 00124 this -> getBaseCellTopology(), 00125 this -> getCardinality() ); 00126 #endif 00127 const int deg = this->getDegree(); 00128 00129 // add more here and put in appropriate extra case statements below to enable higher derivatives. 00130 void (*tabulators[])(ArrayScalar &, const int, const ArrayScalar &) 00131 = { TabulatorTri<Scalar,ArrayScalar,0>::tabulate , 00132 TabulatorTri<Scalar,ArrayScalar,1>::tabulate , 00133 TabulatorTri<Scalar,ArrayScalar,2>::tabulate }; 00134 00135 00136 switch (operatorType) { 00137 case OPERATOR_VALUE: 00138 tabulators[0]( outputValues , deg , inputPoints ); 00139 break; 00140 case OPERATOR_GRAD: 00141 tabulators[1]( outputValues , deg , inputPoints ); 00142 break; 00143 case OPERATOR_D1: 00144 case OPERATOR_D2: 00145 // add more case OPEATOR_Dn statements if you've added more items to the 00146 // array above. 00147 tabulators[operatorType-OPERATOR_D1+1]( outputValues , deg , inputPoints ); 00148 break; 00149 default: 00150 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00151 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid or unsupported operator" ); 00152 00153 } 00154 00155 return; 00156 } 00157 00158 00159 00160 template<class Scalar, class ArrayScalar> 00161 void Basis_HGRAD_TRI_Cn_FEM_ORTH<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00162 const ArrayScalar & inputPoints, 00163 const ArrayScalar & cellVertices, 00164 const EOperator operatorType) const { 00165 TEST_FOR_EXCEPTION( (true), std::logic_error, 00166 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function"); 00167 } 00168 00169 00170 00171 template<typename Scalar, typename ArrayScalar> 00172 void TabulatorTri<Scalar,ArrayScalar,0>::tabulate(ArrayScalar &outputValues , 00173 const int deg , 00174 const ArrayScalar &z ) 00175 { 00176 const int np = z.dimension( 0 ); 00177 00178 // each point needs to be transformed from Pavel's element 00179 // z(i,0) --> (2.0 * z(i,0) - 1.0) 00180 // z(i,1) --> (2.0 * z(i,1) - 1.0) 00181 00182 // set up constant term 00183 int idx_cur = idx(0,0); 00184 int idx_curp1,idx_curm1; 00185 00186 // set D^{0,0} = 1.0 00187 for (int i=0;i<np;i++) { 00188 outputValues(idx_cur,i) = Scalar( 1.0 ) + z(i,0) - z(i,0) + z(i,1) - z(i,1); 00189 } 00190 00191 00192 if (deg > 0) { 00193 Teuchos::Array<Scalar> f1(np),f2(np),f3(np); 00194 00195 for (int i=0;i<np;i++) { 00196 f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0)); 00197 f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0)); 00198 f3[i] = f2[i] * f2[i]; 00199 } 00200 00201 // set D^{1,0} = f1 00202 idx_cur =idx(1,0); 00203 for (int i=0;i<np;i++) { 00204 outputValues(idx_cur,i) = f1[i]; 00205 } 00206 00207 // recurrence in p 00208 for (int p=1;p<deg;p++) { 00209 idx_cur = idx(p,0); 00210 idx_curp1 = idx(p+1,0); 00211 idx_curm1 = idx(p-1,0); 00212 Scalar a = (2.0*p+1.0)/(1.0+p); 00213 Scalar b = p / (p+1.0); 00214 00215 for (int i=0;i<np;i++) { 00216 outputValues(idx_curp1,i) = a * f1[i] * outputValues(idx_cur,i) 00217 - b * f3[i] * outputValues(idx_curm1,i); 00218 } 00219 } 00220 00221 // D^{p,1} 00222 for (int p=0;p<deg;p++) { 00223 int idxp0 = idx(p,0); 00224 int idxp1 = idx(p,1); 00225 for (int i=0;i<np;i++) { 00226 outputValues(idxp1,i) = outputValues(idxp0,i) 00227 *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0)); 00228 } 00229 } 00230 00231 00232 // recurrence in q 00233 for (int p=0;p<deg-1;p++) { 00234 for (int q=1;q<deg-p;q++) { 00235 int idxpqp1=idx(p,q+1); 00236 int idxpq=idx(p,q); 00237 int idxpqm1=idx(p,q-1); 00238 Scalar a,b,c; 00239 jrc<Scalar>((Scalar)(2*p+1),(Scalar)0,q,a,b,c); 00240 for (int i=0;i<np;i++) { 00241 outputValues(idxpqp1,i) 00242 = (a*(2.0*z(i,1)-1.0)+b)*outputValues(idxpq,i) 00243 - c*outputValues(idxpqm1,i); 00244 } 00245 } 00246 } 00247 } 00248 00249 // orthogonalize 00250 for (int p=0;p<=deg;p++) { 00251 for (int q=0;q<=deg-p;q++) { 00252 for (int i=0;i<np;i++) { 00253 outputValues(idx(p,q),i) *= sqrt( (p+0.5)*(p+q+1.0)); 00254 } 00255 } 00256 } 00257 00258 return; 00259 } 00260 00261 00262 00263 template<typename Scalar, typename ArrayScalar> 00264 void TabulatorTri<Scalar,ArrayScalar,1>::tabulate(ArrayScalar &outputValues , 00265 const int deg , 00266 const ArrayScalar &z ) 00267 { 00268 const int np = z.dimension(0); 00269 const int card = outputValues.dimension(0); 00270 FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) ); 00271 for (int i=0;i<np;i++) { 00272 for (int j=0;j<2;j++) { 00273 dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) ); 00274 dZ(i,j).diff(j,2); 00275 } 00276 } 00277 FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np); 00278 00279 TabulatorTri<Sacado::Fad::DFad<Scalar>,FieldContainer<Sacado::Fad::DFad<Scalar> >,0>::tabulate( dResult , 00280 deg , 00281 dZ ); 00282 00283 for (int i=0;i<card;i++) { 00284 for (int j=0;j<np;j++) { 00285 for (int k=0;k<2;k++) { 00286 outputValues(i,j,k) = dResult(i,j).dx(k); 00287 } 00288 } 00289 } 00290 00291 return; 00292 00293 } 00294 00295 00296 00297 template<typename Scalar, typename ArrayScalar, unsigned derivOrder> 00298 void TabulatorTri<Scalar,ArrayScalar,derivOrder>::tabulate( ArrayScalar &outputValues , 00299 const int deg , 00300 const ArrayScalar &z ) 00301 { 00302 const int np = z.dimension(0); 00303 const int card = outputValues.dimension(0); 00304 FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) ); 00305 for (int i=0;i<np;i++) { 00306 for (int j=0;j<2;j++) { 00307 dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) ); 00308 dZ(i,j).diff(j,2); 00309 } 00310 } 00311 FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np,derivOrder+1); 00312 00313 TabulatorTri<Sacado::Fad::DFad<Scalar>,FieldContainer<Sacado::Fad::DFad<Scalar> >,derivOrder-1>::tabulate(dResult , 00314 deg , 00315 dZ ); 00316 00317 for (int i=0;i<card;i++) { 00318 for (int j=0;j<np;j++) { 00319 outputValues(i,j,0) = dResult(i,j,0).dx(0); 00320 for (unsigned k=0;k<derivOrder;k++) { 00321 outputValues(i,j,k+1) = dResult(i,j,k).dx(1); 00322 } 00323 } 00324 } 00325 00326 return; 00327 00328 00329 } 00330 00331 00332 00333 int idx(int p , int q) 00334 { 00335 return (p+q)*(p+q+1)/2+q; 00336 } 00337 00338 00339 00340 template<typename Scalar> 00341 void jrc( const Scalar &alpha , const Scalar &beta , 00342 const int &n , 00343 Scalar &an , Scalar &bn, Scalar &cn ) 00344 { 00345 an = (2.0 * n + 1.0 + alpha + beta) * ( 2.0 * n + 2.0 + alpha + beta ) 00346 / ( 2.0 * ( n + 1 ) * ( n + 1 + alpha + beta ) ); 00347 bn = (alpha*alpha-beta*beta)*(2.0*n+1.0+alpha+beta) 00348 / ( 2.0*(n+1.0)*(2.0*n+alpha+beta)*(n+1.0+alpha+beta) ); 00349 cn = (n+alpha)*(n+beta)*(2.0*n+2.0+alpha+beta) 00350 / ( (n+1.0)*(n+1.0+alpha+beta)*(2.0*n+alpha+beta) ); 00351 00352 return; 00353 } 00354 00355 00356 }// namespace Intrepid
1.7.4