Intrepid
http://trilinos.sandia.gov/packages/docs/r10.6/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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or
00025 //                    Denis Ridzal (dridzal@sandia.gov) or
00026 //                    Robert Kirby (robert.c.kirby@ttu.edu)
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00031 
00037 #include "Intrepid_FieldContainer.hpp"
00038 #include "Teuchos_oblackholestream.hpp"
00039 #include "Teuchos_RCP.hpp"
00040 #include "Teuchos_GlobalMPISession.hpp"
00041 #include "Intrepid_PointTools.hpp"
00042 #include "Intrepid_HCURL_TET_In_FEM.hpp"
00043 #include "Shards_CellTopology.hpp"
00044 
00045 #include <iostream>
00046 using namespace Intrepid;
00047 
00052 int main(int argc, char *argv[]) {
00053 
00054   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00055   
00056   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00057   int iprint     = argc - 1;
00058   
00059   Teuchos::RCP<std::ostream> outStream;
00060   Teuchos::oblackholestream bhs; // outputs nothing
00061   
00062   if (iprint > 0)
00063     outStream = Teuchos::rcp(&std::cout, false);
00064   else
00065     outStream = Teuchos::rcp(&bhs, false);
00066   
00067   // Save the format state of the original std::cout.
00068   Teuchos::oblackholestream oldFormatState;
00069   oldFormatState.copyfmt(std::cout);
00070   
00071   *outStream \
00072     << "===============================================================================\n" \
00073     << "|                                                                             |\n" \
00074     << "|                           Unit Test HCURL_TET_In_FEM                        |\n" \
00075     << "|                                                                             |\n" \
00076     << "|     1) Tests tetrahedral Nedelec basis                                      |\n" \
00077     << "|                                                                             |\n" \
00078     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00079     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00080     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00081     << "|                                                                             |\n" \
00082     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00083     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00084     << "|                                                                             |\n" \
00085     << "===============================================================================\n";
00086   
00087   int errorFlag  = 0;
00088 
00089   // code-code comparison with FIAT 
00090   try {
00091     const int deg = 1;
00092     Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00093 
00094     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00095     FieldContainer<double> lattice( np_lattice , 3 );
00096     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
00097     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00098                                                             myBasis.getBaseCellTopology() , 
00099                                                             deg , 
00100                                                             0 , 
00101                                                             POINTTYPE_EQUISPACED );    
00102 
00103     myBasis.getValues( myBasisValues , lattice , OPERATOR_VALUE );
00104 
00105     const double fiat_vals[] = {
00106 1.000000000000001e+00, -2.498001805406602e-16, -1.665334536937735e-16,
00107 9.999999999999998e-01, 1.000000000000000e+00, 1.000000000000000e+00,
00108 5.828670879282072e-16, 1.110223024625157e-16, 2.498001805406602e-16,
00109 7.771561172376096e-16, 8.326672684688674e-17, 1.110223024625157e-16,
00110 2.081668171172169e-16, -2.914335439641036e-16, 1.280865063236792e-16,
00111 -3.191891195797325e-16, 1.000000000000000e+00, -4.293998586504916e-17,
00112 -9.999999999999994e-01, 2.081668171172169e-16, 2.400576428367544e-16,
00113 2.220446049250313e-16, -5.551115123125783e-17, 1.084013877651281e-16,
00114 3.469446951953614e-16, -1.000000000000000e+00, 1.387778780781446e-16,
00115 -1.804112415015879e-16, 1.942890293094024e-16, -1.387778780781446e-16,
00116 -9.999999999999993e-01, -9.999999999999996e-01, -9.999999999999998e-01,
00117 5.551115123125783e-17, -2.220446049250313e-16, -8.326672684688674e-17,
00118 -2.220446049250313e-16, -5.551115123125783e-17, 9.999999999999999e-01,
00119 1.665334536937735e-16, 1.110223024625157e-16, -6.383782391594650e-16,
00120 1.110223024625157e-16, 1.110223024625157e-16, -1.110223024625157e-16,
00121 9.999999999999990e-01, 9.999999999999994e-01, 9.999999999999996e-01,
00122 1.387778780781446e-16, -2.496931404305374e-17, -1.665334536937735e-16,
00123 -2.498001805406602e-16, -2.149987498083074e-16, 1.000000000000000e+00,
00124 8.326672684688674e-17, -3.769887250591415e-17, 8.326672684688674e-17,
00125 -9.999999999999994e-01, 1.556977698723022e-16, 2.220446049250313e-16,
00126 -9.422703950001342e-18, 1.665334536937735e-16, -2.359223927328458e-16,
00127 -9.422703950001268e-18, -8.326672684688674e-17, 1.387778780781446e-17,
00128 -7.525083148581445e-17, 2.775557561562891e-17, 1.000000000000000e+00,
00129 2.789513560273035e-16, -9.999999999999998e-01, -5.551115123125783e-17
00130     };
00131 
00132     int cur=0;
00133     for (int i=0;i<myBasisValues.dimension(0);i++) {
00134       for (int j=0;j<myBasisValues.dimension(1);j++) {
00135         for (int k=0;k<myBasisValues.dimension(2);k++) {
00136           if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > 10.0*INTREPID_TOL ) {
00137             errorFlag++;
00138             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00139             
00140             // Output the multi-index of the value where the error is:
00141             *outStream << " At multi-index { ";
00142             *outStream << i << " " << j << " " << k;
00143             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00144                        << " but correct value: " << fiat_vals[cur] << "\n";
00145             *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
00146           }
00147           cur++;
00148         }
00149       }
00150     }
00151   }
00152   catch (std::exception err) {
00153     *outStream << err.what() << "\n\n";
00154     errorFlag = -1000;
00155   }
00156 
00157   try {
00158     const int deg = 1;
00159     Basis_HCURL_TET_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00160 
00161     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00162     FieldContainer<double> lattice( np_lattice , 3 );
00163     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 3 );
00164     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00165                                                             myBasis.getBaseCellTopology() , 
00166                                                             deg , 
00167                                                             0 , 
00168                                                             POINTTYPE_EQUISPACED );    
00169 
00170     myBasis.getValues( myBasisValues , lattice , OPERATOR_CURL );
00171 
00172     const double fiat_curls[] = {
00173       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00174       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00175       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00176       -5.551115123125783e-16, -2.000000000000000e+00, 2.000000000000000e+00,
00177       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00178       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00179       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00180       -4.440892098500626e-16, -2.775557561562891e-16, 2.000000000000000e+00,
00181       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00182       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00183       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00184       -2.000000000000000e+00, 5.551115123125783e-17, 2.000000000000000e+00,
00185       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00186       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00187       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00188       -2.000000000000000e+00, 2.000000000000000e+00, 9.861075762086680e-17,
00189       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00190       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00191       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00192       -2.775557561562891e-17, -2.000000000000000e+00, 4.287451790760826e-16,
00193       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00194       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00195       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16,
00196       2.000000000000000e+00, -2.185751579730777e-16, 1.526556658859590e-16
00197     };
00198 
00199     int cur=0;
00200     for (int i=0;i<myBasisValues.dimension(0);i++) {
00201       for (int j=0;j<myBasisValues.dimension(1);j++) {
00202         for (int k=0;k<myBasisValues.dimension(2);k++) {
00203           if (std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) > 10.0*INTREPID_TOL ) {
00204             errorFlag++;
00205             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00206             
00207             // Output the multi-index of the value where the error is:
00208             *outStream << " At multi-index { ";
00209             *outStream << i << " " << j << " " << k;
00210             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00211                        << " but correct value: " << fiat_curls[cur] << "\n";
00212             *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_curls[cur] ) << "\n";
00213           }
00214           cur++;
00215         }
00216       }
00217     }
00218   }
00219   catch (std::exception err) {
00220     *outStream << err.what() << "\n\n";
00221     errorFlag = -1000;
00222   }  
00223 
00224 
00225   if (errorFlag != 0)
00226     std::cout << "End Result: TEST FAILED\n";
00227   else
00228     std::cout << "End Result: TEST PASSED\n";
00229   
00230   // reset format state of std::cout
00231   std::cout.copyfmt(oldFormatState);
00232   
00233   return errorFlag;
00234 }