Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/FunctionSpaceTools/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
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/src/Discretization/FunctionSpaceTools/Intrepid_FunctionSpaceToolsDef.hpp
