Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HCURL_TET_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_HCURL_TET_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 HCURL_TET_In_FEM                        |\n" \
00088     << "|                                                                             |\n" \
00089     << "|     1) Tests tetrahedral Nedelec 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   // code-code comparison with FIAT 
00103   try {
00104     const int deg = 1;
00105     Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00106 
00107     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00108     FieldContainer<double> lattice( np_lattice , 3 );
00109     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
00110     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00111                                                             myBasis.getBaseCellTopology() , 
00112                                                             deg , 
00113                                                             0 , 
00114                                                             POINTTYPE_EQUISPACED );    
00115 
00116     myBasis.getValues( myBasisValues , lattice , OPERATOR_VALUE );
00117 
00118     const double fiat_vals[] = {
00119 1.000000000000001e+00, -2.498001805406602e-16, -1.665334536937735e-16,
00120 9.999999999999998e-01, 1.000000000000000e+00, 1.000000000000000e+00,
00121 5.828670879282072e-16, 1.110223024625157e-16, 2.498001805406602e-16,
00122 7.771561172376096e-16, 8.326672684688674e-17, 1.110223024625157e-16,
00123 2.081668171172169e-16, -2.914335439641036e-16, 1.280865063236792e-16,
00124 -3.191891195797325e-16, 1.000000000000000e+00, -4.293998586504916e-17,
00125 -9.999999999999994e-01, 2.081668171172169e-16, 2.400576428367544e-16,
00126 2.220446049250313e-16, -5.551115123125783e-17, 1.084013877651281e-16,
00127 3.469446951953614e-16, -1.000000000000000e+00, 1.387778780781446e-16,
00128 -1.804112415015879e-16, 1.942890293094024e-16, -1.387778780781446e-16,
00129 -9.999999999999993e-01, -9.999999999999996e-01, -9.999999999999998e-01,
00130 5.551115123125783e-17, -2.220446049250313e-16, -8.326672684688674e-17,
00131 -2.220446049250313e-16, -5.551115123125783e-17, 9.999999999999999e-01,
00132 1.665334536937735e-16, 1.110223024625157e-16, -6.383782391594650e-16,
00133 1.110223024625157e-16, 1.110223024625157e-16, -1.110223024625157e-16,
00134 9.999999999999990e-01, 9.999999999999994e-01, 9.999999999999996e-01,
00135 1.387778780781446e-16, -2.496931404305374e-17, -1.665334536937735e-16,
00136 -2.498001805406602e-16, -2.149987498083074e-16, 1.000000000000000e+00,
00137 8.326672684688674e-17, -3.769887250591415e-17, 8.326672684688674e-17,
00138 -9.999999999999994e-01, 1.556977698723022e-16, 2.220446049250313e-16,
00139 -9.422703950001342e-18, 1.665334536937735e-16, -2.359223927328458e-16,
00140 -9.422703950001268e-18, -8.326672684688674e-17, 1.387778780781446e-17,
00141 -7.525083148581445e-17, 2.775557561562891e-17, 1.000000000000000e+00,
00142 2.789513560273035e-16, -9.999999999999998e-01, -5.551115123125783e-17
00143     };
00144 
00145     int cur=0;
00146     for (int i=0;i<myBasisValues.dimension(0);i++) {
00147       for (int j=0;j<myBasisValues.dimension(1);j++) {
00148         for (int k=0;k<myBasisValues.dimension(2);k++) {
00149           if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > 10.0*INTREPID_TOL ) {
00150             errorFlag++;
00151             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00152             
00153             // Output the multi-index of the value where the error is:
00154             *outStream << " At multi-index { ";
00155             *outStream << i << " " << j << " " << k;
00156             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00157                        << " but correct value: " << fiat_vals[cur] << "\n";
00158             *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
00159           }
00160           cur++;
00161         }
00162       }
00163     }
00164   }
00165   catch (std::exception err) {
00166     *outStream << err.what() << "\n\n";
00167     errorFlag = -1000;
00168   }
00169 
00170   try {
00171     const int deg = 1;
00172     Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00173 
00174     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00175     FieldContainer<double> lattice( np_lattice , 3 );
00176     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
00177     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00178                                                             myBasis.getBaseCellTopology() , 
00179                                                             deg , 
00180                                                             0 , 
00181                                                             POINTTYPE_EQUISPACED );    
00182 
00183     myBasis.getValues( myBasisValues , lattice , OPERATOR_CURL );
00184 
00185     const double fiat_curls[] = {
00186       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00187       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00188       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00189       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00190       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00191       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00192       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00193       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00194       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00195       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00196       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00197       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00198       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00199       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00200       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00201       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00202       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00203       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00204       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00205       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00206       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00207       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00208       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00209       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16
00210     };
00211 
00212     int cur=0;
00213     for (int i=0;i<myBasisValues.dimension(0);i++) {
00214       for (int j=0;j<myBasisValues.dimension(1);j++) {
00215         for (int k=0;k<myBasisValues.dimension(2);k++) {
00216           if (std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) > 10.0*INTREPID_TOL ) {
00217             errorFlag++;
00218             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00219             
00220             // Output the multi-index of the value where the error is:
00221             *outStream << " At multi-index { ";
00222             *outStream << i << " " << j << " " << k;
00223             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00224                        << " but correct value: " << fiat_curls[cur] << "\n";
00225             *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) << "\n";
00226           }
00227           cur++;
00228         }
00229       }
00230     }
00231   }
00232   catch (std::exception err) {
00233     *outStream << err.what() << "\n\n";
00234     errorFlag = -1000;
00235   }  
00236 
00237 
00238   if (errorFlag != 0)
00239     std::cout << "End Result: TEST FAILED\n";
00240   else
00241     std::cout << "End Result: TEST PASSED\n";
00242   
00243   // reset format state of std::cout
00244   std::cout.copyfmt(oldFormatState);
00245   
00246   return errorFlag;
00247 }