|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // 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
1.7.6.1