Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/OrthogonalBases/test_01.cpp
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 
00044 
00051 #include "Intrepid_FieldContainer.hpp"
00052 #include "Teuchos_oblackholestream.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_GlobalMPISession.hpp"
00055 #include "Intrepid_OrthogonalBases.hpp"
00056 #include <Intrepid_CubatureDirectTriDefault.hpp>
00057 #include <iostream>
00058 using namespace Intrepid;
00059 
00064 int main(int argc, char *argv[]) {
00065 
00066   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00067   
00068   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00069   int iprint     = argc - 1;
00070   
00071   Teuchos::RCP<std::ostream> outStream;
00072   Teuchos::oblackholestream bhs; // outputs nothing
00073   
00074   if (iprint > 0)
00075     outStream = Teuchos::rcp(&std::cout, false);
00076   else
00077     outStream = Teuchos::rcp(&bhs, false);
00078   
00079   // Save the format state of the original std::cout.
00080   Teuchos::oblackholestream oldFormatState;
00081   oldFormatState.copyfmt(std::cout);
00082   
00083   *outStream \
00084     << "===============================================================================\n" \
00085     << "|                                                                             |\n" \
00086     << "|                           Unit Test OrthogonalBases                         |\n" \
00087     << "|                                                                             |\n" \
00088     << "|     1) Tests orthogonality of triangular orthogonal basis (Dubiner)         |\n" \
00089     << "|                                                                             |\n" \
00090     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00091     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00092     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00093     << "|                                                                             |\n" \
00094     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00095     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00096     << "|                                                                             |\n" \
00097     << "===============================================================================\n";
00098   
00099   int errorFlag  = 0;
00100   
00101   // First, get a reference quadrature rule
00102 
00103   CubatureDirectTriDefault<double,FieldContainer<double> > myCub(20);
00104   FieldContainer<double> cubPts( myCub.getNumPoints() , 2 );
00105   FieldContainer<double> cubWts( myCub.getNumPoints() );
00106 
00107   myCub.getCubature( cubPts , cubWts );
00108   
00109   // Tabulate the basis functions at the cubature points
00110   const int deg = 10;
00111   const int polydim = (deg+1)*(deg+2)/2;
00112   FieldContainer<double> basisAtCubPts( polydim , myCub.getNumPoints() );
00113   OrthogonalBases::tabulateTriangle<double,FieldContainer<double>,FieldContainer<double> >( cubPts , deg , basisAtCubPts );
00114 
00115   // Now let's compute the mass matrix
00116   for (int i=0;i<polydim;i++) {
00117     for (int j=0;j<polydim;j++) {
00118       double cur = 0;
00119       for (int k=0;k<myCub.getNumPoints();k++) {
00120         cur += cubWts(k) * basisAtCubPts( i , k ) * basisAtCubPts( j , k );
00121       }
00122       if (i != j && fabs( cur ) > INTREPID_TOL) {
00123         errorFlag++;
00124       }
00125       else if (i == j && fabs( cur ) < INTREPID_TOL ) {
00126         errorFlag++;
00127       }
00128 
00129     }
00130   }
00131 
00132   
00133   if (errorFlag != 0)
00134     std::cout << "End Result: TEST FAILED\n";
00135   else
00136     std::cout << "End Result: TEST PASSED\n";
00137   
00138   // reset format state of std::cout
00139   std::cout.copyfmt(oldFormatState);
00140   
00141   return errorFlag;
00142 }