Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Discretization/TensorProductSpaceTools/Intrepid_TensorProductSpaceToolsDef.hpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00050 namespace Intrepid {
00051   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00052            class ArrayTypeBasis>
00053   void TensorProductSpaceTools::evaluate( ArrayTypeOut &vals ,
00054                                           const ArrayTypeCoeffs &coeffs ,
00055                                           const Array<RCP<ArrayTypeBasis> > &bases )
00056   {
00057     const unsigned space_dim = bases.size();
00058 
00059 #ifdef HAVE_INTREPID_DEBUG
00060     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00061                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00062 
00063     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00064                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00065 
00066     for (unsigned d=0;d<space_dim;d++)
00067       {
00068         TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
00069                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00070       }
00071 
00072     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00073                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00074 
00075     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00076                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00077     
00078     int product_of_basis_dimensions = 1;
00079     int product_of_basis_points = 1;
00080     for (unsigned d=0;d<space_dim;d++)
00081       {
00082         product_of_basis_dimensions *= bases[d]->dimension(0);
00083         product_of_basis_points *= bases[d]->dimension(1);
00084       }
00085     
00086     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00087                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00088     
00089     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00090                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
00091 #endif    
00092     
00093 
00094     switch (space_dim)
00095       {
00096       case 2:
00097         evaluate2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00098         break;
00099       case 3:
00100         evaluate3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00101         break;
00102       }
00103 
00104   }
00105 
00106   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00107            class ArrayTypeBasis>
00108   void TensorProductSpaceTools::evaluateGradient( ArrayTypeOut &vals ,
00109                                                   const ArrayTypeCoeffs &coeffs ,
00110                                                   const Array<RCP<ArrayTypeBasis> > &bases ,
00111                                                   const Array<RCP<ArrayTypeBasis> > &Dbases )
00112   {
00113     const unsigned space_dim = bases.size();
00114 
00115 #ifdef HAVE_INTREPID_DEBUG
00116     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
00117                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00118 
00119     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00120                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00121 
00122     TEUCHOS_TEST_FOR_EXCEPTION( Dbases.size() != (int)space_dim , std::invalid_argument ,
00123                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases and Dbases must be same size." );
00124    
00125     for (unsigned d=0;d<space_dim;d++)
00126       {
00127         TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
00128                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00129 
00130         TEUCHOS_TEST_FOR_EXCEPTION( Dbases[d]->rank() != 3 , std::invalid_argument ,
00131                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 3 array." );
00132       }
00133 
00134     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00135                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00136 
00137     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00138                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00139     
00140     int product_of_basis_dimensions = 1;
00141     int product_of_basis_points = 1;
00142     for (unsigned d=0;d<space_dim;d++)
00143       {
00144         product_of_basis_dimensions *= bases[d]->dimension(0);
00145         product_of_basis_points *= bases[d]->dimension(1);
00146       }
00147     
00148     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00149                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00150     
00151     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00152                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
00153 
00154     for (unsigned d=0;d<space_dim;d++)
00155       {
00156         for (unsigned i=0;i<2;i++)
00157           {
00158             TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->dimension(i) != Dbases[d]->dimension(i) , std::invalid_argument ,
00159                                         ">>>ERROR (TensorProductSpaceTools::evaluate): bases and Dbases have incompatible shape." );
00160           }
00161       }
00162 #endif    
00163 
00164     switch (space_dim)
00165       {
00166       case 2:
00167         evaluateGradient2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00168         break;
00169       case 3:
00170         evaluateGradient3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00171         break;
00172       }
00173 
00174   }
00175 
00176   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00177            class ArrayTypeBasis, class ArrayTypeWeights>
00178   void TensorProductSpaceTools::moments( ArrayTypeOut &vals ,
00179                                          const ArrayTypeData &data ,
00180                                          const Array<RCP<ArrayTypeBasis> > &basisVals ,
00181                                          const Array<RCP<ArrayTypeWeights> > &wts )
00182   {
00183     const unsigned space_dim = basisVals.size();
00184 
00185 #ifdef HAVE_INTREPID_DEBUG
00186     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00187                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00188 
00189     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 3 , std::invalid_argument ,
00190                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00191 
00192     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00193                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00194 
00195     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00196                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00197    
00198     for (unsigned d=0;d<space_dim;d++)
00199       {
00200         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00201                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00202 
00203         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00204                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00205       }
00206 
00207     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00208                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00209 
00210     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != data.dimension(1) , std::invalid_argument,
00211                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00212     
00213     int product_of_basis_dimensions = 1;
00214     int product_of_basis_points = 1;
00215     for (unsigned d=0;d<space_dim;d++)
00216       {
00217         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00218         product_of_basis_points *= basisVals[d]->dimension(1);
00219       }
00220     
00221     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00222                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00223     
00224     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
00225                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00226 
00227 #endif    
00228 
00229     switch (space_dim)
00230       {
00231       case 2:
00232         moments2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00233         break;
00234       case 3:
00235         moments3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00236         break;
00237       }
00238 
00239   }
00240 
00241   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00242            class ArrayTypeBasis, class ArrayTypeWeights>
00243   void TensorProductSpaceTools::momentsGrad( ArrayTypeOut &vals ,
00244                                              const ArrayTypeData &data ,
00245                                              const Array<RCP<ArrayTypeBasis> > &basisVals ,
00246                                              const Array<RCP<ArrayTypeBasis> > &basisDVals ,
00247                                              const Array<RCP<ArrayTypeWeights> > &wts )
00248   {
00249     const unsigned space_dim = basisVals.size();
00250 
00251 #ifdef HAVE_INTREPID_DEBUG
00252     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00253                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00254 
00255     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 4 , std::invalid_argument ,
00256                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 4 array." );
00257 
00258     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00259                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00260 
00261     TEUCHOS_TEST_FOR_EXCEPTION( basisDVals.size() != (int)space_dim , std::invalid_argument ,
00262                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00263 
00264     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00265                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00266    
00267     for (unsigned d=0;d<space_dim;d++)
00268       {
00269         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00270                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00271 
00272         TEUCHOS_TEST_FOR_EXCEPTION( basisDVals[d]->rank() != 3 , std::invalid_argument ,
00273                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated derivative basis must be rank 3 array." );
00274 
00275         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00276                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00277       }
00278 
00279     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00280                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00281     
00282     int product_of_basis_dimensions = 1;
00283     int product_of_basis_points = 1;
00284     for (unsigned d=0;d<space_dim;d++)
00285       {
00286         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00287         product_of_basis_points *= basisVals[d]->dimension(1);
00288       }
00289     
00290     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != product_of_basis_dimensions , std::invalid_argument,
00291                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00292     
00293     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(1) != product_of_basis_points , std::invalid_argument,
00294                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00295 
00296 #endif    
00297 
00298     switch (space_dim)
00299       {
00300       case 2:
00301         momentsGrad2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00302         break;
00303       case 3:
00304         momentsGrad3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00305         break;
00306       }
00307 
00308   }
00309 
00310 
00311   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00312            class ArrayTypeBasis>
00313   void TensorProductSpaceTools::evaluate2D( ArrayTypeOut &vals ,
00314                                             const ArrayTypeCoeffs &coeffs ,
00315                                             const Array<RCP<ArrayTypeBasis> > &bases )
00316   {
00317     const int numBfx = bases[0]->dimension(0);
00318     const int numBfy = bases[1]->dimension(0);
00319     
00320     const int numPtsx = bases[0]->dimension(1);
00321     const int numPtsy = bases[1]->dimension(1);
00322     
00323     const int numCells = vals.dimension(0);
00324     const int numFields = vals.dimension(1);
00325 
00326     const ArrayTypeBasis &Phix = *bases[0];
00327     const ArrayTypeBasis &Phiy = *bases[1];
00328     
00329     FieldContainer<double> Xi(numCells,numBfy,numPtsx);
00330     
00331     // sum factorization step
00332 
00333     for (int f=0;f<numFields;f++)
00334       {
00335         for (int cell=0;cell<numCells;cell++)
00336           {
00337             for (int j=0;j<numBfy;j++)
00338               {
00339                 for (int i=0;i<numBfx;i++)
00340                   {
00341                     const int I = j * numBfx + i;
00342                     for (int k=0;k<numPtsx;k++)
00343                       {
00344                         Xi(cell,j,k) += coeffs(cell,f,I) * Phix(i,k);
00345                       }
00346                   }
00347               }
00348           }
00349         
00350         for (int cell=0;cell<numCells;cell++)
00351           {
00352             for (int kpty=0;kpty<numPtsy;kpty++)
00353               {
00354                 for (int kptx=0;kptx<numPtsx;kptx++)
00355                   {
00356                     vals(cell,f,kptx+numPtsx*kpty) = 0.0;
00357                   }
00358               }
00359             
00360             // evaluation step
00361             for (int kpty=0;kpty<numPtsy;kpty++)
00362               {
00363                 for (int kptx=0;kptx<numPtsx;kptx++)
00364                   {
00365                     const int I=kptx+numPtsx*kpty;
00366                     
00367                     for (int j=0;j<numBfy;j++)
00368                       {
00369                         vals(cell,f,I) += Phiy(j,kpty) * Xi(cell,j,kptx);
00370                       }
00371                   }
00372               }
00373           }
00374       }
00375     
00376     return;
00377   }
00378 
00379   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00380            class ArrayTypeBasis>
00381   void TensorProductSpaceTools::evaluate3D( ArrayTypeOut &vals ,
00382                                             const ArrayTypeCoeffs &coeffs ,
00383                                             const Array<RCP<ArrayTypeBasis> > &bases )  
00384 
00385   {
00386     // checks to do:
00387     // vals point dimension is product of sizes of point arrays
00388     // vals
00389     
00390     const int numBfx = bases[0]->dimension(0);
00391     const int numBfy = bases[1]->dimension(0);
00392     const int numBfz = bases[2]->dimension(0);
00393     
00394     const int numPtsx = bases[0]->dimension(1);
00395     const int numPtsy = bases[1]->dimension(1);
00396     const int numPtsz = bases[2]->dimension(1);
00397     
00398     const int numCells = vals.dimension(0);
00399     const int numFields = vals.dimension(1);
00400     
00401     const ArrayTypeBasis &Phix = *bases[0];
00402     const ArrayTypeBasis &Phiy = *bases[1];
00403     const FieldContainer<double> &Phiz = *bases[2];
00404 
00405     FieldContainer<double> Xi(numCells, 
00406                               numBfz, numBfy , numPtsx);
00407     
00408     FieldContainer<double> Theta(numCells, 
00409                                  numBfz , numPtsy, numPtsx );
00410     
00411     
00412     for (int f=0;f<numFields;f++)
00413       {
00414 
00415         // Xi section
00416         for (int c=0;c<numCells;c++)
00417           {
00418             for (int k=0;k<numBfz;k++)
00419               {
00420                 for (int j=0;j<numBfy;j++)
00421                   {
00422                     for (int l=0;l<numPtsx;l++)
00423                       {
00424                         for (int i=0;i<numBfx;i++)
00425                           {
00426                             const int coeffIndex = k*numBfz*numBfx + j * numBfx + i;
00427                             Xi(c,k,j,l) += Phix(i,l) * coeffs(c,f,coeffIndex);
00428                           }
00429                       }
00430                   }
00431               }
00432           }
00433     
00434         // Theta section
00435         for (int c=0;c<numCells;c++)
00436           {
00437             for (int k=0;k<numBfz;k++)
00438               {
00439                 for (int m=0;m<numPtsy;m++)
00440                   {
00441                     for (int l=0;l<numPtsx;l++)
00442                       {
00443                         for (int j=0;j<numBfy;j++)
00444                           {
00445                             Theta(c,k,m,l) += Phiy(j,m) * Xi(c,k,j,l);
00446                           }
00447                       }
00448                   }
00449               }
00450           }
00451     
00452         // final section
00453         for (int c=0;c<numCells;c++)
00454           {
00455             for (int n=0;n<numPtsz;n++)
00456               {
00457                 for (int m=0;m<numPtsy;m++)
00458                   {
00459                     for (int l=0;l<numPtsx;l++)
00460                       {
00461                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) = 0.0;
00462                         for (int k=0;k<numBfz;k++)
00463                           {
00464                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) += Theta(c,k,m,l) * Phiz(k,n);
00465                           }
00466                       }
00467                   }
00468               }
00469           }
00470       }
00471 
00472     return;
00473     
00474   }
00475 
00476 
00477   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00478            class ArrayTypeBasis>
00479   void TensorProductSpaceTools::evaluateGradient2D( ArrayTypeOut &vals ,
00480                                                     const ArrayTypeCoeffs &coeffs ,
00481                                                     const Array<RCP<ArrayTypeBasis> > &bases ,
00482                                                     const Array<RCP<ArrayTypeBasis> > &Dbases )
00483   {
00484     
00485     const int numBfx = bases[0]->dimension(0);
00486     const int numBfy = bases[1]->dimension(0);
00487     const int numPtsx = bases[0]->dimension(1);
00488     const int numPtsy = bases[1]->dimension(1);
00489     const int numCells = vals.dimension(0);
00490     const int numFields = vals.dimension(1);
00491     const ArrayTypeBasis &Phix = *bases[0];
00492     const ArrayTypeBasis &Phiy = *bases[1];
00493     const ArrayTypeBasis &DPhix = *Dbases[0];
00494     const ArrayTypeBasis &DPhiy = *Dbases[1];
00495     
00496     FieldContainer<double> Xi(numBfy,numPtsx);
00497 
00498     for (int f=0;f<numFields;f++)
00499       {
00500     
00501         for (int cell=0;cell<numCells;cell++)
00502           {
00503             // x derivative
00504         
00505             // sum factorization step
00506             for (int j=0;j<numBfy;j++)
00507               {
00508                 for (int k=0;k<numPtsx;k++)
00509                   {
00510                     Xi(j,k) = 0.0;
00511                   }
00512               }
00513         
00514             for (int j=0;j<numBfy;j++)
00515               {
00516                 for (int i=0;i<numBfx;i++)
00517                   {
00518                     const int I = j * numBfx + i;
00519                     for (int k=0;k<numPtsx;k++)
00520                       {
00521                         Xi(j,k) += coeffs(cell,f,I) * DPhix(i,k,0);
00522                       }
00523                   }
00524               }
00525         
00526             for (int kpty=0;kpty<numPtsy;kpty++)
00527               {
00528                 for (int kptx=0;kptx<numPtsx;kptx++)
00529                   {
00530                     vals(cell,f,kptx+numPtsx*kpty,0) = 0.0;
00531                   }
00532               }
00533         
00534             // evaluation step
00535             for (int kpty=0;kpty<numPtsy;kpty++)
00536               {
00537                 for (int kptx=0;kptx<numPtsx;kptx++)
00538                   {
00539                     const int I=kptx+numPtsx*kpty;
00540                 
00541                     for (int j=0;j<numBfy;j++)
00542                       {
00543                         vals(cell,f,I,0) += Phiy(j,kpty) * Xi(j,kptx);
00544                       }
00545                   }
00546               }
00547         
00548             // y derivative
00549         
00550             // sum factorization step
00551             for (int j=0;j<numBfy;j++)
00552               {
00553                 for (int k=0;k<numPtsx;k++)
00554                   {
00555                     Xi(j,k) = 0.0;
00556                   }
00557               }
00558         
00559             for (int j=0;j<numBfy;j++)
00560               {
00561                 for (int i=0;i<numBfx;i++)
00562                   {
00563                     const int I = j * numBfx + i;
00564                     for (int k=0;k<numPtsx;k++)
00565                       {
00566                         Xi(j,k) += coeffs(cell,f,I) * Phix(i,k);
00567                       }
00568                   }
00569               }
00570         
00571             for (int kpty=0;kpty<numPtsy;kpty++)
00572               {
00573                 for (int kptx=0;kptx<numPtsx;kptx++)
00574                   {
00575                     vals(cell,f,kptx+numPtsx*kpty,1) = 0.0;
00576                   }
00577               }
00578         
00579             // evaluation step
00580             for (int kpty=0;kpty<numPtsy;kpty++)
00581               {
00582                 for (int kptx=0;kptx<numPtsx;kptx++)
00583                   {
00584                     const int I=kptx+numPtsx*kpty;
00585                 
00586                     for (int j=0;j<numBfy;j++)
00587                       {
00588                         vals(cell,f,I,1) += DPhiy(j,kpty,0) * Xi(j,kptx);
00589                       }
00590                   }
00591               }
00592           }
00593       }    
00594     return;
00595   }
00596 
00597   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00598            class ArrayTypeBasis>
00599   void TensorProductSpaceTools::evaluateGradient3D( ArrayTypeOut &vals ,
00600                                                     const ArrayTypeCoeffs &coeffs ,
00601                                                     const Array<RCP<ArrayTypeBasis> > &bases ,
00602                                                     const Array<RCP<ArrayTypeBasis> > &Dbases )
00603   
00604   {
00605     // checks to do:
00606     // vals point dimension is product of sizes of point arrays
00607     // vals
00608 
00609     const int numBfx = bases[0]->dimension(0);
00610     const int numBfy = bases[1]->dimension(0);
00611     const int numBfz = bases[2]->dimension(0);
00612     const int numPtsx = bases[0]->dimension(1);
00613     const int numPtsy = bases[1]->dimension(1);
00614     const int numPtsz = bases[2]->dimension(1);
00615     const int numCells = vals.dimension(0);
00616     const int numFields = vals.dimension(1);
00617     const ArrayTypeBasis &Phix = *bases[0];
00618     const ArrayTypeBasis &Phiy = *bases[1];
00619     const ArrayTypeBasis &Phiz = *bases[2];
00620     const ArrayTypeBasis &DPhix = *Dbases[0];
00621     const ArrayTypeBasis &DPhiy = *Dbases[1];
00622     const ArrayTypeBasis &DPhiz = *Dbases[2];
00623 
00624     FieldContainer<double> Xi(numCells, 
00625                               numBfz, numBfy , numPtsx , 3) ;
00626 
00627     FieldContainer<double> Theta(numCells, 
00628                                  numBfz , numPtsy, numPtsx , 3);
00629   
00630     for (int f=0;f<numFields;f++)
00631       {
00632 
00633         // Xi section
00634         for (int c=0;c<numCells;c++)
00635           {
00636             for (int k=0;k<numBfz;k++)
00637               {
00638                 for (int j=0;j<numBfy;j++)
00639                   {
00640                     for (int l=0;l<numPtsx;l++)
00641                       {
00642                         for (int i=0;i<numBfx;i++)
00643                           {
00644                             const int coeffIndex = k*numBfz*numBfx + j * numBfx + i;
00645                             Xi(c,k,j,l,0) += DPhix(i,l,0) * coeffs(c,f,coeffIndex);
00646                             Xi(c,k,j,l,1) += Phix(i,l) * coeffs(c,f,coeffIndex);
00647                             Xi(c,k,j,l,2) += Phix(i,l) * coeffs(c,f,coeffIndex);
00648                           }
00649                       }
00650                   }
00651               }
00652           }
00653 
00654         // Theta section
00655         for (int c=0;c<numCells;c++)
00656           {
00657             for (int k=0;k<numBfz;k++)
00658               {
00659                 for (int m=0;m<numPtsy;m++)
00660                   {
00661                     for (int l=0;l<numPtsx;l++)
00662                       {
00663                         for (int j=0;j<numBfy;j++)
00664                           {
00665                             Theta(c,k,m,l,0) += Phiy(j,m) * Xi(c,k,j,l,0);
00666                             Theta(c,k,m,l,1) += DPhiy(j,m,0) * Xi(c,k,j,l,1);
00667                             Theta(c,k,m,l,2) += Phiy(j,m) * Xi(c,k,j,l,2);
00668                           }
00669                       }
00670                   }
00671               }
00672           }
00673 
00674         // final section
00675         for (int c=0;c<numCells;c++)
00676           {
00677             for (int n=0;n<numPtsz;n++)
00678               {
00679                 for (int m=0;m<numPtsy;m++)
00680                   {
00681                     for (int l=0;l<numPtsx;l++)
00682                       {
00683                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) = 0.0;
00684                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) = 0.0;
00685                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) = 0.0;
00686                         for (int k=0;k<numBfz;k++)
00687                           {
00688                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) += Theta(c,k,m,l,0) * Phiz(k,n);
00689                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) += Theta(c,k,m,l,1) * Phiz(k,n);
00690                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) += Theta(c,k,m,l,2) * DPhiz(k,n,0);
00691 
00692                           }
00693                       }
00694                   }
00695               }
00696           }
00697       }
00698     return;
00699   }
00700 
00701   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00702            class ArrayTypeBasis, class ArrayTypeWeights>
00703   void TensorProductSpaceTools::moments2D( ArrayTypeOut &vals ,
00704                                            const ArrayTypeData &data ,
00705                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
00706                                            const Array<RCP<ArrayTypeWeights> > &wts )
00707   {
00708     const int numBfx = basisVals[0]->dimension(0);
00709     const int numBfy = basisVals[1]->dimension(0);
00710     const int numPtsx = basisVals[0]->dimension(1);
00711     const int numPtsy = basisVals[1]->dimension(1);
00712     const int numCells = vals.dimension(0);
00713     const int numFields = vals.dimension(1);
00714     const ArrayTypeBasis &Phix = *basisVals[0];
00715     const ArrayTypeBasis &Phiy = *basisVals[1];
00716     
00717     FieldContainer<double> Xi(numBfx,numPtsy);
00718     
00719     const ArrayTypeWeights &wtsx = *wts[0];
00720     const ArrayTypeWeights &wtsy = *wts[1];
00721 
00722     for (int f=0;f<numFields;f++)
00723       {
00724         for (int cell=0;cell<numCells;cell++)
00725           {
00726             // sum factorization step
00727             for (int i=0;i<numBfx;i++)
00728               {
00729                 for (int k=0;k<numPtsy;k++)
00730                   {
00731                     Xi(i,k) = 0.0;
00732                   }
00733               }
00734         
00735             for (int i=0;i<numBfx;i++)
00736               {
00737                 for (int l=0;l<numPtsy;l++)
00738                   {
00739                     for (int k=0;k<numPtsx;k++)
00740                       {
00741                         Xi(i,l) += wtsx(k) * data(cell,f,l*numPtsx+k) * Phix(i,k);
00742                       }
00743                   }
00744               }
00745         
00746             for (int j=0;j<numBfy;j++)
00747               {
00748                 for (int i=0;i<numBfx;i++)
00749                   {
00750                     vals(cell,f,j*numBfx+i) = 0.0;
00751                   }
00752               }
00753 
00754             // evaluate moments with sum factorization
00755             for (int j=0;j<numBfy;j++)
00756               {
00757                 for (int i=0;i<numBfx;i++)
00758                   {
00759                     for (int l=0;l<numPtsy;l++)
00760                       {
00761                         vals(cell,f,j*numBfx+i) += wtsy(l) * Phiy(j,l) * Xi(i,l);
00762                       }
00763                   }
00764               }
00765           }
00766       }
00767     return;
00768 
00769   }
00770 
00771   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00772            class ArrayTypeBasis, class ArrayTypeWeights>
00773   void TensorProductSpaceTools::moments3D( ArrayTypeOut &vals ,
00774                                            const ArrayTypeData &data ,
00775                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
00776                                            const Array<RCP<ArrayTypeWeights> > &wts )
00777   {
00778     const int numBfx = basisVals[0]->dimension(0);
00779     const int numBfy = basisVals[1]->dimension(0);
00780     const int numBfz = basisVals[2]->dimension(0);
00781 
00782     const int numPtsx = basisVals[0]->dimension(1);
00783     const int numPtsy = basisVals[1]->dimension(1);
00784     const int numPtsz = basisVals[2]->dimension(1);
00785 
00786     const int numCells = vals.dimension(0);
00787     
00788     const ArrayTypeBasis &Phix = *basisVals[0];
00789     const ArrayTypeBasis &Phiy = *basisVals[1];
00790     const ArrayTypeBasis &Phiz = *basisVals[2];
00791 
00792     const ArrayTypeWeights &Wtx = *wts[0];
00793     const ArrayTypeWeights &Wty = *wts[1];
00794     const ArrayTypeWeights &Wtz = *wts[2];
00795 
00796     FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy);
00797     FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz);
00798 
00799     for (int f=0;f<numCells;f++)
00800       {
00801         // Xi phase
00802         for (int c=0;c<numCells;c++)
00803           {
00804             for (int i=0;i<numBfx;i++)
00805               {
00806                 for (int n=0;n<numPtsz;n++)
00807                   {
00808                     for (int m=0;m<numPtsy;m++)
00809                       {
00810                         for (int l=0;l<numPtsx;l++)
00811                           {
00812                             Xi(c,i,n,m) += Wtx(l) * Phix(i,l) * data(c,f,n*numPtsy*numPtsx+m*numPtsx+l);
00813                           }
00814                       }
00815                   }
00816               }
00817           }
00818 
00819         // Theta phase
00820         for (int c=0;c<numCells;c++)
00821           {
00822             for (int j=0;j<numBfy;j++)
00823               {
00824                 for (int i=0;i<numBfx;i++)
00825                   {
00826                     for (int n=0;n<numPtsz;n++)
00827                       {
00828                         for (int m=0;m<numPtsy;m++)
00829                           {
00830                             Theta(c,j,i,n) += Wty(m) * Phiy(j,m) * Xi(c,i,n,m);
00831                           }
00832                       }
00833                   }
00834               }
00835           }
00836 
00837         // Final phase
00838         for (int c=0;c<numCells;c++)
00839           {
00840             for (int k=0;k<numBfz;k++)
00841               {
00842                 for (int j=0;j<numBfy;j++)
00843                   {
00844                     for (int i=0;i<numBfx;i++)
00845                       {
00846                         const int momIdx = k*numBfx*numBfy+j*numBfx+i;
00847                         for (int n=0;n<numPtsz;n++)
00848                           {
00849                             vals(c,f,momIdx) += Wtz(n) * Phiz(k,n) * Theta(c,j,i,n);
00850                           }
00851                       }
00852                   }
00853               }
00854           }
00855 
00856       }
00857     return;
00858 
00859   }
00860 
00861   // data is now (C,P,D)
00862   // want to compute the moments against the gradients of the basis 
00863   // functions.
00864 
00865   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00866            class ArrayTypeBasis, class ArrayTypeWeights>
00867   void TensorProductSpaceTools::momentsGrad2D( ArrayTypeOut &vals ,
00868                                                const ArrayTypeData &data ,
00869                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
00870                                                const Array<RCP<ArrayTypeBasis> > &Dbases ,
00871                                                const Array<RCP<ArrayTypeWeights> > &wts )
00872   {
00873     
00874     const int numBfx = basisVals[0]->dimension(0);
00875     const int numBfy = basisVals[1]->dimension(0);
00876     const int numPtsx = basisVals[0]->dimension(1);
00877     const int numPtsy = basisVals[1]->dimension(1);
00878     const int numCells = vals.dimension(0);
00879     const int numFields = vals.dimension(1);
00880     const ArrayTypeBasis &Phix = *basisVals[0];
00881     const ArrayTypeBasis &Phiy = *basisVals[1];
00882     const ArrayTypeBasis &DPhix = *Dbases[0];
00883     const ArrayTypeBasis &DPhiy = *Dbases[1];
00884     const ArrayTypeWeights &wtsx = *wts[0];
00885     const ArrayTypeWeights &wtsy = *wts[1];
00886     
00887     FieldContainer<double> Xi(numBfx,numPtsy,2);
00888 
00889     for (int f=0;f<numFields;f++)
00890       {
00891         // Xi phase
00892         for (int c=0;c<numCells;c++)
00893           {
00894             for (int i=0;i<numBfx;i++)
00895               {
00896                 for (int m=0;m<numPtsy;m++)
00897                   {
00898                     for (int l=0;l<numPtsx;l++)
00899                       {
00900                         Xi(i,m,0) += wtsx(l) * data(c,f,m*numPtsy+l,0) * DPhix(i,l);
00901                         Xi(i,m,1) += wtsx(l) * data(c,f,m*numPtsy+l,1) * Phix(i,l);
00902                       }
00903                   }
00904               }
00905           }
00906 
00907         // main phase
00908         for (int c=0;c<numCells;c++)
00909           {
00910             for (int j=0;j<numBfy;j++)
00911               {
00912                 for (int i=0;i<numBfx;i++)
00913                   {
00914                     const int bfIdx = j*numBfx+i;
00915                     vals(c,f,bfIdx) = 0.0;
00916                     for (int m=0;m<numPtsy;m++)
00917                       {
00918                         vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,0) * Phiy(j,m);
00919                         vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,1) * DPhiy(j,m);
00920                       }
00921                   }
00922               }
00923           }
00924       }
00925     return;
00926 
00927   }
00928 
00929   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00930            class ArrayTypeBasis, class ArrayTypeWeights>
00931   void TensorProductSpaceTools::momentsGrad3D( ArrayTypeOut &vals ,
00932                                                const ArrayTypeData &data ,
00933                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
00934                                                const Array<RCP<ArrayTypeBasis> > &basisDVals ,
00935                                                const Array<RCP<ArrayTypeWeights> > &wts )
00936   {
00937 
00938     const int numBfx = basisVals[0]->dimension(0);
00939     const int numBfy = basisVals[1]->dimension(0);
00940     const int numBfz = basisVals[2]->dimension(0);
00941     const int numPtsx = basisVals[0]->dimension(1);
00942     const int numPtsy = basisVals[1]->dimension(1);
00943     const int numPtsz = basisVals[2]->dimension(1);
00944     const int numCells = vals.dimension(0);
00945     const int numFields = vals.dimension(1);
00946     const ArrayTypeBasis &Phix = *basisVals[0];
00947     const ArrayTypeBasis &Phiy = *basisVals[1];
00948     const ArrayTypeBasis &Phiz = *basisVals[2];
00949     const ArrayTypeBasis &DPhix = *basisDVals[0];
00950     const ArrayTypeBasis &DPhiy = *basisDVals[1];
00951     const ArrayTypeBasis &DPhiz = *basisDVals[2];
00952     const ArrayTypeWeights &wtsx = *wts[0];
00953     const ArrayTypeWeights &wtsy = *wts[1];
00954     const ArrayTypeWeights &wtsz = *wts[2];
00955 
00956     FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy,3);
00957     FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz,3);
00958   
00959     // Xi phase
00960     for (int f=0;f<numFields;f++)
00961       {
00962         for (int c=0;c<numCells;c++)
00963           {
00964             for (int i=0;i<numBfx;i++)
00965               {
00966                 for (int n=0;n<numPtsz;n++)
00967                   {
00968                     for (int m=0;m<numPtsy;m++)
00969                       {
00970                         for (int l=0;l<numPtsx;l++)
00971                           {
00972                             const int dataIdx = n * numPtsy * numPtsx + m * numPtsx + l;
00973                             Xi(c,i,n,m,0) += wtsx(l) * DPhix(i,l,0) * data(c,f,dataIdx);
00974                             Xi(c,i,n,m,1) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx);
00975                             Xi(c,i,n,m,2) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx);
00976                           }
00977                       }
00978                   }
00979               }
00980           }
00981 
00982         // Theta phase
00983         for (int c=0;c<numCells;c++)
00984           {
00985             for (int j=0;j<numBfy;j++)
00986               {
00987                 for (int i=0;i<numBfx;i++)
00988                   {
00989                     for (int n=0;n<numPtsz;n++)
00990                       {
00991                         for (int m=0;m<numPtsy;m++)
00992                           {
00993                             Theta(c,j,i,n,0) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,0);
00994                             Theta(c,j,i,n,1) += wtsy(j) * DPhiy(j,m,0) * Xi(c,i,n,m,1);
00995                             Theta(c,j,i,n,2) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,2);
00996                           }
00997                       }
00998                   }
00999               }
01000           }
01001 
01002         // last phase
01003         for (int c=0;c<numCells;c++)
01004           {
01005             for (int k=0;k<numBfz;k++)
01006               {
01007                 for (int j=0;j<numBfy;j++)
01008                   {
01009                     for (int i=0;i<numBfx;i++)
01010                       {
01011                         const int momIdx = k * numBfx * numBfy + j * numBfx + i;
01012                         for (int n=0;n<numPtsz;n++)
01013                           {
01014                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,0) * Phiz(k,n);
01015                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,1) * Phiz(k,n);
01016                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,2) * DPhiz(k,n,0);
01017                           }
01018                       }
01019                   }
01020               }
01021           }
01022       }
01023 
01024 }
01025 
01026 
01027 
01028 } // end namespace Intrepid