|
Intrepid
|
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 }
1.7.4