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