Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/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::evaluateCollocated( ArrayTypeOut &vals ,
00109                                                     const ArrayTypeCoeffs &coeffs ,
00110                                                     const Array<RCP<ArrayTypeBasis> > &bases )
00111   {
00112     const unsigned space_dim = bases.size();
00113 
00114 #ifdef HAVE_INTREPID_DEBUG
00115     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00116                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): vals must be rank 3 array." );
00117 
00118     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00119                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): coeffs must be rank 3 array." );
00120 
00121     for (unsigned d=0;d<space_dim;d++)
00122       {
00123         TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
00124                                     ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): each tabulated basis must be rank 2 array." );
00125       }
00126 
00127     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00128                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of cells for vals and coeffs must match." );
00129 
00130     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00131                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of fields for vals and coeffs must match." );
00132     
00133     int product_of_basis_dimensions = 1;
00134     int product_of_basis_points = 1;
00135     for (unsigned d=0;d<space_dim;d++)
00136       {
00137         product_of_basis_dimensions *= bases[d]->dimension(0);
00138         product_of_basis_points *= bases[d]->dimension(1);
00139       }
00140     
00141     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00142                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of points in vals and bases ." );
00143     
00144     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00145                                 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of basis functions in coeffs and bases ." );
00146 #endif    
00147     
00148 
00149     switch (space_dim)
00150       {
00151       case 2:
00152         evaluateCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00153         break;
00154       case 3:
00155         evaluateCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00156         break;
00157       }
00158 
00159   }
00160 
00161 //   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00162 //         class ArrayTypeBasis>
00163 //   void TensorProductSpaceTools::evaluateCollocated( ArrayTypeOut &vals ,
00164 //                                                  const ArrayTypeCoeffs &coeffs ,
00165 //                                                  const Array<Array<RCP<ArrayTypeBasis> > > &bases )
00166 //   {
00167 //     const unsigned space_dim = bases.size();
00168 
00169 // #ifdef HAVE_INTREPID_DEBUG
00170 //     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
00171 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): vals must be rank 3 array." );
00172 
00173 //     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00174 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): coeffs must be rank 3 array." );
00175 
00176 //     for (unsigned d0=0;d0<bases.size();d0++)
00177 //       {
00178 //      for (unsigned d1=1;d1<space_dim;d1++)
00179 //        {
00180 //          TEUCHOS_TEST_FOR_EXCEPTION( bases[d0][d1]->rank() != 2 , std::invalid_argument ,
00181 //                                      ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): each tabulated basis must be rank 2 array." );
00182 //        }
00183 //       }
00184 
00185 //     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00186 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of cells for vals and coeffs must match." );
00187 
00188 //     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00189 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of fields for vals and coeffs must match." );
00190     
00191 //     int product_of_basis_dimensions = 1;
00192 //     int product_of_basis_points = 1;
00193 //     for (unsigned d=0;d<space_dim;d++)
00194 //       {
00195 //      product_of_basis_dimensions *= bases[0][d]->dimension(0);
00196 //      product_of_basis_points *= bases[0][d]->dimension(1);
00197 //       }
00198     
00199 //     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00200 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of points in vals and bases ." );
00201     
00202 //     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00203 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of basis functions in coeffs and bases ." );
00204 // #endif    
00205 
00206 //     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank(3) != bases.size() , std::invalid_argument ,
00207 //                              ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): wrong dimensions for vals");
00208     
00209 
00210 //     switch (space_dim)
00211 //       {
00212 //       case 2:
00213 //      evaluateCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00214 //      break;
00215 //       case 3:
00216 //      evaluateCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
00217 //      break;
00218 //       }
00219 
00220 //   }
00221 
00222   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00223            class ArrayTypeBasis>
00224   void TensorProductSpaceTools::evaluateGradient( ArrayTypeOut &vals ,
00225                                                   const ArrayTypeCoeffs &coeffs ,
00226                                                   const Array<RCP<ArrayTypeBasis> > &bases ,
00227                                                   const Array<RCP<ArrayTypeBasis> > &Dbases )
00228   {
00229     const unsigned space_dim = bases.size();
00230 
00231 #ifdef HAVE_INTREPID_DEBUG
00232     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
00233                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00234 
00235     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00236                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00237 
00238     TEUCHOS_TEST_FOR_EXCEPTION( Dbases.size() != (int)space_dim , std::invalid_argument ,
00239                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases and Dbases must be same size." );
00240    
00241     for (unsigned d=0;d<space_dim;d++)
00242       {
00243         TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
00244                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00245 
00246         TEUCHOS_TEST_FOR_EXCEPTION( Dbases[d]->rank() != 3 , std::invalid_argument ,
00247                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 3 array." );
00248       }
00249 
00250     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00251                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00252 
00253     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00254                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00255     
00256     int product_of_basis_dimensions = 1;
00257     int product_of_basis_points = 1;
00258     for (unsigned d=0;d<space_dim;d++)
00259       {
00260         product_of_basis_dimensions *= bases[d]->dimension(0);
00261         product_of_basis_points *= bases[d]->dimension(1);
00262       }
00263     
00264     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00265                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00266     
00267     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00268                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
00269 
00270     for (unsigned d=0;d<space_dim;d++)
00271       {
00272         for (unsigned i=0;i<2;i++)
00273           {
00274             TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->dimension(i) != Dbases[d]->dimension(i) , std::invalid_argument ,
00275                                         ">>>ERROR (TensorProductSpaceTools::evaluate): bases and Dbases have incompatible shape." );
00276           }
00277       }
00278 #endif    
00279 
00280     switch (space_dim)
00281       {
00282       case 2:
00283         TensorProductSpaceTools::evaluateGradient2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00284         break;
00285       case 3:
00286         TensorProductSpaceTools::evaluateGradient3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00287         break;
00288       }
00289 
00290   }
00291 
00292   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00293            class ArrayTypeBasis>
00294   void TensorProductSpaceTools::evaluateGradientCollocated( ArrayTypeOut &vals ,
00295                                                             const ArrayTypeCoeffs &coeffs ,
00296                                                             const Array<RCP<ArrayTypeBasis> > &bases ,
00297                                                             const Array<RCP<ArrayTypeBasis> > &Dbases )
00298   {
00299     const unsigned space_dim = bases.size();
00300 
00301 #ifdef HAVE_INTREPID_DEBUG
00302     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
00303                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00304 
00305     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
00306                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00307 
00308     TEUCHOS_TEST_FOR_EXCEPTION( Dbases.size() != (int)space_dim , std::invalid_argument ,
00309                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases and Dbases must be same size." );
00310    
00311     for (unsigned d=0;d<space_dim;d++)
00312       {
00313         TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
00314                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00315 
00316         TEUCHOS_TEST_FOR_EXCEPTION( Dbases[d]->rank() != 3 , std::invalid_argument ,
00317                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 3 array." );
00318       }
00319 
00320     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
00321                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00322 
00323     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
00324                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00325     
00326     int product_of_basis_dimensions = 1;
00327     int product_of_basis_points = 1;
00328     for (unsigned d=0;d<space_dim;d++)
00329       {
00330         product_of_basis_dimensions *= bases[d]->dimension(0);
00331         product_of_basis_points *= bases[d]->dimension(1);
00332       }
00333     
00334     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
00335                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00336     
00337     TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00338                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
00339 
00340     for (unsigned d=0;d<space_dim;d++)
00341       {
00342         for (unsigned i=0;i<2;i++)
00343           {
00344             TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->dimension(i) != Dbases[d]->dimension(i) , std::invalid_argument ,
00345                                         ">>>ERROR (TensorProductSpaceTools::evaluate): bases and Dbases have incompatible shape." );
00346           }
00347       }
00348 #endif    
00349 
00350     switch (space_dim)
00351       {
00352       case 2:
00353         evaluateGradientCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00354         break;
00355       case 3:
00356         evaluateGradientCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
00357         break;
00358       }
00359 
00360   }
00361 
00362   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00363            class ArrayTypeBasis, class ArrayTypeWeights>
00364   void TensorProductSpaceTools::moments( ArrayTypeOut &vals ,
00365                                          const ArrayTypeData &data ,
00366                                          const Array<RCP<ArrayTypeBasis> > &basisVals ,
00367                                          const Array<RCP<ArrayTypeWeights> > &wts )
00368   {
00369     const unsigned space_dim = basisVals.size();
00370 
00371 #ifdef HAVE_INTREPID_DEBUG
00372     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00373                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00374 
00375     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 3 , std::invalid_argument ,
00376                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00377 
00378     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00379                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00380 
00381     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00382                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00383    
00384     for (unsigned d=0;d<space_dim;d++)
00385       {
00386         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00387                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00388 
00389         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00390                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00391       }
00392 
00393     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00394                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00395 
00396     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != data.dimension(1) , std::invalid_argument,
00397                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00398     
00399     int product_of_basis_dimensions = 1;
00400     int product_of_basis_points = 1;
00401     for (unsigned d=0;d<space_dim;d++)
00402       {
00403         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00404         product_of_basis_points *= basisVals[d]->dimension(1);
00405       }
00406     
00407     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00408                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00409     
00410     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
00411                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00412 
00413 #endif    
00414 
00415     switch (space_dim)
00416       {
00417       case 2:
00418         moments2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00419         break;
00420       case 3:
00421         moments3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00422         break;
00423       }
00424 
00425   }
00426 
00427   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00428            class ArrayTypeBasis, class ArrayTypeWeights>
00429   void TensorProductSpaceTools::momentsCollocated( ArrayTypeOut &vals ,
00430                                                    const ArrayTypeData &data ,
00431                                                    const Array<RCP<ArrayTypeBasis> > &basisVals ,
00432                                                    const Array<RCP<ArrayTypeWeights> > &wts )
00433   {
00434     const unsigned space_dim = basisVals.size();
00435 
00436 #ifdef HAVE_INTREPID_DEBUG
00437     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00438                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00439 
00440     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 3 , std::invalid_argument ,
00441                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
00442 
00443     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00444                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00445 
00446     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00447                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00448    
00449     for (unsigned d=0;d<space_dim;d++)
00450       {
00451         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00452                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00453 
00454         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00455                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00456       }
00457 
00458     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00459                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00460 
00461     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != data.dimension(1) , std::invalid_argument,
00462                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
00463     
00464     int product_of_basis_dimensions = 1;
00465     int product_of_basis_points = 1;
00466     for (unsigned d=0;d<space_dim;d++)
00467       {
00468         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00469         product_of_basis_points *= basisVals[d]->dimension(1);
00470       }
00471     
00472     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00473                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00474     
00475     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
00476                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00477 
00478 #endif    
00479 
00480     switch (space_dim)
00481       {
00482       case 2:
00483         moments2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00484         break;
00485       case 3:
00486         moments3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
00487         break;
00488       }
00489 
00490   }
00491 
00492   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00493            class ArrayTypeBasis, class ArrayTypeWeights>
00494   void TensorProductSpaceTools::momentsGrad( ArrayTypeOut &vals ,
00495                                              const ArrayTypeData &data ,
00496                                              const Array<RCP<ArrayTypeBasis> > &basisVals ,
00497                                              const Array<RCP<ArrayTypeBasis> > &basisDVals ,
00498                                              const Array<RCP<ArrayTypeWeights> > &wts )
00499   {
00500     const unsigned space_dim = basisVals.size();
00501 
00502 #ifdef HAVE_INTREPID_DEBUG
00503     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00504                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00505 
00506     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 4 , std::invalid_argument ,
00507                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 4 array." );
00508 
00509     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00510                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00511 
00512     TEUCHOS_TEST_FOR_EXCEPTION( basisDVals.size() != (int)space_dim , std::invalid_argument ,
00513                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00514 
00515     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00516                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00517    
00518     for (unsigned d=0;d<space_dim;d++)
00519       {
00520         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00521                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00522 
00523         TEUCHOS_TEST_FOR_EXCEPTION( basisDVals[d]->rank() != 3 , std::invalid_argument ,
00524                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated derivative basis must be rank 3 array." );
00525 
00526         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00527                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00528       }
00529 
00530     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00531                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00532     
00533     int product_of_basis_dimensions = 1;
00534     int product_of_basis_points = 1;
00535     for (unsigned d=0;d<space_dim;d++)
00536       {
00537         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00538         product_of_basis_points *= basisVals[d]->dimension(1);
00539       }
00540     
00541     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00542                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00543     
00544     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
00545                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00546 
00547 #endif    
00548 
00549     switch (space_dim)
00550       {
00551       case 2:
00552         momentsGrad2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
00553         break;
00554       case 3:
00555         momentsGrad3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
00556         break;
00557       }
00558 
00559   }
00560 
00561   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
00562            class ArrayTypeBasis, class ArrayTypeWeights>
00563   void TensorProductSpaceTools::momentsGradCollocated( ArrayTypeOut &vals ,
00564                                              const ArrayTypeData &data ,
00565                                              const Array<RCP<ArrayTypeBasis> > &basisVals ,
00566                                              const Array<RCP<ArrayTypeBasis> > &basisDVals ,
00567                                              const Array<RCP<ArrayTypeWeights> > &wts )
00568   {
00569     const unsigned space_dim = basisVals.size();
00570 
00571 #ifdef HAVE_INTREPID_DEBUG
00572     TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
00573                                 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
00574 
00575     TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 4 , std::invalid_argument ,
00576                                 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 4 array." );
00577 
00578     TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
00579                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00580 
00581     TEUCHOS_TEST_FOR_EXCEPTION( basisDVals.size() != (int)space_dim , std::invalid_argument ,
00582                                 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
00583 
00584     TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
00585                                 ">>> ERROR (TensorProductSpaceTools::evaluate):  quadrature weights ill-sized." );
00586    
00587     for (unsigned d=0;d<space_dim;d++)
00588       {
00589         TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
00590                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
00591 
00592         TEUCHOS_TEST_FOR_EXCEPTION( basisDVals[d]->rank() != 3 , std::invalid_argument ,
00593                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated derivative basis must be rank 3 array." );
00594 
00595         TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
00596                                     ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
00597       }
00598 
00599     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
00600                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
00601     
00602     int product_of_basis_dimensions = 1;
00603     int product_of_basis_points = 1;
00604     for (unsigned d=0;d<space_dim;d++)
00605       {
00606         product_of_basis_dimensions *= basisVals[d]->dimension(0);
00607         product_of_basis_points *= basisVals[d]->dimension(1);
00608       }
00609     
00610     TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
00611                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
00612     
00613     TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
00614                                 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
00615 
00616 #endif    
00617 
00618     switch (space_dim)
00619       {
00620       case 2:
00621         momentsGradCollocated2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
00622         break;
00623       case 3:
00624         momentsGradCollocated3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
00625         break;
00626       }
00627 
00628   }
00629 
00630 
00631   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00632            class ArrayTypeBasis>
00633   void TensorProductSpaceTools::evaluate2D( ArrayTypeOut &vals ,
00634                                             const ArrayTypeCoeffs &coeffs ,
00635                                             const Array<RCP<ArrayTypeBasis> > &bases )
00636   {
00637     const int numBfx = bases[0]->dimension(0);
00638     const int numBfy = bases[1]->dimension(0);
00639     
00640     const int numPtsx = bases[0]->dimension(1);
00641     const int numPtsy = bases[1]->dimension(1);
00642     
00643     const int numCells = vals.dimension(0);
00644     const int numFields = vals.dimension(1);
00645 
00646     const ArrayTypeBasis &Phix = *bases[0];
00647     const ArrayTypeBasis &Phiy = *bases[1];
00648     
00649     FieldContainer<double> Xi(numCells,numBfy,numPtsx);
00650     
00651     // sum factorization step
00652 
00653     for (int f=0;f<numFields;f++)
00654       {
00655         for (int cell=0;cell<numCells;cell++)
00656           {
00657             for (int j=0;j<numBfy;j++)
00658               {
00659                 for (int i=0;i<numBfx;i++)
00660                   {
00661                     const int I = j * numBfx + i;
00662                     for (int k=0;k<numPtsx;k++)
00663                       {
00664                         Xi(cell,j,k) += coeffs(cell,f,I) * Phix(i,k);
00665                       }
00666                   }
00667               }
00668           }
00669         
00670         for (int cell=0;cell<numCells;cell++)
00671           {
00672             for (int kpty=0;kpty<numPtsy;kpty++)
00673               {
00674                 for (int kptx=0;kptx<numPtsx;kptx++)
00675                   {
00676                     vals(cell,f,kptx+numPtsx*kpty) = 0.0;
00677                   }
00678               }
00679             
00680             // evaluation step
00681             for (int kpty=0;kpty<numPtsy;kpty++)
00682               {
00683                 for (int kptx=0;kptx<numPtsx;kptx++)
00684                   {
00685                     const int I=kptx+numPtsx*kpty;
00686                     
00687                     for (int j=0;j<numBfy;j++)
00688                       {
00689                         vals(cell,f,I) += Phiy(j,kpty) * Xi(cell,j,kptx);
00690                       }
00691                   }
00692               }
00693           }
00694       }
00695     
00696     return;
00697   }
00698 
00699   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00700            class ArrayTypeBasis>
00701   void TensorProductSpaceTools::evaluateCollocated2D( ArrayTypeOut &vals ,
00702                                                       const ArrayTypeCoeffs &coeffs ,
00703                                                       const Array<RCP<ArrayTypeBasis> > &bases )
00704   {
00705     // just copy coeffs to vals!
00706 
00707     const int numBfx = bases[0]->dimension(0);
00708     const int numBfy = bases[1]->dimension(0);
00709     
00710     const int numCells = vals.dimension(0);
00711     const int numFields = vals.dimension(1);
00712 
00713     for (int cell=0;cell<numCells;cell++)
00714       {
00715         for (int f=0;f<numFields;f++)
00716           {
00717             for (int j=0;j<numBfy;j++)
00718               {
00719                 for (int i=0;i<numBfx;i++)
00720                   {
00721                     const int I = j * numBfx + i;
00722                     vals(cell,f,I) = coeffs(cell,f,I);
00723                   }
00724               }
00725           }
00726       }
00727 
00728     return;
00729   }
00730 
00731   // template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00732   //       class ArrayTypeBasis>
00733   // void TensorProductSpaceTools::evaluateCollocated2D( ArrayTypeOut &vals ,
00734   //                                                  const ArrayTypeCoeffs &coeffs ,
00735   //                                                  const Array<Array<RCP<ArrayTypeBasis> > > &bases )
00736   // {
00737   //   // just copy coeffs to vals!
00738   //   const int numCells = vals.dimension(0);
00739   //   const int numFields = vals.dimension(1);
00740     
00741   //   const int numBfx = bases[comp][0]->dimension(0);
00742   //   const int numBfy = bases[comp][1]->dimension(0);
00743     
00744     
00745   //   for (int cell=0;cell<numCells;cell++)
00746   //     {
00747   //    for (int f=0;f<numFields;f++)
00748   //      {
00749   //        for (int j=0;j<numBfy;j++)
00750   //          {
00751   //            for (int i=0;i<numBfx;i++)
00752   //              {                     const int I = j * numBfx + i;
00753   //                vals(cell,f,I,comp) = coeffs(cell,f,I);
00754   //              }
00755   //          }
00756   //      }
00757   //     }
00758 
00759   //   return;
00760   // }
00761 
00762   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00763            class ArrayTypeBasis>
00764   void TensorProductSpaceTools::evaluate3D( ArrayTypeOut &vals ,
00765                                             const ArrayTypeCoeffs &coeffs ,
00766                                             const Array<RCP<ArrayTypeBasis> > &bases )  
00767 
00768   {
00769     // checks to do:
00770     // vals point dimension is product of sizes of point arrays
00771     // vals
00772     
00773     const int numBfx = bases[0]->dimension(0);
00774     const int numBfy = bases[1]->dimension(0);
00775     const int numBfz = bases[2]->dimension(0);
00776     
00777     const int numPtsx = bases[0]->dimension(1);
00778     const int numPtsy = bases[1]->dimension(1);
00779     const int numPtsz = bases[2]->dimension(1);
00780     
00781     const int numCells = vals.dimension(0);
00782     const int numFields = vals.dimension(1);
00783     
00784     const ArrayTypeBasis &Phix = *bases[0];
00785     const ArrayTypeBasis &Phiy = *bases[1];
00786     const FieldContainer<double> &Phiz = *bases[2];
00787 
00788     FieldContainer<double> Xi(numCells, 
00789                               numBfz, numBfy , numPtsx);
00790     
00791     FieldContainer<double> Theta(numCells, 
00792                                  numBfz , numPtsy, numPtsx );
00793     
00794     
00795     for (int f=0;f<numFields;f++)
00796       {
00797 
00798         // Xi section
00799         for (int c=0;c<numCells;c++)
00800           {
00801             for (int k=0;k<numBfz;k++)
00802               {
00803                 for (int j=0;j<numBfy;j++)
00804                   {
00805                     for (int l=0;l<numPtsx;l++)
00806                       {
00807                         for (int i=0;i<numBfx;i++)
00808                           {
00809                             const int coeffIndex = k*numBfy*numBfx + j * numBfx + i;
00810                             Xi(c,k,j,l) += Phix(i,l) * coeffs(c,f,coeffIndex);
00811                           }
00812                       }
00813                   }
00814               }
00815           }
00816     
00817         // Theta section
00818         for (int c=0;c<numCells;c++)
00819           {
00820             for (int k=0;k<numBfz;k++)
00821               {
00822                 for (int m=0;m<numPtsy;m++)
00823                   {
00824                     for (int l=0;l<numPtsx;l++)
00825                       {
00826                         for (int j=0;j<numBfy;j++)
00827                           {
00828                             Theta(c,k,m,l) += Phiy(j,m) * Xi(c,k,j,l);
00829                           }
00830                       }
00831                   }
00832               }
00833           }
00834     
00835         // final section
00836         for (int c=0;c<numCells;c++)
00837           {
00838             for (int n=0;n<numPtsz;n++)
00839               {
00840                 for (int m=0;m<numPtsy;m++)
00841                   {
00842                     for (int l=0;l<numPtsx;l++)
00843                       {
00844                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) = 0.0;
00845                         for (int k=0;k<numBfz;k++)
00846                           {
00847                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) += Theta(c,k,m,l) * Phiz(k,n);
00848                           }
00849                       }
00850                   }
00851               }
00852           }
00853       }
00854 
00855     return;
00856     
00857   }
00858 
00859   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00860            class ArrayTypeBasis>
00861   void TensorProductSpaceTools::evaluateCollocated3D( ArrayTypeOut &vals ,
00862                                                       const ArrayTypeCoeffs &coeffs ,
00863                                                       const Array<RCP<ArrayTypeBasis> > &bases )  
00864 
00865   {
00866     // copy coeffs to vals
00867     
00868     const int numBfx = bases[0]->dimension(0);
00869     const int numBfy = bases[1]->dimension(0);
00870     const int numBfz = bases[2]->dimension(0);
00871     
00872     const int numCells = vals.dimension(0);
00873     const int numFields = vals.dimension(1);
00874     
00875     for (int cell=0;cell<numCells;cell++)
00876       {
00877         for (int field=0;field<numFields;field++)
00878           {     
00879             for (int k=0;k<numBfz;k++)
00880               {
00881                 for (int j=0;j<numBfy;j++)
00882                   {
00883                     for (int i=0;i<numBfx;i++)
00884                       {
00885                         const int I = k*numBfy*numBfx + j * numBfx + i;
00886                         vals(cell,field,I) = coeffs(cell,field,I);
00887                       }
00888                   }
00889               }
00890           }
00891       }
00892 
00893     return;
00894     
00895   }
00896 
00897 
00898   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
00899            class ArrayTypeBasis>
00900   void TensorProductSpaceTools::evaluateGradient2D( ArrayTypeOut &vals ,
00901                                                     const ArrayTypeCoeffs &coeffs ,
00902                                                     const Array<RCP<ArrayTypeBasis> > &bases ,
00903                                                     const Array<RCP<ArrayTypeBasis> > &Dbases )
00904   {
00905     
00906     const int numBfx = bases[0]->dimension(0);
00907     const int numBfy = bases[1]->dimension(0);
00908     const int numPtsx = bases[0]->dimension(1);
00909     const int numPtsy = bases[1]->dimension(1);
00910     const int numCells = vals.dimension(0);
00911     const int numFields = vals.dimension(1);
00912     const ArrayTypeBasis &Phix = *bases[0];
00913     const ArrayTypeBasis &Phiy = *bases[1];
00914     const ArrayTypeBasis &DPhix = *Dbases[0];
00915     const ArrayTypeBasis &DPhiy = *Dbases[1];
00916     
00917     FieldContainer<double> Xi(numBfy,numPtsx);
00918 
00919     for (int f=0;f<numFields;f++)
00920       {
00921     
00922         for (int cell=0;cell<numCells;cell++)
00923           {
00924             // x derivative
00925         
00926             // sum factorization step
00927             for (int j=0;j<numBfy;j++)
00928               {
00929                 for (int k=0;k<numPtsx;k++)
00930                   {
00931                     Xi(j,k) = 0.0;
00932                   }
00933               }
00934         
00935             for (int j=0;j<numBfy;j++)
00936               {
00937                 for (int i=0;i<numBfx;i++)
00938                   {
00939                     const int I = j * numBfx + i;
00940                     for (int k=0;k<numPtsx;k++)
00941                       {
00942                         Xi(j,k) += coeffs(cell,f,I) * DPhix(i,k,0);
00943                       }
00944                   }
00945               }
00946         
00947             for (int kpty=0;kpty<numPtsy;kpty++)
00948               {
00949                 for (int kptx=0;kptx<numPtsx;kptx++)
00950                   {
00951                     vals(cell,f,kptx+numPtsx*kpty,0) = 0.0;
00952                   }
00953               }
00954         
00955             // evaluation step
00956             for (int kpty=0;kpty<numPtsy;kpty++)
00957               {
00958                 for (int kptx=0;kptx<numPtsx;kptx++)
00959                   {
00960                     const int I=kptx+numPtsx*kpty;
00961                 
00962                     for (int j=0;j<numBfy;j++)
00963                       {
00964                         vals(cell,f,I,0) += Phiy(j,kpty) * Xi(j,kptx);
00965                       }
00966                   }
00967               }
00968         
00969             // y derivative
00970         
00971             // sum factorization step
00972             for (int j=0;j<numBfy;j++)
00973               {
00974                 for (int k=0;k<numPtsx;k++)
00975                   {
00976                     Xi(j,k) = 0.0;
00977                   }
00978               }
00979         
00980             for (int j=0;j<numBfy;j++)
00981               {
00982                 for (int i=0;i<numBfx;i++)
00983                   {
00984                     const int I = j * numBfx + i;
00985                     for (int k=0;k<numPtsx;k++)
00986                       {
00987                         Xi(j,k) += coeffs(cell,f,I) * Phix(i,k);
00988                       }
00989                   }
00990               }
00991         
00992             for (int kpty=0;kpty<numPtsy;kpty++)
00993               {
00994                 for (int kptx=0;kptx<numPtsx;kptx++)
00995                   {
00996                     vals(cell,f,kptx+numPtsx*kpty,1) = 0.0;
00997                   }
00998               }
00999         
01000             // evaluation step
01001             for (int kpty=0;kpty<numPtsy;kpty++)
01002               {
01003                 for (int kptx=0;kptx<numPtsx;kptx++)
01004                   {
01005                     const int I=kptx+numPtsx*kpty;
01006                 
01007                     for (int j=0;j<numBfy;j++)
01008                       {
01009                         vals(cell,f,I,1) += DPhiy(j,kpty,0) * Xi(j,kptx);
01010                       }
01011                   }
01012               }
01013           }
01014       }    
01015     return;
01016   }
01017 
01018   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
01019            class ArrayTypeBasis>
01020   void TensorProductSpaceTools::evaluateGradientCollocated2D( ArrayTypeOut &vals ,
01021                                                               const ArrayTypeCoeffs &coeffs ,
01022                                                               const Array<RCP<ArrayTypeBasis> > &bases ,
01023                                                               const Array<RCP<ArrayTypeBasis> > &Dbases )
01024   {
01025     
01026     const int numBfx = bases[0]->dimension(0);
01027     const int numBfy = bases[1]->dimension(0);
01028     const int numPtsx = bases[0]->dimension(1);
01029     const int numPtsy = bases[1]->dimension(1);
01030     const int numCells = vals.dimension(0);
01031     const int numFields = vals.dimension(1);
01032     const ArrayTypeBasis &DPhix = *Dbases[0];
01033     const ArrayTypeBasis &DPhiy = *Dbases[1];
01034     
01035     for (int cell=0;cell<numCells;cell++)
01036       {
01037         for (int field=0;field<numFields;field++)
01038           {
01039             for (int j=0;j<numPtsy;j++)
01040               {
01041                 for (int i=0;i<numPtsx;i++)
01042                   {
01043                     const int I = j * numPtsx + i;
01044                     vals(cell,field,I,0) = 0.0;
01045                     vals(cell,field,I,1) = 0.0;
01046                   }
01047               }
01048           }
01049       }
01050     
01051     // x derivative
01052     for (int cell=0;cell<numCells;cell++)
01053       {
01054         for (int field=0;field<numFields;field++)
01055           {
01056             for (int j=0;j<numPtsy;j++)
01057               {
01058                 for (int i=0;i<numPtsx;i++)
01059                   {
01060                     const int I = j * numPtsx + i;
01061                     for (int ell=0;ell<numBfx;ell++)
01062                       {
01063                         const int Itmp = j*numPtsx + ell;
01064                         vals(cell,field,I,0) +=  coeffs(cell,field,Itmp) * DPhix( ell , i );
01065                       }
01066                   }
01067               }
01068           }
01069       }
01070 
01071     // y derivative
01072     for (int cell=0;cell<numCells;cell++)
01073       {
01074         for (int field=0;field<numFields;field++)
01075           {
01076             for (int j=0;j<numPtsy;j++)
01077               {
01078                 for (int i=0;i<numPtsx;i++)
01079                   {
01080                     const int I = j * numPtsx + i;
01081                     for (int m=0;m<numBfy;m++)
01082                       {
01083                         const int Itmp = m*numPtsx + i;
01084                         vals(cell,field,I,1) +=  coeffs(cell,field,Itmp) * DPhiy( m , j );
01085                       }
01086                   }
01087               }
01088           }
01089       }
01090 
01091   }
01092 
01093   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
01094            class ArrayTypeBasis>
01095   void TensorProductSpaceTools::evaluateGradient3D( ArrayTypeOut &vals ,
01096                                                     const ArrayTypeCoeffs &coeffs ,
01097                                                     const Array<RCP<ArrayTypeBasis> > &bases ,
01098                                                     const Array<RCP<ArrayTypeBasis> > &Dbases )
01099   
01100   {
01101     // checks to do:
01102     // vals point dimension is product of sizes of point arrays
01103     // vals
01104 
01105     const int numBfx = bases[0]->dimension(0);
01106     const int numBfy = bases[1]->dimension(0);
01107     const int numBfz = bases[2]->dimension(0);
01108     const int numPtsx = bases[0]->dimension(1);
01109     const int numPtsy = bases[1]->dimension(1);
01110     const int numPtsz = bases[2]->dimension(1);
01111     const int numCells = vals.dimension(0);
01112     const int numFields = vals.dimension(1);
01113     const ArrayTypeBasis &Phix = *bases[0];
01114     const ArrayTypeBasis &Phiy = *bases[1];
01115     const ArrayTypeBasis &Phiz = *bases[2];
01116     const ArrayTypeBasis &DPhix = *Dbases[0];
01117     const ArrayTypeBasis &DPhiy = *Dbases[1];
01118     const ArrayTypeBasis &DPhiz = *Dbases[2];
01119 
01120     FieldContainer<double> Xi(numCells, 
01121                               numBfz, numBfy , numPtsx , 3) ;
01122 
01123     FieldContainer<double> Theta(numCells, 
01124                                  numBfz , numPtsy, numPtsx , 3);
01125   
01126     for (int f=0;f<numFields;f++)
01127       {
01128 
01129         // Xi section
01130         for (int c=0;c<numCells;c++)
01131           {
01132             for (int k=0;k<numBfz;k++)
01133               {
01134                 for (int j=0;j<numBfy;j++)
01135                   {
01136                     for (int l=0;l<numPtsx;l++)
01137                       {
01138                         for (int i=0;i<numBfx;i++)
01139                           {
01140                             const int coeffIndex = k*numBfz*numBfx + j * numBfx + i;
01141                             Xi(c,k,j,l,0) += DPhix(i,l,0) * coeffs(c,f,coeffIndex);
01142                             Xi(c,k,j,l,1) += Phix(i,l) * coeffs(c,f,coeffIndex);
01143                             Xi(c,k,j,l,2) += Phix(i,l) * coeffs(c,f,coeffIndex);
01144                           }
01145                       }
01146                   }
01147               }
01148           }
01149 
01150         // Theta section
01151         for (int c=0;c<numCells;c++)
01152           {
01153             for (int k=0;k<numBfz;k++)
01154               {
01155                 for (int m=0;m<numPtsy;m++)
01156                   {
01157                     for (int l=0;l<numPtsx;l++)
01158                       {
01159                         for (int j=0;j<numBfy;j++)
01160                           {
01161                             Theta(c,k,m,l,0) += Phiy(j,m) * Xi(c,k,j,l,0);
01162                             Theta(c,k,m,l,1) += DPhiy(j,m,0) * Xi(c,k,j,l,1);
01163                             Theta(c,k,m,l,2) += Phiy(j,m) * Xi(c,k,j,l,2);
01164                           }
01165                       }
01166                   }
01167               }
01168           }
01169 
01170         // final section
01171         for (int c=0;c<numCells;c++)
01172           {
01173             for (int n=0;n<numPtsz;n++)
01174               {
01175                 for (int m=0;m<numPtsy;m++)
01176                   {
01177                     for (int l=0;l<numPtsx;l++)
01178                       {
01179                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) = 0.0;
01180                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) = 0.0;
01181                         vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) = 0.0;
01182                         for (int k=0;k<numBfz;k++)
01183                           {
01184                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) += Theta(c,k,m,l,0) * Phiz(k,n);
01185                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) += Theta(c,k,m,l,1) * Phiz(k,n);
01186                             vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) += Theta(c,k,m,l,2) * DPhiz(k,n,0);
01187 
01188                           }
01189                       }
01190                   }
01191               }
01192           }
01193       }
01194     return;
01195   }
01196 
01197   template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
01198            class ArrayTypeBasis>
01199   void TensorProductSpaceTools::evaluateGradientCollocated3D( ArrayTypeOut &vals ,
01200                                                     const ArrayTypeCoeffs &coeffs ,
01201                                                     const Array<RCP<ArrayTypeBasis> > &bases ,
01202                                                     const Array<RCP<ArrayTypeBasis> > &Dbases )
01203   
01204   {
01205     const int numBfx = bases[0]->dimension(0);
01206     const int numBfy = bases[1]->dimension(0);
01207     const int numBfz = bases[2]->dimension(0);
01208     const int numPtsx = bases[0]->dimension(1);
01209     const int numPtsy = bases[1]->dimension(1);
01210     const int numPtsz = bases[2]->dimension(1);
01211     const int numCells = vals.dimension(0);
01212     const int numFields = vals.dimension(1);
01213     const ArrayTypeBasis &Phix = *bases[0];
01214     const ArrayTypeBasis &Phiy = *bases[1];
01215     const ArrayTypeBasis &Phiz = *bases[2];
01216     const ArrayTypeBasis &DPhix = *Dbases[0];
01217     const ArrayTypeBasis &DPhiy = *Dbases[1];
01218     const ArrayTypeBasis &DPhiz = *Dbases[2];
01219 
01220     for (int cell=0;cell<numCells;cell++)
01221       {
01222         for (int field=0;field<numFields;field++)
01223           {
01224             for (int k=0;k<numPtsz;k++)
01225               {
01226                 for (int j=0;j<numPtsy;j++)
01227                   {
01228                     for (int i=0;i<numPtsx;i++)
01229                       {
01230                         const int I = k * numPtsx * numPtsy + j * numPtsx + i;
01231                         vals(cell,field,I,0) = 0.0;
01232                         vals(cell,field,I,1) = 0.0;
01233                         vals(cell,field,I,2) = 0.0;
01234                       }
01235                   }
01236               }
01237           }
01238       }
01239     
01240     // x derivative
01241     for (int cell=0;cell<numCells;cell++)
01242       {
01243         for (int field=0;field<numFields;field++)
01244           {
01245             for (int k=0;k<numPtsz;k++)
01246               {
01247                 for (int j=0;j<numPtsy;j++)
01248                   {
01249                     for (int i=0;i<numPtsx;i++)
01250                       {
01251                         const int I = k * numPtsx * numPtsy + j * numPtsx + i;
01252                         for (int ell=0;ell<numBfx;ell++)
01253                           {
01254                             const int Itmp = k * numPtsx * numPtsy + j*numPtsx + ell;
01255                             vals(cell,field,I,0) +=  coeffs(cell,field,Itmp) * DPhix( ell , i );
01256                           }
01257                       }
01258                   }
01259               }
01260           }
01261       }
01262 
01263     // y derivative
01264     for (int cell=0;cell<numCells;cell++)
01265       {
01266         for (int field=0;field<numFields;field++)
01267           {
01268             for (int k=0;k<numPtsz;k++)
01269               {
01270                 for (int j=0;j<numPtsy;j++)
01271                   {
01272                     for (int i=0;i<numPtsx;i++)
01273                       {
01274                         const int I = k * numPtsx * numPtsy + j * numPtsx + i;
01275                         for (int m=0;m<numBfy;m++)
01276                           {
01277                             const int Itmp = k * numPtsx * numPtsy + m * numPtsx + i;
01278                             vals(cell,field,I,1) +=  coeffs(cell,field,Itmp) * DPhiy( m , j );
01279                           }
01280                       }
01281                   }
01282               }
01283           }
01284       }
01285 
01286 
01287     // z derivative
01288     for (int cell=0;cell<numCells;cell++)
01289       {
01290         for (int field=0;field<numFields;field++)
01291           {
01292             for (int k=0;k<numPtsz;k++)
01293               {
01294                 for (int j=0;j<numPtsy;j++)
01295                   {
01296                     for (int i=0;i<numPtsx;i++)
01297                       {
01298                         const int I = k * numPtsx * numPtsy + j * numPtsx + i;
01299                         for (int n=0;n<numBfz;n++)
01300                           {
01301                             const int Itmp = n * numPtsx * numPtsy + j * numPtsx + i;
01302                             vals(cell,field,I,2) +=  coeffs(cell,field,Itmp) * DPhiz( n , k );
01303                           }
01304                       }
01305                   }
01306               }
01307           }
01308       }
01309 
01310 
01311   
01312     return;
01313   }
01314 
01315   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01316            class ArrayTypeBasis, class ArrayTypeWeights>
01317   void TensorProductSpaceTools::moments2D( ArrayTypeOut &vals ,
01318                                            const ArrayTypeData &data ,
01319                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
01320                                            const Array<RCP<ArrayTypeWeights> > &wts )
01321   {
01322     const int numBfx = basisVals[0]->dimension(0);
01323     const int numBfy = basisVals[1]->dimension(0);
01324     const int numPtsx = basisVals[0]->dimension(1);
01325     const int numPtsy = basisVals[1]->dimension(1);
01326     const int numCells = vals.dimension(0);
01327     const int numFields = vals.dimension(1);
01328     const ArrayTypeBasis &Phix = *basisVals[0];
01329     const ArrayTypeBasis &Phiy = *basisVals[1];
01330     
01331     FieldContainer<double> Xi(numBfx,numPtsy);
01332     
01333     const ArrayTypeWeights &wtsx = *wts[0];
01334     const ArrayTypeWeights &wtsy = *wts[1];
01335 
01336     for (int f=0;f<numFields;f++)
01337       {
01338         for (int cell=0;cell<numCells;cell++)
01339           {
01340             // sum factorization step
01341             for (int i=0;i<numBfx;i++)
01342               {
01343                 for (int k=0;k<numPtsy;k++)
01344                   {
01345                     Xi(i,k) = 0.0;
01346                   }
01347               }
01348         
01349             for (int i=0;i<numBfx;i++)
01350               {
01351                 for (int l=0;l<numPtsy;l++)
01352                   {
01353                     for (int k=0;k<numPtsx;k++)
01354                       {
01355                         Xi(i,l) += wtsx(k) * data(cell,f,l*numPtsx+k) * Phix(i,k);
01356                       }
01357                   }
01358               }
01359         
01360             for (int j=0;j<numBfy;j++)
01361               {
01362                 for (int i=0;i<numBfx;i++)
01363                   {
01364                     vals(cell,f,j*numBfx+i) = 0.0;
01365                   }
01366               }
01367 
01368             // evaluate moments with sum factorization
01369             for (int j=0;j<numBfy;j++)
01370               {
01371                 for (int i=0;i<numBfx;i++)
01372                   {
01373                     for (int l=0;l<numPtsy;l++)
01374                       {
01375                         vals(cell,f,j*numBfx+i) += wtsy(l) * Phiy(j,l) * Xi(i,l);
01376                       }
01377                   }
01378               }
01379           }
01380       }
01381     return;
01382 
01383   }
01384 
01385   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01386            class ArrayTypeBasis, class ArrayTypeWeights>
01387   void TensorProductSpaceTools::momentsCollocated2D( ArrayTypeOut &vals ,
01388                                            const ArrayTypeData &data ,
01389                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
01390                                            const Array<RCP<ArrayTypeWeights> > &wts )
01391   {
01392     const int numBfx = basisVals[0]->dimension(0);
01393     const int numBfy = basisVals[1]->dimension(0);
01394     const int numPtsx = basisVals[0]->dimension(1);
01395     const int numPtsy = basisVals[1]->dimension(1);
01396     const int numCells = vals.dimension(0);
01397     const int numFields = vals.dimension(1);
01398     
01399     FieldContainer<double> Xi(numBfx,numPtsy);
01400     
01401     const ArrayTypeWeights &wtsx = *wts[0];
01402     const ArrayTypeWeights &wtsy = *wts[1];
01403 
01404     for (int cell=0;cell<numCells;cell++)
01405       {
01406         for (int field=0;field<numFields;field++)
01407           {
01408             for (int i=0;i<numBfx;i++)
01409               {
01410                 for (int j=0;j<numBfy;j++)
01411                   {
01412                     const int I = numBfy * i + j;
01413                     vals(cell,field,I) = wtsx(i) * wtsy(j) * data(cell,field,I); 
01414                   }
01415               }
01416           }
01417       }
01418 
01419     return;
01420 
01421   }
01422 
01423   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01424            class ArrayTypeBasis, class ArrayTypeWeights>
01425   void TensorProductSpaceTools::moments3D( ArrayTypeOut &vals ,
01426                                            const ArrayTypeData &data ,
01427                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
01428                                            const Array<RCP<ArrayTypeWeights> > &wts )
01429   {
01430     const int numBfx = basisVals[0]->dimension(0);
01431     const int numBfy = basisVals[1]->dimension(0);
01432     const int numBfz = basisVals[2]->dimension(0);
01433 
01434     const int numPtsx = basisVals[0]->dimension(1);
01435     const int numPtsy = basisVals[1]->dimension(1);
01436     const int numPtsz = basisVals[2]->dimension(1);
01437 
01438     const int numCells = vals.dimension(0);
01439     const int numFields = vals.dimension(1);
01440 
01441     const ArrayTypeBasis &Phix = *basisVals[0];
01442     const ArrayTypeBasis &Phiy = *basisVals[1];
01443     const ArrayTypeBasis &Phiz = *basisVals[2];
01444 
01445     const ArrayTypeWeights &Wtx = *wts[0];
01446     const ArrayTypeWeights &Wty = *wts[1];
01447     const ArrayTypeWeights &Wtz = *wts[2];
01448 
01449     FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy);
01450     FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz);
01451 
01452     for (int f=0;f<numFields;f++)
01453       {
01454         // Xi phase
01455         for (int c=0;c<numCells;c++)
01456           {
01457             for (int i=0;i<numBfx;i++)
01458               {
01459                 for (int n=0;n<numPtsz;n++)
01460                   {
01461                     for (int m=0;m<numPtsy;m++)
01462                       {
01463                         for (int l=0;l<numPtsx;l++)
01464                           {
01465                             Xi(c,i,n,m) += Wtx(l) * Phix(i,l) * data(c,f,n*numPtsy*numPtsx+m*numPtsx+l);
01466                           }
01467                       }
01468                   }
01469               }
01470           }
01471 
01472         // Theta phase
01473         for (int c=0;c<numCells;c++)
01474           {
01475             for (int j=0;j<numBfy;j++)
01476               {
01477                 for (int i=0;i<numBfx;i++)
01478                   {
01479                     for (int n=0;n<numPtsz;n++)
01480                       {
01481                         for (int m=0;m<numPtsy;m++)
01482                           {
01483                             Theta(c,j,i,n) += Wty(m) * Phiy(j,m) * Xi(c,i,n,m);
01484                           }
01485                       }
01486                   }
01487               }
01488           }
01489 
01490         // Final phase
01491         for (int c=0;c<numCells;c++)
01492           {
01493             for (int k=0;k<numBfz;k++)
01494               {
01495                 for (int j=0;j<numBfy;j++)
01496                   {
01497                     for (int i=0;i<numBfx;i++)
01498                       {
01499                         const int momIdx = k*numBfx*numBfy+j*numBfx+i;
01500                         for (int n=0;n<numPtsz;n++)
01501                           {
01502                             vals(c,f,momIdx) += Wtz(n) * Phiz(k,n) * Theta(c,j,i,n);
01503                           }
01504                       }
01505                   }
01506               }
01507           }
01508 
01509       }
01510     return;
01511 
01512   }
01513 
01514   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01515            class ArrayTypeBasis, class ArrayTypeWeights>
01516   void TensorProductSpaceTools::momentsCollocated3D( ArrayTypeOut &vals ,
01517                                            const ArrayTypeData &data ,
01518                                            const Array<RCP<ArrayTypeBasis> > &basisVals ,
01519                                            const Array<RCP<ArrayTypeWeights> > &wts )
01520   {
01521     const int numBfx = basisVals[0]->dimension(0);
01522     const int numBfy = basisVals[1]->dimension(0);
01523     const int numBfz = basisVals[2]->dimension(0);
01524 
01525     const int numPtsx = basisVals[0]->dimension(1);
01526     const int numPtsy = basisVals[1]->dimension(1);
01527     const int numPtsz = basisVals[2]->dimension(1);
01528 
01529     const int numCells = vals.dimension(0);
01530     const int numFields = vals.dimension(1);
01531     
01532     const ArrayTypeWeights &Wtx = *wts[0];
01533     const ArrayTypeWeights &Wty = *wts[1];
01534     const ArrayTypeWeights &Wtz = *wts[2];
01535 
01536     for (int cell=0;cell<numCells;cell++)
01537       {
01538         for (int field=0;field<numFields;field++)
01539           {
01540             for (int k=0;k<numBfz;k++)
01541               {
01542                 for (int j=0;j<numBfy;j++)
01543                   {
01544                     for (int i=0;i<numBfx;i++)
01545                       {
01546                         const int I = k * numBfy * numBfx + j * numBfx + i;
01547                         vals(cell,field,I) = Wtx(i) * Wty(j) * Wtz(k) * data(cell,field,I);
01548                       }
01549                   }
01550               }
01551           }
01552       }
01553 
01554     return;
01555 
01556   }
01557 
01558   // data is now (C,P,D)
01559   // want to compute the moments against the gradients of the basis 
01560   // functions.
01561 
01562   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01563            class ArrayTypeBasis, class ArrayTypeWeights>
01564   void TensorProductSpaceTools::momentsGrad2D( ArrayTypeOut &vals ,
01565                                                const ArrayTypeData &data ,
01566                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
01567                                                const Array<RCP<ArrayTypeBasis> > &Dbases ,
01568                                                const Array<RCP<ArrayTypeWeights> > &wts )
01569   {
01570     
01571     const int numBfx = basisVals[0]->dimension(0);
01572     const int numBfy = basisVals[1]->dimension(0);
01573     const int numPtsx = basisVals[0]->dimension(1);
01574     const int numPtsy = basisVals[1]->dimension(1);
01575     const int numCells = vals.dimension(0);
01576     const int numFields = vals.dimension(1);
01577     const ArrayTypeBasis &Phix = *basisVals[0];
01578     const ArrayTypeBasis &Phiy = *basisVals[1];
01579     const ArrayTypeBasis &DPhix = *Dbases[0];
01580     const ArrayTypeBasis &DPhiy = *Dbases[1];
01581     const ArrayTypeWeights &wtsx = *wts[0];
01582     const ArrayTypeWeights &wtsy = *wts[1];
01583     
01584     FieldContainer<double> Xi(numBfx,numPtsy,2);
01585 
01586     for (int f=0;f<numFields;f++)
01587       {
01588         // Xi phase
01589         for (int c=0;c<numCells;c++)
01590           {
01591             for (int i=0;i<numBfx;i++)
01592               {
01593                 for (int m=0;m<numPtsy;m++)
01594                   {
01595                     for (int l=0;l<numPtsx;l++)
01596                       {
01597                         Xi(i,m,0) += wtsx(l) * data(c,f,m*numPtsy+l,0) * DPhix(i,l);
01598                         Xi(i,m,1) += wtsx(l) * data(c,f,m*numPtsy+l,1) * Phix(i,l);
01599                       }
01600                   }
01601               }
01602           }
01603 
01604         // main phase
01605         for (int c=0;c<numCells;c++)
01606           {
01607             for (int j=0;j<numBfy;j++)
01608               {
01609                 for (int i=0;i<numBfx;i++)
01610                   {
01611                     const int bfIdx = j*numBfx+i;
01612                     vals(c,f,bfIdx) = 0.0;
01613                     for (int m=0;m<numPtsy;m++)
01614                       {
01615                         vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,0) * Phiy(j,m);
01616                         vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,1) * DPhiy(j,m);
01617                       }
01618                   }
01619               }
01620           }
01621       }
01622     return;
01623 
01624   }
01625 
01626   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01627            class ArrayTypeBasis, class ArrayTypeWeights>
01628   void TensorProductSpaceTools::momentsGradCollocated2D( ArrayTypeOut &vals ,
01629                                                const ArrayTypeData &data ,
01630                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
01631                                                const Array<RCP<ArrayTypeBasis> > &Dbases ,
01632                                                const Array<RCP<ArrayTypeWeights> > &wts )
01633   {
01634     
01635     const int numBfx = basisVals[0]->dimension(0);
01636     const int numBfy = basisVals[1]->dimension(0);
01637     const int numPtsx = basisVals[0]->dimension(1);
01638     const int numPtsy = basisVals[1]->dimension(1);
01639     const int numCells = vals.dimension(0);
01640     const int numFields = vals.dimension(1);
01641     const ArrayTypeBasis &Phix = *basisVals[0];
01642     const ArrayTypeBasis &Phiy = *basisVals[1];
01643     const ArrayTypeBasis &DPhix = *Dbases[0];
01644     const ArrayTypeBasis &DPhiy = *Dbases[1];
01645     const ArrayTypeWeights &wtsx = *wts[0];
01646     const ArrayTypeWeights &wtsy = *wts[1];
01647 
01648     for (int cell=0;cell<numCells;cell++)
01649       {
01650         for (int field=0;field<numFields;field++)
01651           {
01652             for (int j=0;j<numBfy;j++)
01653               {
01654                 for (int i=0;i<numBfx;i++)
01655                   {
01656                     const int I=j*numBfx+i;
01657                     vals(cell,field,I) = 0.0;
01658                   }
01659               }
01660           }
01661       }
01662 
01663     for (int cell=0;cell<numCells;cell++)
01664       {
01665         for (int field=0;field<numFields;field++)
01666           {
01667             for (int j=0;j<numBfy;j++)
01668               {
01669                 for (int i=0;i<numBfx;i++)
01670                   {
01671                     const int I=j*numBfx+i;
01672                     for (int m=0;m<numPtsx;m++)
01673                       {
01674                         const int Itmp = j * numBfy + m;
01675                         vals(cell,field,Itmp) += wtsx(m) * wtsy(j) * DPhix(i,m);
01676                       }
01677                   }
01678               }
01679           }
01680       }
01681 
01682     for (int cell=0;cell<numCells;cell++)
01683       {
01684         for (int field=0;field<numFields;field++)
01685           {
01686             for (int j=0;j<numBfy;j++)
01687               {
01688                 for (int i=0;i<numBfx;i++)
01689                   {
01690                     const int I=j*numBfx+i;
01691                     for (int n=0;n<numPtsy;n++)
01692                       {
01693                         const int Itmp = n * numBfy + i;
01694                         vals(cell,field,Itmp) += wtsx(i) * wtsy(n) * DPhiy(j,n);
01695                       }
01696                   }
01697               }
01698           }
01699       }
01700     
01701   }
01702 
01703   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01704            class ArrayTypeBasis, class ArrayTypeWeights>
01705   void TensorProductSpaceTools::momentsGrad3D( ArrayTypeOut &vals ,
01706                                                const ArrayTypeData &data ,
01707                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
01708                                                const Array<RCP<ArrayTypeBasis> > &basisDVals ,
01709                                                const Array<RCP<ArrayTypeWeights> > &wts )
01710   {
01711 
01712     const int numBfx = basisVals[0]->dimension(0);
01713     const int numBfy = basisVals[1]->dimension(0);
01714     const int numBfz = basisVals[2]->dimension(0);
01715     const int numPtsx = basisVals[0]->dimension(1);
01716     const int numPtsy = basisVals[1]->dimension(1);
01717     const int numPtsz = basisVals[2]->dimension(1);
01718     const int numCells = vals.dimension(0);
01719     const int numFields = vals.dimension(1);
01720     const ArrayTypeBasis &Phix = *basisVals[0];
01721     const ArrayTypeBasis &Phiy = *basisVals[1];
01722     const ArrayTypeBasis &Phiz = *basisVals[2];
01723     const ArrayTypeBasis &DPhix = *basisDVals[0];
01724     const ArrayTypeBasis &DPhiy = *basisDVals[1];
01725     const ArrayTypeBasis &DPhiz = *basisDVals[2];
01726     const ArrayTypeWeights &wtsx = *wts[0];
01727     const ArrayTypeWeights &wtsy = *wts[1];
01728     const ArrayTypeWeights &wtsz = *wts[2];
01729 
01730     FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy,3);
01731     FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz,3);
01732   
01733     // Xi phase
01734     for (int f=0;f<numFields;f++)
01735       {
01736         for (int c=0;c<numCells;c++)
01737           {
01738             for (int i=0;i<numBfx;i++)
01739               {
01740                 for (int n=0;n<numPtsz;n++)
01741                   {
01742                     for (int m=0;m<numPtsy;m++)
01743                       {
01744                         for (int l=0;l<numPtsx;l++)
01745                           {
01746                             const int dataIdx = n * numPtsy * numPtsx + m * numPtsx + l;
01747                             Xi(c,i,n,m,0) += wtsx(l) * DPhix(i,l,0) * data(c,f,dataIdx,0);
01748                             Xi(c,i,n,m,1) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,1);
01749                             Xi(c,i,n,m,2) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,2);
01750                           }
01751                       }
01752                   }
01753               }
01754           }
01755 
01756         // Theta phase
01757         for (int c=0;c<numCells;c++)
01758           {
01759             for (int j=0;j<numBfy;j++)
01760               {
01761                 for (int i=0;i<numBfx;i++)
01762                   {
01763                     for (int n=0;n<numPtsz;n++)
01764                       {
01765                         for (int m=0;m<numPtsy;m++)
01766                           {
01767                             Theta(c,j,i,n,0) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,0);
01768                             Theta(c,j,i,n,1) += wtsy(j) * DPhiy(j,m,0) * Xi(c,i,n,m,1);
01769                             Theta(c,j,i,n,2) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,2);
01770                           }
01771                       }
01772                   }
01773               }
01774           }
01775 
01776         // last phase
01777         for (int c=0;c<numCells;c++)
01778           {
01779             for (int k=0;k<numBfz;k++)
01780               {
01781                 for (int j=0;j<numBfy;j++)
01782                   {
01783                     for (int i=0;i<numBfx;i++)
01784                       {
01785                         const int momIdx = k * numBfx * numBfy + j * numBfx + i;
01786                         for (int n=0;n<numPtsz;n++)
01787                           {
01788                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,0) * Phiz(k,n);
01789                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,1) * Phiz(k,n);
01790                             vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,2) * DPhiz(k,n,0);
01791                           }
01792                       }
01793                   }
01794               }
01795           }
01796       }
01797 
01798 }
01799 
01800   template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
01801            class ArrayTypeBasis, class ArrayTypeWeights>
01802   void TensorProductSpaceTools::momentsGradCollocated3D( ArrayTypeOut &vals ,
01803                                                const ArrayTypeData &data ,
01804                                                const Array<RCP<ArrayTypeBasis> > &basisVals ,
01805                                                const Array<RCP<ArrayTypeBasis> > &basisDVals ,
01806                                                const Array<RCP<ArrayTypeWeights> > &wts )
01807   {
01808 
01809     const int numBfx = basisVals[0]->dimension(0);
01810     const int numBfy = basisVals[1]->dimension(0);
01811     const int numBfz = basisVals[2]->dimension(0);
01812     const int numPtsx = basisVals[0]->dimension(1);
01813     const int numPtsy = basisVals[1]->dimension(1);
01814     const int numPtsz = basisVals[2]->dimension(1);
01815     const int numCells = vals.dimension(0);
01816     const int numFields = vals.dimension(1);
01817     const ArrayTypeBasis &Phix = *basisVals[0];
01818     const ArrayTypeBasis &Phiy = *basisVals[1];
01819     const ArrayTypeBasis &Phiz = *basisVals[2];
01820     const ArrayTypeBasis &DPhix = *basisDVals[0];
01821     const ArrayTypeBasis &DPhiy = *basisDVals[1];
01822     const ArrayTypeBasis &DPhiz = *basisDVals[2];
01823     const ArrayTypeWeights &wtsx = *wts[0];
01824     const ArrayTypeWeights &wtsy = *wts[1];
01825     const ArrayTypeWeights &wtsz = *wts[2];
01826 
01827     for (int cell=0;cell<numCells;cell++)
01828       {
01829         for (int field=0;field<numFields;field++)
01830           {
01831             // x component of data versus x derivative of bases
01832             for (int k=0;k<numBfz;k++)
01833               {
01834                 for (int j=0;j<numBfy;j++)
01835                   {
01836                     for (int i=0;i<numBfx;i++)
01837                       {
01838                         const int I = numBfy * numBfx * k + numBfy * j + i;
01839                         for (int ell=0;ell<numPtsx;ell++)
01840                           {
01841                             const int Itmp = numBfy * numBfx * k + numBfy * j + ell;
01842                             vals(cell,field,I) += wtsx(ell) * wtsy(j) * wtsz(k) * DPhix( i , ell );
01843                           }
01844                       }
01845                   }
01846               }
01847             // y component of data versus y derivative of bases
01848             for (int k=0;k<numBfz;k++)
01849               {
01850                 for (int j=0;j<numBfy;j++)
01851                   {
01852                     for (int i=0;i<numBfx;i++)
01853                       {
01854                         const int I = numBfy * numBfx * k + numBfy * j + i;
01855                         for (int m=0;m<numPtsy;m++)
01856                           {
01857                             const int Itmp = numBfy * numBfx * k + numBfy * m + i;
01858                             vals(cell,field,I) += wtsx(i) * wtsy(m) * wtsz(k) * DPhiy( j , m );
01859                           }
01860                       }
01861                   }
01862               }
01863             // z component of data versus z derivative of bases
01864             for (int k=0;k<numBfz;k++)
01865               {
01866                 for (int j=0;j<numBfy;j++)
01867                   {
01868                     for (int i=0;i<numBfx;i++)
01869                       {
01870                         const int I = numBfy * numBfx * k + numBfy * j + i;
01871                         for (int n=0;n<numPtsz;n++)
01872                           {
01873                             const int Itmp = numBfy * numBfx * n + numBfy * j + i;
01874                             vals(cell,field,I) += wtsx(i) * wtsy(j) * wtsz(n) * DPhiz( k , n );
01875                           }
01876                       }
01877                   }
01878               }
01879           }
01880       }
01881   
01882 }
01883 
01884 
01885 
01886 } // end namespace Intrepid