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