|
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 ArrayTypeOut, class ArrayTypeIn> 00053 void FunctionSpaceTools::HGRADtransformVALUE(ArrayTypeOut & outVals, 00054 const ArrayTypeIn & inVals) { 00055 00056 ArrayTools::cloneFields<Scalar>(outVals, inVals); 00057 00058 } 00059 00060 00061 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn> 00062 void FunctionSpaceTools::HGRADtransformGRAD(ArrayTypeOut & outVals, 00063 const ArrayTypeJac & jacobianInverse, 00064 const ArrayTypeIn & inVals, 00065 const char transpose) { 00066 00067 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobianInverse, inVals, transpose); 00068 00069 } 00070 00071 00072 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeIn> 00073 void FunctionSpaceTools::HCURLtransformVALUE(ArrayTypeOut & outVals, 00074 const ArrayTypeJac & jacobianInverse, 00075 const ArrayTypeIn & inVals, 00076 const char transpose) { 00077 00078 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobianInverse, inVals, transpose); 00079 00080 } 00081 00082 00083 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn> 00084 void FunctionSpaceTools::HCURLtransformCURL(ArrayTypeOut & outVals, 00085 const ArrayTypeJac & jacobian, 00086 const ArrayTypeDet & jacobianDet, 00087 const ArrayTypeIn & inVals, 00088 const char transpose) { 00089 00090 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobian, inVals, transpose); 00091 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, outVals, true); 00092 00093 } 00094 00095 00096 template<class Scalar, class ArrayTypeOut, class ArrayTypeJac, class ArrayTypeDet, class ArrayTypeIn> 00097 void FunctionSpaceTools::HDIVtransformVALUE(ArrayTypeOut & outVals, 00098 const ArrayTypeJac & jacobian, 00099 const ArrayTypeDet & jacobianDet, 00100 const ArrayTypeIn & inVals, 00101 const char transpose) { 00102 00103 ArrayTools::matvecProductDataField<Scalar>(outVals, jacobian, inVals, transpose); 00104 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, outVals, true); 00105 00106 } 00107 00108 00109 template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn> 00110 void FunctionSpaceTools::HDIVtransformDIV(ArrayTypeOut & outVals, 00111 const ArrayTypeDet & jacobianDet, 00112 const ArrayTypeIn & inVals) { 00113 00114 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true); 00115 00116 } 00117 00118 00119 template<class Scalar, class ArrayTypeOut, class ArrayTypeDet, class ArrayTypeIn> 00120 void FunctionSpaceTools::HVOLtransformVALUE(ArrayTypeOut & outVals, 00121 const ArrayTypeDet & jacobianDet, 00122 const ArrayTypeIn & inVals) { 00123 00124 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, jacobianDet, inVals, true); 00125 00126 } 00127 00128 00129 template<class Scalar, class ArrayOut, class ArrayInLeft, class ArrayInRight> 00130 void FunctionSpaceTools::integrate(ArrayOut & outputValues, 00131 const ArrayInLeft & leftValues, 00132 const ArrayInRight & rightValues, 00133 const ECompEngine compEngine, 00134 const bool sumInto) { 00135 int outRank = outputValues.rank(); 00136 00137 switch (outRank) { 00138 case 1: 00139 dataIntegral<Scalar>(outputValues, leftValues, rightValues, compEngine, sumInto); 00140 break; 00141 case 2: 00142 functionalIntegral<Scalar>(outputValues, leftValues, rightValues, compEngine, sumInto); 00143 break; 00144 case 3: 00145 operatorIntegral<Scalar>(outputValues, leftValues, rightValues, compEngine, sumInto); 00146 break; 00147 default: 00148 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 1) && (outRank != 2) && (outRank != 3)), std::invalid_argument, 00149 ">>> ERROR (FunctionSpaceTools::integrate): Output container must have rank 1, 2 or 3."); 00150 } 00151 00152 } // integrate 00153 00154 00155 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight> 00156 void FunctionSpaceTools::operatorIntegral(ArrayOutFields & outputFields, 00157 const ArrayInFieldsLeft & leftFields, 00158 const ArrayInFieldsRight & rightFields, 00159 const ECompEngine compEngine, 00160 const bool sumInto) { 00161 int lRank = leftFields.rank(); 00162 00163 switch (lRank) { 00164 case 3: 00165 ArrayTools::contractFieldFieldScalar<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto); 00166 break; 00167 case 4: 00168 ArrayTools::contractFieldFieldVector<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto); 00169 break; 00170 case 5: 00171 ArrayTools::contractFieldFieldTensor<Scalar>(outputFields, leftFields, rightFields, compEngine, sumInto); 00172 break; 00173 default: 00174 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 3) && (lRank != 4) && (lRank != 5)), std::invalid_argument, 00175 ">>> ERROR (FunctionSpaceTools::operatorIntegral): Left fields input container must have rank 3, 4 or 5."); 00176 } 00177 00178 } // operatorIntegral 00179 00180 00181 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00182 void FunctionSpaceTools::functionalIntegral(ArrayOutFields & outputFields, 00183 const ArrayInData & inputData, 00184 const ArrayInFields & inputFields, 00185 const ECompEngine compEngine, 00186 const bool sumInto) { 00187 int dRank = inputData.rank(); 00188 00189 switch (dRank) { 00190 case 2: 00191 ArrayTools::contractDataFieldScalar<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto); 00192 break; 00193 case 3: 00194 ArrayTools::contractDataFieldVector<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto); 00195 break; 00196 case 4: 00197 ArrayTools::contractDataFieldTensor<Scalar>(outputFields, inputData, inputFields, compEngine, sumInto); 00198 break; 00199 default: 00200 TEUCHOS_TEST_FOR_EXCEPTION( ((dRank != 2) && (dRank != 3) && (dRank != 4)), std::invalid_argument, 00201 ">>> ERROR (FunctionSpaceTools::functionalIntegral): Data input container must have rank 2, 3 or 4."); 00202 } 00203 00204 } // functionalIntegral 00205 00206 00207 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00208 void FunctionSpaceTools::dataIntegral(ArrayOutData & outputData, 00209 const ArrayInDataLeft & inputDataLeft, 00210 const ArrayInDataRight & inputDataRight, 00211 const ECompEngine compEngine, 00212 const bool sumInto) { 00213 int lRank = inputDataLeft.rank(); 00214 00215 switch (lRank) { 00216 case 2: 00217 ArrayTools::contractDataDataScalar<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto); 00218 break; 00219 case 3: 00220 ArrayTools::contractDataDataVector<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto); 00221 break; 00222 case 4: 00223 ArrayTools::contractDataDataTensor<Scalar>(outputData, inputDataLeft, inputDataRight, compEngine, sumInto); 00224 break; 00225 default: 00226 TEUCHOS_TEST_FOR_EXCEPTION( ((lRank != 2) && (lRank != 3) && (lRank != 4)), std::invalid_argument, 00227 ">>> ERROR (FunctionSpaceTools::dataIntegral): Left data input container must have rank 2, 3 or 4."); 00228 } 00229 00230 } // dataIntegral 00231 00232 00233 00234 template<class Scalar, class ArrayOut, class ArrayDet, class ArrayWeights> 00235 inline void FunctionSpaceTools::computeCellMeasure(ArrayOut & outVals, 00236 const ArrayDet & inDet, 00237 const ArrayWeights & inWeights) { 00238 00239 #ifdef HAVE_INTREPID_DEBUG 00240 TEUCHOS_TEST_FOR_EXCEPTION( (inDet.rank() != 2), std::invalid_argument, 00241 ">>> ERROR (FunctionSpaceTools::computeCellMeasure): Input determinants container must have rank 2."); 00242 #endif 00243 00244 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, inDet, inWeights); 00245 // must use absolute value of inDet, so flip sign where needed 00246 for (int cell=0; cell<outVals.dimension(0); cell++) { 00247 if (inDet(cell,0) < 0.0) { 00248 for (int point=0; point<outVals.dimension(1); point++) { 00249 outVals(cell, point) *= -1.0; 00250 } 00251 } 00252 } 00253 00254 } // computeCellMeasure 00255 00256 00257 00258 template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights> 00259 void FunctionSpaceTools::computeFaceMeasure(ArrayOut & outVals, 00260 const ArrayJac & inJac, 00261 const ArrayWeights & inWeights, 00262 const int whichFace, 00263 const shards::CellTopology & parentCell) { 00264 00265 #ifdef HAVE_INTREPID_DEBUG 00266 TEUCHOS_TEST_FOR_EXCEPTION( (inJac.rank() != 4), std::invalid_argument, 00267 ">>> ERROR (FunctionSpaceTools::computeFaceMeasure): Input Jacobian container must have rank 4."); 00268 #endif 00269 00270 // temporary storage for face normals 00271 FieldContainer<Scalar> faceNormals(inJac.dimension(0), inJac.dimension(1), inJac.dimension(2)); 00272 00273 // compute normals 00274 CellTools<Scalar>::getPhysicalFaceNormals(faceNormals, inJac, whichFace, parentCell); 00275 00276 // compute lenghts of normals 00277 RealSpaceTools<Scalar>::vectorNorm(outVals, faceNormals, NORM_TWO); 00278 00279 // multiply with weights 00280 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, outVals, inWeights); 00281 00282 } 00283 00284 00285 00286 template<class Scalar, class ArrayOut, class ArrayJac, class ArrayWeights> 00287 void FunctionSpaceTools::computeEdgeMeasure(ArrayOut & outVals, 00288 const ArrayJac & inJac, 00289 const ArrayWeights & inWeights, 00290 const int whichEdge, 00291 const shards::CellTopology & parentCell) { 00292 00293 #ifdef HAVE_INTREPID_DEBUG 00294 TEUCHOS_TEST_FOR_EXCEPTION( (inJac.rank() != 4), std::invalid_argument, 00295 ">>> ERROR (FunctionSpaceTools::computeEdgeMeasure): Input Jacobian container must have rank 4."); 00296 #endif 00297 00298 // temporary storage for edge tangents 00299 FieldContainer<Scalar> edgeTangents(inJac.dimension(0), inJac.dimension(1), inJac.dimension(2)); 00300 00301 // compute normals 00302 CellTools<Scalar>::getPhysicalEdgeTangents(edgeTangents, inJac, whichEdge, parentCell); 00303 00304 // compute lenghts of tangents 00305 RealSpaceTools<Scalar>::vectorNorm(outVals, edgeTangents, NORM_TWO); 00306 00307 // multiply with weights 00308 ArrayTools::scalarMultiplyDataData<Scalar>(outVals, outVals, inWeights); 00309 00310 } 00311 00312 00313 00314 template<class Scalar, class ArrayTypeOut, class ArrayTypeMeasure, class ArrayTypeIn> 00315 void FunctionSpaceTools::multiplyMeasure(ArrayTypeOut & outVals, 00316 const ArrayTypeMeasure & inMeasure, 00317 const ArrayTypeIn & inVals) { 00318 00319 ArrayTools::scalarMultiplyDataField<Scalar>(outVals, inMeasure, inVals); 00320 00321 } // multiplyMeasure 00322 00323 00324 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00325 void FunctionSpaceTools::scalarMultiplyDataField(ArrayOutFields & outputFields, 00326 ArrayInData & inputData, 00327 ArrayInFields & inputFields, 00328 const bool reciprocal) { 00329 00330 ArrayTools::scalarMultiplyDataField<Scalar>(outputFields, inputData, inputFields, reciprocal); 00331 00332 } // scalarMultiplyDataField 00333 00334 00335 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00336 void FunctionSpaceTools::scalarMultiplyDataData(ArrayOutData & outputData, 00337 ArrayInDataLeft & inputDataLeft, 00338 ArrayInDataRight & inputDataRight, 00339 const bool reciprocal) { 00340 00341 ArrayTools::scalarMultiplyDataData<Scalar>(outputData, inputDataLeft, inputDataRight, reciprocal); 00342 00343 } // scalarMultiplyDataData 00344 00345 00346 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00347 void FunctionSpaceTools::dotMultiplyDataField(ArrayOutFields & outputFields, 00348 const ArrayInData & inputData, 00349 const ArrayInFields & inputFields) { 00350 00351 ArrayTools::dotMultiplyDataField<Scalar>(outputFields, inputData, inputFields); 00352 00353 } // dotMultiplyDataField 00354 00355 00356 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00357 void FunctionSpaceTools::dotMultiplyDataData(ArrayOutData & outputData, 00358 const ArrayInDataLeft & inputDataLeft, 00359 const ArrayInDataRight & inputDataRight) { 00360 00361 ArrayTools::dotMultiplyDataData<Scalar>(outputData, inputDataLeft, inputDataRight); 00362 00363 } // dotMultiplyDataData 00364 00365 00366 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00367 void FunctionSpaceTools::vectorMultiplyDataField(ArrayOutFields & outputFields, 00368 const ArrayInData & inputData, 00369 const ArrayInFields & inputFields) { 00370 00371 int outRank = outputFields.rank(); 00372 00373 switch (outRank) { 00374 case 3: 00375 case 4: 00376 ArrayTools::crossProductDataField<Scalar>(outputFields, inputData, inputFields); 00377 break; 00378 case 5: 00379 ArrayTools::outerProductDataField<Scalar>(outputFields, inputData, inputFields); 00380 break; 00381 default: 00382 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 3) && (outRank != 4) && (outRank != 5)), std::invalid_argument, 00383 ">>> ERROR (FunctionSpaceTools::vectorMultiplyDataField): Output container must have rank 3, 4 or 5."); 00384 } 00385 00386 } // vectorMultiplyDataField 00387 00388 00389 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00390 void FunctionSpaceTools::vectorMultiplyDataData(ArrayOutData & outputData, 00391 const ArrayInDataLeft & inputDataLeft, 00392 const ArrayInDataRight & inputDataRight) { 00393 00394 int outRank = outputData.rank(); 00395 00396 switch (outRank) { 00397 case 2: 00398 case 3: 00399 ArrayTools::crossProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight); 00400 break; 00401 case 4: 00402 ArrayTools::outerProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight); 00403 break; 00404 default: 00405 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 2) && (outRank != 3) && (outRank != 4)), std::invalid_argument, 00406 ">>> ERROR (FunctionSpaceTools::vectorMultiplyDataData): Output container must have rank 2, 3 or 4."); 00407 } 00408 00409 } // vectorMultiplyDataData 00410 00411 00412 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00413 void FunctionSpaceTools::tensorMultiplyDataField(ArrayOutFields & outputFields, 00414 const ArrayInData & inputData, 00415 const ArrayInFields & inputFields, 00416 const char transpose) { 00417 00418 int outRank = outputFields.rank(); 00419 00420 switch (outRank) { 00421 case 4: 00422 ArrayTools::matvecProductDataField<Scalar>(outputFields, inputData, inputFields, transpose); 00423 break; 00424 case 5: 00425 ArrayTools::matmatProductDataField<Scalar>(outputFields, inputData, inputFields, transpose); 00426 break; 00427 default: 00428 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 4) && (outRank != 5)), std::invalid_argument, 00429 ">>> ERROR (FunctionSpaceTools::tensorMultiplyDataField): Output container must have rank 4 or 5."); 00430 } 00431 00432 } // tensorMultiplyDataField 00433 00434 00435 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00436 void FunctionSpaceTools::tensorMultiplyDataData(ArrayOutData & outputData, 00437 const ArrayInDataLeft & inputDataLeft, 00438 const ArrayInDataRight & inputDataRight, 00439 const char transpose) { 00440 00441 int outRank = outputData.rank(); 00442 00443 switch (outRank) { 00444 case 3: 00445 ArrayTools::matvecProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose); 00446 break; 00447 case 4: 00448 ArrayTools::matmatProductDataData<Scalar>(outputData, inputDataLeft, inputDataRight, transpose); 00449 break; 00450 default: 00451 TEUCHOS_TEST_FOR_EXCEPTION( ((outRank != 3) && (outRank != 4)), std::invalid_argument, 00452 ">>> ERROR (FunctionSpaceTools::tensorMultiplyDataData): Output container must have rank 3 or 4."); 00453 } 00454 00455 } // tensorMultiplyDataData 00456 00457 00458 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign> 00459 void FunctionSpaceTools::applyLeftFieldSigns(ArrayTypeInOut & inoutOperator, 00460 const ArrayTypeSign & fieldSigns) { 00461 #ifdef HAVE_INTREPID_DEBUG 00462 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.rank() != 3), std::invalid_argument, 00463 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Input operator container must have rank 3."); 00464 TEUCHOS_TEST_FOR_EXCEPTION( (fieldSigns.rank() != 2), std::invalid_argument, 00465 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Input field signs container must have rank 2."); 00466 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument, 00467 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): Zeroth dimensions (number of cells) of the operator and field signs containers must agree!"); 00468 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(1) != fieldSigns.dimension(1) ), std::invalid_argument, 00469 ">>> ERROR (FunctionSpaceTools::applyLeftFieldSigns): First dimensions (number of left fields) of the operator and field signs containers must agree!"); 00470 #endif 00471 00472 for (int cell=0; cell<inoutOperator.dimension(0); cell++) { 00473 for (int lbf=0; lbf<inoutOperator.dimension(1); lbf++) { 00474 for (int rbf=0; rbf<inoutOperator.dimension(2); rbf++) { 00475 inoutOperator(cell, lbf, rbf) *= fieldSigns(cell, lbf); 00476 } 00477 } 00478 } 00479 00480 } // applyLeftFieldSigns 00481 00482 00483 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign> 00484 void FunctionSpaceTools::applyRightFieldSigns(ArrayTypeInOut & inoutOperator, 00485 const ArrayTypeSign & fieldSigns) { 00486 #ifdef HAVE_INTREPID_DEBUG 00487 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.rank() != 3), std::invalid_argument, 00488 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Input operator container must have rank 3."); 00489 TEUCHOS_TEST_FOR_EXCEPTION( (fieldSigns.rank() != 2), std::invalid_argument, 00490 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Input field signs container must have rank 2."); 00491 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument, 00492 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Zeroth dimensions (number of cells) of the operator and field signs containers must agree!"); 00493 TEUCHOS_TEST_FOR_EXCEPTION( (inoutOperator.dimension(2) != fieldSigns.dimension(1) ), std::invalid_argument, 00494 ">>> ERROR (FunctionSpaceTools::applyRightFieldSigns): Second dimension of the operator container and first dimension of the field signs container (number of right fields) must agree!"); 00495 #endif 00496 00497 for (int cell=0; cell<inoutOperator.dimension(0); cell++) { 00498 for (int lbf=0; lbf<inoutOperator.dimension(1); lbf++) { 00499 for (int rbf=0; rbf<inoutOperator.dimension(2); rbf++) { 00500 inoutOperator(cell, lbf, rbf) *= fieldSigns(cell, rbf); 00501 } 00502 } 00503 } 00504 00505 } // applyRightFieldSigns 00506 00507 00508 template<class Scalar, class ArrayTypeInOut, class ArrayTypeSign> 00509 void FunctionSpaceTools::applyFieldSigns(ArrayTypeInOut & inoutFunction, 00510 const ArrayTypeSign & fieldSigns) { 00511 00512 #ifdef HAVE_INTREPID_DEBUG 00513 TEUCHOS_TEST_FOR_EXCEPTION( ((inoutFunction.rank() < 2) || (inoutFunction.rank() > 5)), std::invalid_argument, 00514 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Input function container must have rank 2, 3, 4, or 5."); 00515 TEUCHOS_TEST_FOR_EXCEPTION( (fieldSigns.rank() != 2), std::invalid_argument, 00516 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Input field signs container must have rank 2."); 00517 TEUCHOS_TEST_FOR_EXCEPTION( (inoutFunction.dimension(0) != fieldSigns.dimension(0) ), std::invalid_argument, 00518 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Zeroth dimensions (number of integration domains) of the function and field signs containers must agree!"); 00519 TEUCHOS_TEST_FOR_EXCEPTION( (inoutFunction.dimension(1) != fieldSigns.dimension(1) ), std::invalid_argument, 00520 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): First dimensions (number of fields) of the function and field signs containers must agree!"); 00521 #endif 00522 00523 int numCells = inoutFunction.dimension(0); 00524 int numFields = inoutFunction.dimension(1); 00525 int fRank = inoutFunction.rank(); 00526 00527 switch (fRank) { 00528 case 2: { 00529 for (int cell=0; cell<numCells; cell++) { 00530 for (int bf=0; bf<numFields; bf++) { 00531 inoutFunction(cell, bf) *= fieldSigns(cell, bf); 00532 } 00533 } 00534 } 00535 break; 00536 00537 case 3: { 00538 int numPoints = inoutFunction.dimension(2); 00539 for (int cell=0; cell<numCells; cell++) { 00540 for (int bf=0; bf<numFields; bf++) { 00541 for (int pt=0; pt<numPoints; pt++) { 00542 inoutFunction(cell, bf, pt) *= fieldSigns(cell, bf); 00543 } 00544 } 00545 } 00546 } 00547 break; 00548 00549 case 4: { 00550 int numPoints = inoutFunction.dimension(2); 00551 int spaceDim1 = inoutFunction.dimension(3); 00552 for (int cell=0; cell<numCells; cell++) { 00553 for (int bf=0; bf<numFields; bf++) { 00554 for (int pt=0; pt<numPoints; pt++) { 00555 for (int d1=0; d1<spaceDim1; d1++) { 00556 inoutFunction(cell, bf, pt, d1) *= fieldSigns(cell, bf); 00557 } 00558 } 00559 } 00560 } 00561 } 00562 break; 00563 00564 case 5: { 00565 int numPoints = inoutFunction.dimension(2); 00566 int spaceDim1 = inoutFunction.dimension(3); 00567 int spaceDim2 = inoutFunction.dimension(4); 00568 for (int cell=0; cell<numCells; cell++) { 00569 for (int bf=0; bf<numFields; bf++) { 00570 for (int pt=0; pt<numPoints; pt++) { 00571 for (int d1=0; d1<spaceDim1; d1++) { 00572 for (int d2=0; d2<spaceDim2; d2++) { 00573 inoutFunction(cell, bf, pt, d1, d2) *= fieldSigns(cell, bf); 00574 } 00575 } 00576 } 00577 } 00578 } 00579 } 00580 break; 00581 00582 default: 00583 TEUCHOS_TEST_FOR_EXCEPTION( !( (fRank == 2) || (fRank == 3) || (fRank == 4) || (fRank == 5)), std::invalid_argument, 00584 ">>> ERROR (FunctionSpaceTools::applyFieldSigns): Method defined only for rank-2, 3, 4, or 5 input function containers."); 00585 00586 } // end switch fRank 00587 00588 } // applyFieldSigns 00589 00590 00591 template<class Scalar, class ArrayOutPointVals, class ArrayInCoeffs, class ArrayInFields> 00592 void FunctionSpaceTools::evaluate(ArrayOutPointVals & outPointVals, 00593 const ArrayInCoeffs & inCoeffs, 00594 const ArrayInFields & inFields) { 00595 00596 #ifdef HAVE_INTREPID_DEBUG 00597 TEUCHOS_TEST_FOR_EXCEPTION( ((inFields.rank() < 3) || (inFields.rank() > 5)), std::invalid_argument, 00598 ">>> ERROR (FunctionSpaceTools::evaluate): Input fields container must have rank 3, 4, or 5."); 00599 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.rank() != 2), std::invalid_argument, 00600 ">>> ERROR (FunctionSpaceTools::evaluate): Input coefficient container must have rank 2."); 00601 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.rank() != inFields.rank()-1), std::invalid_argument, 00602 ">>> ERROR (FunctionSpaceTools::evaluate): Output values container must have rank one less than the rank of the input fields container."); 00603 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.dimension(0) != inFields.dimension(0) ), std::invalid_argument, 00604 ">>> ERROR (FunctionSpaceTools::evaluate): Zeroth dimensions (number of cells) of the coefficient and fields input containers must agree!"); 00605 TEUCHOS_TEST_FOR_EXCEPTION( (inCoeffs.dimension(1) != inFields.dimension(1) ), std::invalid_argument, 00606 ">>> ERROR (FunctionSpaceTools::evaluate): First dimensions (number of fields) of the coefficient and fields input containers must agree!"); 00607 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.dimension(0) != inFields.dimension(0) ), std::invalid_argument, 00608 ">>> ERROR (FunctionSpaceTools::evaluate): Zeroth dimensions (number of cells) of the input fields container and the output values container must agree!"); 00609 for (int i=1; i<outPointVals.rank(); i++) { 00610 std::string errmsg = ">>> ERROR (FunctionSpaceTools::evaluate): Dimensions "; 00611 errmsg += (char)(48+i); 00612 errmsg += " and "; 00613 errmsg += (char)(48+i+1); 00614 errmsg += " of the output values and input fields containers must agree!"; 00615 TEUCHOS_TEST_FOR_EXCEPTION( (outPointVals.dimension(i) != inFields.dimension(i+1)), std::invalid_argument, errmsg ); 00616 } 00617 #endif 00618 00619 int numCells = inFields.dimension(0); 00620 int numFields = inFields.dimension(1); 00621 int numPoints = inFields.dimension(2); 00622 int fRank = inFields.rank(); 00623 00624 switch (fRank) { 00625 case 3: { 00626 for (int cell=0; cell<numCells; cell++) { 00627 for (int pt=0; pt<numPoints; pt++) { 00628 for (int bf=0; bf<numFields; bf++) { 00629 outPointVals(cell, pt) += inCoeffs(cell, bf) * inFields(cell, bf, pt); 00630 } 00631 } 00632 } 00633 } 00634 break; 00635 00636 case 4: { 00637 int spaceDim1 = inFields.dimension(3); 00638 for (int cell=0; cell<numCells; cell++) { 00639 for (int pt=0; pt<numPoints; pt++) { 00640 for (int d1=0; d1<spaceDim1; d1++) { 00641 for (int bf=0; bf<numFields; bf++) { 00642 outPointVals(cell, pt, d1) += inCoeffs(cell, bf) * inFields(cell, bf, pt, d1); 00643 } 00644 } 00645 } 00646 } 00647 } 00648 break; 00649 00650 case 5: { 00651 int spaceDim1 = inFields.dimension(3); 00652 int spaceDim2 = inFields.dimension(4); 00653 for (int cell=0; cell<numCells; cell++) { 00654 for (int pt=0; pt<numPoints; pt++) { 00655 for (int d1=0; d1<spaceDim1; d1++) { 00656 for (int d2=0; d2<spaceDim2; d2++) { 00657 for (int bf=0; bf<numFields; bf++) { 00658 outPointVals(cell, pt, d1, d2) += inCoeffs(cell, bf) * inFields(cell, bf, pt, d1, d2); 00659 } 00660 } 00661 } 00662 } 00663 } 00664 } 00665 break; 00666 00667 default: 00668 TEUCHOS_TEST_FOR_EXCEPTION( !( (fRank == 3) || (fRank == 4) || (fRank == 5)), std::invalid_argument, 00669 ">>> ERROR (FunctionSpaceTools::evaluate): Method defined only for rank-3, 4, or 5 input fields containers."); 00670 00671 } // end switch fRank 00672 00673 } // evaluate 00674 00675 } // end namespace Intrepid
1.7.6.1