|
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 // 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 }
1.7.6.1