|
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_HDIV_TRI_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 HDIV_TRI_In_FEM |\n" \ 00075 << "| |\n" \ 00076 << "| 1) Tests triangular Raviart-Thomas 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 // test for basis values, compare against fiat 00090 try { 00091 const int deg = 2; 00092 Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED ); 00093 00094 // Get a lattice 00095 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 ); 00096 FieldContainer<double> lattice( np_lattice , 2 ); 00097 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 ); 00098 PointTools::getLattice<double,FieldContainer<double> >( lattice , 00099 myBasis.getBaseCellTopology() , 00100 deg , 00101 0 , 00102 POINTTYPE_EQUISPACED ); 00103 00104 myBasis.getValues( myBasisValues , lattice , OPERATOR_VALUE ); 00105 00106 const double fiat_vals[] = { 00107 0.000000000000000e+00, -2.000000000000000e+00, 00108 2.500000000000000e-01, -5.000000000000000e-01, 00109 -1.000000000000000e+00, 1.000000000000000e+00, 00110 0.000000000000000e+00, -2.500000000000000e-01, 00111 -5.000000000000000e-01, 5.000000000000000e-01, 00112 0.000000000000000e+00, 0.000000000000000e+00, 00113 0.000000000000000e+00, 1.000000000000000e+00, 00114 2.500000000000000e-01, -5.000000000000000e-01, 00115 2.000000000000000e+00, -2.000000000000000e+00, 00116 0.000000000000000e+00, 5.000000000000000e-01, 00117 2.500000000000000e-01, -2.500000000000000e-01, 00118 0.000000000000000e+00, 0.000000000000000e+00, 00119 0.000000000000000e+00, 0.000000000000000e+00, 00120 2.500000000000000e-01, 0.000000000000000e+00, 00121 2.000000000000000e+00, 0.000000000000000e+00, 00122 0.000000000000000e+00, -5.000000000000000e-01, 00123 2.500000000000000e-01, 2.500000000000000e-01, 00124 0.000000000000000e+00, -1.000000000000000e+00, 00125 0.000000000000000e+00, 0.000000000000000e+00, 00126 -5.000000000000000e-01, 0.000000000000000e+00, 00127 -1.000000000000000e+00, 0.000000000000000e+00, 00128 0.000000000000000e+00, 2.500000000000000e-01, 00129 2.500000000000000e-01, 2.500000000000000e-01, 00130 0.000000000000000e+00, 2.000000000000000e+00, 00131 1.000000000000000e+00, 0.000000000000000e+00, 00132 5.000000000000000e-01, 0.000000000000000e+00, 00133 0.000000000000000e+00, 0.000000000000000e+00, 00134 -5.000000000000000e-01, 2.500000000000000e-01, 00135 -2.500000000000000e-01, 2.500000000000000e-01, 00136 -2.000000000000000e+00, 2.000000000000000e+00, 00137 -2.000000000000000e+00, 0.000000000000000e+00, 00138 -2.500000000000000e-01, 0.000000000000000e+00, 00139 0.000000000000000e+00, 0.000000000000000e+00, 00140 -5.000000000000000e-01, 2.500000000000000e-01, 00141 5.000000000000000e-01, -5.000000000000000e-01, 00142 1.000000000000000e+00, -1.000000000000000e+00, 00143 0.000000000000000e+00, 0.000000000000000e+00, 00144 1.500000000000000e+00, 0.000000000000000e+00, 00145 0.000000000000000e+00, 0.000000000000000e+00, 00146 0.000000000000000e+00, 7.500000000000000e-01, 00147 7.500000000000000e-01, -7.500000000000000e-01, 00148 0.000000000000000e+00, 0.000000000000000e+00, 00149 0.000000000000000e+00, 0.000000000000000e+00, 00150 7.500000000000000e-01, 0.000000000000000e+00, 00151 0.000000000000000e+00, 0.000000000000000e+00, 00152 0.000000000000000e+00, 1.500000000000000e+00, 00153 -7.500000000000000e-01, 7.500000000000000e-01, 00154 0.000000000000000e+00, 0.000000000000000e+00 00155 }; 00156 00157 int cur=0; 00158 for (int i=0;i<myBasisValues.dimension(0);i++) { 00159 for (int j=0;j<myBasisValues.dimension(1);j++) { 00160 for (int k=0;k<myBasisValues.dimension(2);k++) { 00161 if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) { 00162 errorFlag++; 00163 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00164 00165 // Output the multi-index of the value where the error is: 00166 *outStream << " At multi-index { "; 00167 *outStream << i << " " << j << " " << k; 00168 *outStream << "} computed value: " << myBasisValues(i,j,k) 00169 << " but correct value: " << fiat_vals[cur] << "\n"; 00170 *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n"; 00171 } 00172 cur++; 00173 } 00174 } 00175 } 00176 } 00177 catch (std::exception err) { 00178 *outStream << err.what() << "\n\n"; 00179 errorFlag = -1000; 00180 } 00181 try { 00182 const int deg = 2; 00183 Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED ); 00184 00185 // Get a lattice 00186 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 ); 00187 FieldContainer<double> lattice( np_lattice , 2 ); 00188 FieldContainer<double> myBasisDivs( myBasis.getCardinality() , np_lattice ); 00189 PointTools::getLattice<double,FieldContainer<double> >( lattice , 00190 myBasis.getBaseCellTopology() , 00191 deg , 00192 0 , 00193 POINTTYPE_EQUISPACED ); 00194 00195 myBasis.getValues( myBasisDivs , lattice , OPERATOR_DIV ); 00196 00197 00198 const double fiat_divs[] = { 00199 7.000000000000000e+00, 00200 2.500000000000000e+00, 00201 -2.000000000000000e+00, 00202 2.500000000000000e+00, 00203 -2.000000000000000e+00, 00204 -2.000000000000000e+00, 00205 -2.000000000000000e+00, 00206 2.500000000000000e+00, 00207 7.000000000000000e+00, 00208 -2.000000000000000e+00, 00209 2.500000000000000e+00, 00210 -2.000000000000000e+00, 00211 -2.000000000000000e+00, 00212 2.500000000000000e+00, 00213 7.000000000000000e+00, 00214 -2.000000000000000e+00, 00215 2.500000000000000e+00, 00216 -2.000000000000000e+00, 00217 -2.000000000000000e+00, 00218 -2.000000000000000e+00, 00219 -2.000000000000000e+00, 00220 2.500000000000000e+00, 00221 2.500000000000000e+00, 00222 7.000000000000000e+00, 00223 -2.000000000000000e+00, 00224 -2.000000000000000e+00, 00225 -2.000000000000000e+00, 00226 2.500000000000000e+00, 00227 2.500000000000000e+00, 00228 7.000000000000000e+00, 00229 7.000000000000000e+00, 00230 2.500000000000000e+00, 00231 -2.000000000000000e+00, 00232 2.500000000000000e+00, 00233 -2.000000000000000e+00, 00234 -2.000000000000000e+00, 00235 9.000000000000000e+00, 00236 0.000000000000000e+00, 00237 -9.000000000000000e+00, 00238 4.500000000000000e+00, 00239 -4.500000000000000e+00, 00240 0.000000000000000e+00, 00241 9.000000000000000e+00, 00242 4.500000000000000e+00, 00243 0.000000000000000e+00, 00244 0.000000000000000e+00, 00245 -4.500000000000000e+00, 00246 -9.000000000000000e+00 00247 }; 00248 00249 int cur=0; 00250 for (int i=0;i<myBasisDivs.dimension(0);i++) { 00251 for (int j=0;j<myBasisDivs.dimension(1);j++) { 00252 if (std::abs( myBasisDivs(i,j) - fiat_divs[cur] ) > INTREPID_TOL ) { 00253 errorFlag++; 00254 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00255 00256 // Output the multi-index of the value where the error is: 00257 *outStream << " At multi-index { "; 00258 *outStream << i << " " << j; 00259 *outStream << "} computed value: " << myBasisDivs(i,j) 00260 << " but correct value: " << fiat_divs[cur] << "\n"; 00261 *outStream << "Difference: " << std::abs( myBasisDivs(i,j) - fiat_divs[cur] ) << "\n"; 00262 } 00263 cur++; 00264 } 00265 } 00266 } 00267 catch (std::exception err) { 00268 *outStream << err.what() << "\n\n"; 00269 errorFlag = -1000; 00270 } 00271 00272 00273 if (errorFlag != 0) 00274 std::cout << "End Result: TEST FAILED\n"; 00275 else 00276 std::cout << "End Result: TEST PASSED\n"; 00277 00278 // reset format state of std::cout 00279 std::cout.copyfmt(oldFormatState); 00280 00281 return errorFlag; 00282 }
1.7.4