Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HDIV_TRI_In_FEM/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 
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Teuchos_oblackholestream.hpp"
00052 #include "Teuchos_RCP.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 #include "Intrepid_PointTools.hpp"
00055 #include "Intrepid_HDIV_TRI_In_FEM.hpp"
00056 #include "Shards_CellTopology.hpp"
00057 
00058 #include <iostream>
00059 using namespace Intrepid;
00060 
00065 int main(int argc, char *argv[]) {
00066 
00067   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00068   
00069   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00070   int iprint     = argc - 1;
00071   
00072   Teuchos::RCP<std::ostream> outStream;
00073   Teuchos::oblackholestream bhs; // outputs nothing
00074   
00075   if (iprint > 0)
00076     outStream = Teuchos::rcp(&std::cout, false);
00077   else
00078     outStream = Teuchos::rcp(&bhs, false);
00079   
00080   // Save the format state of the original std::cout.
00081   Teuchos::oblackholestream oldFormatState;
00082   oldFormatState.copyfmt(std::cout);
00083   
00084   *outStream \
00085     << "===============================================================================\n" \
00086     << "|                                                                             |\n" \
00087     << "|                           Unit Test HDIV_TRI_In_FEM                         |\n" \
00088     << "|                                                                             |\n" \
00089     << "|     1) Tests triangular Raviart-Thomas basis                                |\n" \
00090     << "|                                                                             |\n" \
00091     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00092     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00093     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00094     << "|                                                                             |\n" \
00095     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00096     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00097     << "|                                                                             |\n" \
00098     << "===============================================================================\n";
00099   
00100   int errorFlag  = 0;
00101 
00102   // test for basis values, compare against fiat
00103   try {
00104     const int deg = 2;
00105     Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >  myBasis( deg , POINTTYPE_EQUISPACED );
00106 
00107     // Get a lattice
00108     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00109     FieldContainer<double> lattice( np_lattice , 2 );
00110     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 );
00111     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00112                                                             myBasis.getBaseCellTopology() , 
00113                                                             deg , 
00114                                                             0 , 
00115                                                             POINTTYPE_EQUISPACED );    
00116 
00117     myBasis.getValues( myBasisValues , lattice , OPERATOR_VALUE );
00118 
00119     const double fiat_vals[] = {
00120       0.000000000000000e+00, -2.000000000000000e+00,
00121       2.500000000000000e-01, -5.000000000000000e-01,
00122       -1.000000000000000e+00, 1.000000000000000e+00,
00123       0.000000000000000e+00, -2.500000000000000e-01,
00124       -5.000000000000000e-01, 5.000000000000000e-01,
00125       0.000000000000000e+00, 0.000000000000000e+00,
00126       0.000000000000000e+00, 1.000000000000000e+00,
00127       2.500000000000000e-01, -5.000000000000000e-01,
00128       2.000000000000000e+00, -2.000000000000000e+00,
00129       0.000000000000000e+00, 5.000000000000000e-01,
00130       2.500000000000000e-01, -2.500000000000000e-01,
00131       0.000000000000000e+00, 0.000000000000000e+00,
00132       0.000000000000000e+00, 0.000000000000000e+00,
00133       2.500000000000000e-01, 0.000000000000000e+00,
00134       2.000000000000000e+00, 0.000000000000000e+00,
00135       0.000000000000000e+00, -5.000000000000000e-01,
00136       2.500000000000000e-01, 2.500000000000000e-01,
00137       0.000000000000000e+00, -1.000000000000000e+00,
00138       0.000000000000000e+00, 0.000000000000000e+00,
00139       -5.000000000000000e-01, 0.000000000000000e+00,
00140       -1.000000000000000e+00, 0.000000000000000e+00,
00141       0.000000000000000e+00, 2.500000000000000e-01,
00142       2.500000000000000e-01, 2.500000000000000e-01,
00143       0.000000000000000e+00, 2.000000000000000e+00,
00144       1.000000000000000e+00, 0.000000000000000e+00,
00145       5.000000000000000e-01, 0.000000000000000e+00,
00146       0.000000000000000e+00, 0.000000000000000e+00,
00147       -5.000000000000000e-01, 2.500000000000000e-01,
00148       -2.500000000000000e-01, 2.500000000000000e-01,
00149       -2.000000000000000e+00, 2.000000000000000e+00,
00150       -2.000000000000000e+00, 0.000000000000000e+00,
00151       -2.500000000000000e-01, 0.000000000000000e+00,
00152       0.000000000000000e+00, 0.000000000000000e+00,
00153       -5.000000000000000e-01, 2.500000000000000e-01,
00154       5.000000000000000e-01, -5.000000000000000e-01,
00155       1.000000000000000e+00, -1.000000000000000e+00,
00156       0.000000000000000e+00, 0.000000000000000e+00,
00157       1.500000000000000e+00, 0.000000000000000e+00,
00158       0.000000000000000e+00, 0.000000000000000e+00,
00159       0.000000000000000e+00, 7.500000000000000e-01,
00160       7.500000000000000e-01, -7.500000000000000e-01,
00161       0.000000000000000e+00, 0.000000000000000e+00,
00162       0.000000000000000e+00, 0.000000000000000e+00,
00163       7.500000000000000e-01, 0.000000000000000e+00,
00164       0.000000000000000e+00, 0.000000000000000e+00,
00165       0.000000000000000e+00, 1.500000000000000e+00,
00166       -7.500000000000000e-01, 7.500000000000000e-01,
00167       0.000000000000000e+00, 0.000000000000000e+00
00168     };
00169 
00170     int cur=0;
00171     for (int i=0;i<myBasisValues.dimension(0);i++) {
00172       for (int j=0;j<myBasisValues.dimension(1);j++) {
00173         for (int k=0;k<myBasisValues.dimension(2);k++) {
00174           if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) {
00175             errorFlag++;
00176             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00177             
00178             // Output the multi-index of the value where the error is:
00179             *outStream << " At multi-index { ";
00180             *outStream << i << " " << j << " " << k;
00181             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00182                       << " but correct value: " << fiat_vals[cur] << "\n";
00183             *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
00184           }
00185           cur++;
00186         }
00187       }
00188     }
00189   }
00190   catch (std::exception err) {
00191     *outStream << err.what() << "\n\n";
00192     errorFlag = -1000;
00193   }
00194   try {
00195     const int deg = 2;
00196     Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >  myBasis( deg , POINTTYPE_EQUISPACED );
00197 
00198     // Get a lattice
00199     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00200     FieldContainer<double> lattice( np_lattice , 2 );
00201     FieldContainer<double> myBasisDivs( myBasis.getCardinality() , np_lattice );
00202     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00203                                                             myBasis.getBaseCellTopology() , 
00204                                                             deg , 
00205                                                             0 , 
00206                                                             POINTTYPE_EQUISPACED );    
00207 
00208     myBasis.getValues( myBasisDivs , lattice , OPERATOR_DIV );
00209 
00210 
00211     const double fiat_divs[] = {
00212       7.000000000000000e+00,
00213       2.500000000000000e+00,
00214       -2.000000000000000e+00,
00215       2.500000000000000e+00,
00216       -2.000000000000000e+00,
00217       -2.000000000000000e+00,
00218       -2.000000000000000e+00,
00219       2.500000000000000e+00,
00220       7.000000000000000e+00,
00221       -2.000000000000000e+00,
00222       2.500000000000000e+00,
00223       -2.000000000000000e+00,
00224       -2.000000000000000e+00,
00225       2.500000000000000e+00,
00226       7.000000000000000e+00,
00227       -2.000000000000000e+00,
00228       2.500000000000000e+00,
00229       -2.000000000000000e+00,
00230       -2.000000000000000e+00,
00231       -2.000000000000000e+00,
00232       -2.000000000000000e+00,
00233       2.500000000000000e+00,
00234       2.500000000000000e+00,
00235       7.000000000000000e+00,
00236       -2.000000000000000e+00,
00237       -2.000000000000000e+00,
00238       -2.000000000000000e+00,
00239       2.500000000000000e+00,
00240       2.500000000000000e+00,
00241       7.000000000000000e+00,
00242       7.000000000000000e+00,
00243       2.500000000000000e+00,
00244       -2.000000000000000e+00,
00245       2.500000000000000e+00,
00246       -2.000000000000000e+00,
00247       -2.000000000000000e+00,
00248       9.000000000000000e+00,
00249       0.000000000000000e+00,
00250       -9.000000000000000e+00,
00251       4.500000000000000e+00,
00252       -4.500000000000000e+00,
00253       0.000000000000000e+00,
00254       9.000000000000000e+00,
00255       4.500000000000000e+00,
00256       0.000000000000000e+00,
00257       0.000000000000000e+00,
00258       -4.500000000000000e+00,
00259       -9.000000000000000e+00
00260     };
00261 
00262     int cur=0;
00263     for (int i=0;i<myBasisDivs.dimension(0);i++) {
00264       for (int j=0;j<myBasisDivs.dimension(1);j++) {
00265         if (std::abs( myBasisDivs(i,j) - fiat_divs[cur] ) > INTREPID_TOL ) {
00266           errorFlag++;
00267           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00268           
00269           // Output the multi-index of the value where the error is:
00270           *outStream << " At multi-index { ";
00271           *outStream << i << " " << j;
00272           *outStream << "}  computed value: " << myBasisDivs(i,j)
00273                     << " but correct value: " << fiat_divs[cur] << "\n";
00274           *outStream << "Difference: " << std::abs( myBasisDivs(i,j) - fiat_divs[cur] ) << "\n";
00275         }
00276         cur++;
00277       }
00278     }
00279   }
00280   catch (std::exception err) {
00281     *outStream << err.what() << "\n\n";
00282     errorFlag = -1000;
00283   }
00284 
00285 
00286   if (errorFlag != 0)
00287     std::cout << "End Result: TEST FAILED\n";
00288   else
00289     std::cout << "End Result: TEST PASSED\n";
00290   
00291   // reset format state of std::cout
00292   std::cout.copyfmt(oldFormatState);
00293   
00294   return errorFlag;
00295 }