Intrepid
http://trilinos.sandia.gov/packages/docs/r11.0/packages/intrepid/src/Discretization/FunctionSpaceTools/Intrepid_FunctionSpaceToolsDef.hpp
Go to the documentation of this file.
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