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