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