|
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 00052 template<class Scalar, class ArrayType> 00053 void FunctionSpaceToolsInPlace::HGRADtransformVALUE(ArrayType & inOutVals ) 00054 { 00055 return; 00056 } 00057 00058 template<class Scalar, class ArrayType> 00059 void FunctionSpaceToolsInPlace::HGRADtransformVALUEDual(ArrayType & inOutVals ) 00060 { 00061 return; 00062 } 00063 00064 template<class Scalar, class ArrayType, class ArrayTypeJac> 00065 void FunctionSpaceToolsInPlace::HGRADtransformGRAD(ArrayType & inOutVals, 00066 const ArrayTypeJac & jacobianInverse, 00067 const char transpose) 00068 { 00069 FieldContainer<Scalar> tmp(inOutVals.dimension(3)); 00070 00071 // test for transpose direction, one loop nest for each direction 00072 if (transpose == 'T') 00073 { 00074 for (int c=0;c<inOutVals.dimension(0);c++) 00075 { 00076 for (int f=0;f<inOutVals.dimension(1);f++) 00077 { 00078 for (int p=0;p<inOutVals.dimension(2);p++) 00079 { 00080 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00081 { 00082 tmp(d1) = 0.0; 00083 for (int d2=0;d2<inOutVals.dimension(3);d2++) 00084 { 00085 tmp(d1) += jacobianInverse(c,p,d2,d1) * inOutVals(c,f,p,d2); 00086 } 00087 } 00088 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00089 { 00090 inOutVals(c,f,p,d1) = tmp(d1); 00091 } 00092 } 00093 } 00094 } 00095 } 00096 else if (transpose == 'N') 00097 { 00098 for (int c=0;c<inOutVals.dimension(0);c++) 00099 { 00100 for (int f=0;f<inOutVals.dimension(1);f++) 00101 { 00102 for (int p=0;p<inOutVals.dimension(2);p++) 00103 { 00104 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00105 { 00106 tmp(d1) = 0.0; 00107 for (int d2=0;d2<inOutVals.dimension(3);d2++) 00108 { 00109 tmp(d1) += jacobianInverse(c,p,d1,d2) * inOutVals(c,f,p,d2); 00110 } 00111 } 00112 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00113 { 00114 inOutVals(c,f,p,d1) = tmp(d1); 00115 } 00116 } 00117 } 00118 } 00119 } 00120 else 00121 { 00122 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00123 "Intrepid:FunctionSpaceToolsInPlace::HGRADtransformGRAD::Unknown transpose type" ); 00124 } 00125 } 00126 00127 template<class Scalar, class ArrayType, class ArrayTypeJac> 00128 void FunctionSpaceToolsInPlace::HGRADtransformGRADDual(ArrayType & inOutVals, 00129 const ArrayTypeJac & jacobianInverse, 00130 const char transpose) 00131 { 00132 char t_loc; 00133 if (transpose == 'T') 00134 { 00135 t_loc = 'N'; 00136 } 00137 else if (transpose == 'N') 00138 { 00139 t_loc = 'T'; 00140 } 00141 else 00142 { 00143 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00144 "Intrepid:FunctionSpaceToolsInPlace::HGRADtransformGRADDual::Unknown transpose type" ); 00145 } 00146 FunctionSpaceToolsInPlace::HGRADtransformGRAD<Scalar,ArrayType,ArrayTypeJac>( inOutVals, 00147 jacobianInverse, 00148 t_loc); 00149 } 00150 00151 00152 template<class Scalar, class ArrayType, class ArrayTypeJac> 00153 void FunctionSpaceToolsInPlace::HCURLtransformVALUE(ArrayType & inOutVals, 00154 const ArrayTypeJac & jacobianInverse, 00155 const char transpose) 00156 { 00157 FunctionSpaceToolsInPlace::HGRADtransformGRAD<Scalar,ArrayType,ArrayTypeJac>( inOutVals , jacobianInverse, transpose ); 00158 } 00159 00160 template<class Scalar, class ArrayType, class ArrayTypeJac> 00161 void FunctionSpaceToolsInPlace::HCURLtransformVALUEDual(ArrayType & inOutVals, 00162 const ArrayTypeJac & jacobianInverse, 00163 const char transpose) 00164 { 00165 char t_loc; 00166 if (transpose == 'T') 00167 { 00168 t_loc = 'N'; 00169 } 00170 else if (transpose == 'N') 00171 { 00172 t_loc = 'T'; 00173 } 00174 else 00175 { 00176 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00177 "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformVALUEDual::Unknown transpose type" ); 00178 } 00179 FunctionSpaceToolsInPlace::HCURLtransformVALUEDual<Scalar,ArrayType,ArrayTypeJac>(inOutVals, 00180 jacobianInverse, 00181 t_loc); 00182 00183 } 00184 00185 template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet> 00186 void FunctionSpaceToolsInPlace::HCURLtransformCURL(ArrayType & inOutVals, 00187 const ArrayTypeJac & jacobian, 00188 const ArrayTypeDet & jacobianDet, 00189 const char transpose) 00190 { 00191 FieldContainer<Scalar> tmp(inOutVals.dimension(3)); 00192 00193 // test for transpose direction, one loop nest for each direction 00194 if (transpose == 'T') 00195 { 00196 for (int c=0;c<inOutVals.dimension(0);c++) 00197 { 00198 for (int f=0;f<inOutVals.dimension(1);f++) 00199 { 00200 for (int p=0;p<inOutVals.dimension(2);p++) 00201 { 00202 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00203 { 00204 tmp(d1) = 0.0; 00205 for (int d2=0;d2<inOutVals.dimension(3);d2++) 00206 { 00207 tmp(d1) += jacobian(c,p,d2,d1) * inOutVals(c,f,p,d2); 00208 } 00209 } 00210 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00211 { 00212 inOutVals(c,f,p,d1) = tmp(d1) / jacobianDet(c,p); 00213 } 00214 } 00215 } 00216 } 00217 } 00218 else if (transpose == 'N') 00219 { 00220 for (int c=0;c<inOutVals.dimension(0);c++) 00221 { 00222 for (int f=0;f<inOutVals.dimension(1);f++) 00223 { 00224 for (int p=0;p<inOutVals.dimension(2);p++) 00225 { 00226 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00227 { 00228 tmp(d1) = 0.0; 00229 for (int d2=0;d2<inOutVals.dimension(3);d2++) 00230 { 00231 tmp(d1) += jacobian(c,p,d1,d2) * inOutVals(c,f,p,d2); 00232 } 00233 } 00234 for (int d1=0;d1<inOutVals.dimension(3);d1++) 00235 { 00236 inOutVals(c,f,p,d1) = tmp(d1) / jacobianDet(c,p); 00237 } 00238 } 00239 } 00240 } 00241 } 00242 else 00243 { 00244 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00245 "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformCURL::Unknown transpose type" ); 00246 } 00247 00248 } 00249 00250 template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet> 00251 void FunctionSpaceToolsInPlace::HCURLtransformCURLDual(ArrayType & inOutVals, 00252 const ArrayTypeJac & jacobian, 00253 const ArrayTypeDet & jacobianDet, 00254 const char transpose) 00255 { 00256 char t_loc; 00257 if (transpose == 'T') 00258 { 00259 t_loc = 'N'; 00260 } 00261 else if (transpose == 'N') 00262 { 00263 t_loc = 'T'; 00264 } 00265 else 00266 { 00267 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00268 "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformCURLDual::Unknown transpose type" ); 00269 } 00270 FunctionSpaceToolsInPlace::HCURLtransformCURLDual<Scalar,ArrayType,ArrayTypeJac>(inOutVals, 00271 jacobian, 00272 jacobianDet, 00273 t_loc); 00274 } 00275 00276 template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet> 00277 void FunctionSpaceToolsInPlace::HDIVtransformVALUE(ArrayType & inOutVals, 00278 const ArrayTypeJac & jacobian, 00279 const ArrayTypeDet & jacobianDet, 00280 const char transpose) 00281 { 00282 FunctionSpaceToolsInPlace::HCURLtransformCURL<Scalar,ArrayType,ArrayTypeJac,ArrayTypeDet>( inOutVals,jacobian,jacobianDet,transpose); 00283 } 00284 00285 template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet> 00286 void FunctionSpaceToolsInPlace::HDIVtransformVALUEDual(ArrayType & inOutVals, 00287 const ArrayTypeJac & jacobian, 00288 const ArrayTypeDet & jacobianDet, 00289 const char transpose) 00290 { 00291 char t_loc; 00292 if (transpose == 'T') 00293 { 00294 t_loc = 'N'; 00295 } 00296 else if (transpose == 'N') 00297 { 00298 t_loc = 'T'; 00299 } 00300 else 00301 { 00302 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , 00303 "FunctionSpaceToolsInPlace::HDIVtransformVALUEDual: invalid transpose character"); 00304 } 00305 FunctionSpaceToolsInPlace::HDIVtransformVALUE<Scalar,ArrayType,ArrayTypeJac,ArrayTypeDet>( inOutVals, 00306 jacobian, 00307 jacobianDet , 00308 t_loc ); 00309 } 00310 00311 template<class Scalar, class ArrayType, class ArrayTypeDet> 00312 void FunctionSpaceToolsInPlace::HDIVtransformDIV(ArrayType & inOutVals, 00313 const ArrayTypeDet & jacobianDet) 00314 { 00315 for (int c=0;c<inOutVals.dimension(0);c++) 00316 { 00317 for (int f=0;f<inOutVals.dimension(1);f++) 00318 { 00319 for (int p=0;p<inOutVals.dimension(2);p++) 00320 { 00321 inOutVals(c,f,p) /= jacobianDet(c,p); 00322 } 00323 } 00324 } 00325 } 00326 00327 template<class Scalar, class ArrayType, class ArrayTypeDet> 00328 void FunctionSpaceToolsInPlace::HDIVtransformDIVDual(ArrayType & inOutVals, 00329 const ArrayTypeDet & jacobianDet) 00330 { 00331 FunctionSpaceToolsInPlace::HDIVtransformDIV<Scalar,ArrayType,ArrayTypeDet>( inOutVals , 00332 jacobianDet ); 00333 } 00334 00335 template<class Scalar, class ArrayType, class ArrayTypeDet> 00336 void FunctionSpaceToolsInPlace::HVOLtransformVALUE(ArrayType & inOutVals, 00337 const ArrayTypeDet & jacobianDet) 00338 { 00339 FunctionSpaceToolsInPlace::HDIVtransformDIV<Scalar,ArrayType,ArrayTypeDet>( inOutVals,jacobianDet); 00340 } 00341 00342 template<class Scalar, class ArrayType, class ArrayTypeDet> 00343 void FunctionSpaceToolsInPlace::HVOLtransformVALUEDual(ArrayType & inOutVals, 00344 const ArrayTypeDet & jacobianDet) 00345 { 00346 FunctionSpaceToolsInPlace::HVOLtransformVALUEDual( inOutVals ,jacobianDet ); 00347 } 00348 00349 00350 00351 template<class Scalar, class ArrayType, class ArrayTypeMeasure> 00352 void FunctionSpaceToolsInPlace::multiplyMeasure(ArrayType & inOutVals, 00353 const ArrayTypeMeasure & inMeasure) 00354 { 00355 if (inOutVals.rank() == 2) // inOutVals is (C,P) 00356 { 00357 for (int c=0;c<inOutVals.dimension(0);c++) 00358 { 00359 for (int p=0;p<inOutVals.dimension(0);p++) 00360 { 00361 inOutVals(c,p) *= inMeasure(c,p); 00362 } 00363 } 00364 } 00365 else if (inOutVals.rank() == 3) // inOutVals is (C,F,P) 00366 { 00367 for (int c=0;c<inOutVals.dimension(0);c++) 00368 { 00369 for (int f=0;f<inOutVals.dimension(1);f++) 00370 { 00371 for (int p=0;p<inOutVals.dimension(2);p++) 00372 { 00373 inOutVals(c,f,p) *= inMeasure(c,p); 00374 } 00375 } 00376 } 00377 } 00378 else if (inOutVals.rank() == 4) // inOutVals is (C,F,P) 00379 { 00380 for (int c=0;c<inOutVals.dimension(0);c++) 00381 { 00382 for (int f=0;f<inOutVals.dimension(1);f++) 00383 { 00384 for (int p=0;p<inOutVals.dimension(2);p++) 00385 { 00386 for (int d=0;d<inOutVals.dimension(3);d++) 00387 { 00388 inOutVals(c,f,p,d) *= inMeasure(c,p); 00389 } 00390 } 00391 } 00392 } 00393 } 00394 } 00395 00396 } // end namespace Intrepid
1.7.6.1