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