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