|
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 "Intrepid_HCURL_TRI_In_FEM.hpp" 00044 #include "Shards_CellTopology.hpp" 00045 00046 #include <iostream> 00047 using namespace Intrepid; 00048 00053 int main(int argc, char *argv[]) { 00054 00055 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00056 00057 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided. 00058 int iprint = argc - 1; 00059 00060 Teuchos::RCP<std::ostream> outStream; 00061 Teuchos::oblackholestream bhs; // outputs nothing 00062 00063 if (iprint > 0) 00064 outStream = Teuchos::rcp(&std::cout, false); 00065 else 00066 outStream = Teuchos::rcp(&bhs, false); 00067 00068 // Save the format state of the original std::cout. 00069 Teuchos::oblackholestream oldFormatState; 00070 oldFormatState.copyfmt(std::cout); 00071 00072 *outStream \ 00073 << "===============================================================================\n" \ 00074 << "| |\n" \ 00075 << "| Unit Test HCURL_TRI_In_FEM |\n" \ 00076 << "| |\n" \ 00077 << "| 1) Tests triangular Nedelec-Thomas basis |\n" \ 00078 << "| |\n" \ 00079 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00080 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \ 00081 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \ 00082 << "| |\n" \ 00083 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00084 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00085 << "| |\n" \ 00086 << "===============================================================================\n"; 00087 00088 int errorFlag = 0; 00089 00090 // test for basis values, compare against H(div) basis, for they should be related by rotation 00091 // in the lowest order case 00092 try { 00093 const int deg = 1; 00094 Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> > myBasisDIV( deg , POINTTYPE_EQUISPACED ); 00095 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasisCURL( deg , POINTTYPE_EQUISPACED ); 00096 00097 // Get a lattice 00098 const int np_lattice = PointTools::getLatticeSize( myBasisDIV.getBaseCellTopology() , deg , 0 ); 00099 FieldContainer<double> lattice( np_lattice , 2 ); 00100 00101 FieldContainer<double> myBasisValuesDIV( myBasisDIV.getCardinality() , np_lattice , 2 ); 00102 FieldContainer<double> myBasisValuesCURL( myBasisDIV.getCardinality() , np_lattice , 2 ); 00103 PointTools::getLattice<double,FieldContainer<double> >( lattice , 00104 myBasisDIV.getBaseCellTopology() , 00105 deg , 00106 0 , 00107 POINTTYPE_EQUISPACED ); 00108 00109 myBasisDIV.getValues( myBasisValuesDIV , lattice , OPERATOR_VALUE ); 00110 myBasisCURL.getValues( myBasisValuesCURL, lattice , OPERATOR_VALUE ); 00111 00112 for (int i=0;i<myBasisValuesDIV.dimension(0);i++) { 00113 for (int j=0;j<myBasisValuesDIV.dimension(1);j++) { 00114 if (std::abs( myBasisValuesDIV(i,j,1) + myBasisValuesCURL(i,j,0) ) > INTREPID_TOL ) { 00115 errorFlag++; 00116 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00117 // Output the multi-index of the value where the error is: 00118 *outStream << " At multi-index { "; 00119 *outStream << i << " " << j << " and component 0"; 00120 *outStream << "} computed value: " << myBasisValuesCURL(i,j,0) 00121 << " but correct value: " << -myBasisValuesDIV(i,j,1) << "\n"; 00122 *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,0) + myBasisValuesDIV(i,j,1) ) << "\n"; 00123 } 00124 if (std::abs( myBasisValuesDIV(i,j,0) - myBasisValuesCURL(i,j,1) ) > INTREPID_TOL ) { 00125 errorFlag++; 00126 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00127 // Output the multi-index of the value where the error is: 00128 *outStream << " At multi-index { "; 00129 *outStream << i << " " << j << " and component 1"; 00130 *outStream << "} computed value: " << myBasisValuesCURL(i,j,1) 00131 << " but correct value: " << myBasisValuesDIV(i,j,0) << "\n"; 00132 *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,1) - myBasisValuesDIV(i,j,1) ) << "\n"; 00133 } 00134 } 00135 } 00136 } 00137 catch (std::exception err) { 00138 *outStream << err.what() << "\n\n"; 00139 errorFlag = -1000; 00140 } 00141 00142 // next test: code-code comparison with FIAT 00143 try { 00144 const int deg = 2; 00145 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED ); 00146 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 ); 00147 FieldContainer<double> lattice( np_lattice , 2 ); 00148 PointTools::getLattice<double,FieldContainer<double> >( lattice , 00149 myBasis.getBaseCellTopology() , 00150 deg , 00151 0 , 00152 POINTTYPE_EQUISPACED ); 00153 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 ); 00154 00155 00156 myBasis.getValues( myBasisValues, lattice , OPERATOR_VALUE ); 00157 00158 const double fiat_vals[] = { 00159 2.000000000000000e+00, 0.000000000000000e+00, 00160 5.000000000000000e-01, 2.500000000000000e-01, 00161 -1.000000000000000e+00, -1.000000000000000e+00, 00162 2.500000000000000e-01, 0.000000000000000e+00, 00163 -5.000000000000000e-01, -5.000000000000000e-01, 00164 0.000000000000000e+00, 0.000000000000000e+00, 00165 -1.000000000000000e+00, 0.000000000000000e+00, 00166 5.000000000000000e-01, 2.500000000000000e-01, 00167 2.000000000000000e+00, 2.000000000000000e+00, 00168 -5.000000000000000e-01, 0.000000000000000e+00, 00169 2.500000000000000e-01, 2.500000000000000e-01, 00170 0.000000000000000e+00, 0.000000000000000e+00, 00171 0.000000000000000e+00, 0.000000000000000e+00, 00172 0.000000000000000e+00, 2.500000000000000e-01, 00173 0.000000000000000e+00, 2.000000000000000e+00, 00174 5.000000000000000e-01, 0.000000000000000e+00, 00175 -2.500000000000000e-01, 2.500000000000000e-01, 00176 1.000000000000000e+00, 0.000000000000000e+00, 00177 0.000000000000000e+00, 0.000000000000000e+00, 00178 0.000000000000000e+00, -5.000000000000000e-01, 00179 0.000000000000000e+00, -1.000000000000000e+00, 00180 -2.500000000000000e-01, 0.000000000000000e+00, 00181 -2.500000000000000e-01, 2.500000000000000e-01, 00182 -2.000000000000000e+00, 0.000000000000000e+00, 00183 0.000000000000000e+00, 1.000000000000000e+00, 00184 0.000000000000000e+00, 5.000000000000000e-01, 00185 0.000000000000000e+00, 0.000000000000000e+00, 00186 -2.500000000000000e-01, -5.000000000000000e-01, 00187 -2.500000000000000e-01, -2.500000000000000e-01, 00188 -2.000000000000000e+00, -2.000000000000000e+00, 00189 0.000000000000000e+00, -2.000000000000000e+00, 00190 0.000000000000000e+00, -2.500000000000000e-01, 00191 0.000000000000000e+00, 0.000000000000000e+00, 00192 -2.500000000000000e-01, -5.000000000000000e-01, 00193 5.000000000000000e-01, 5.000000000000000e-01, 00194 1.000000000000000e+00, 1.000000000000000e+00, 00195 0.000000000000000e+00, 0.000000000000000e+00, 00196 0.000000000000000e+00, -7.500000000000000e-01, 00197 0.000000000000000e+00, 0.000000000000000e+00, 00198 1.500000000000000e+00, 0.000000000000000e+00, 00199 7.500000000000000e-01, 7.500000000000000e-01, 00200 0.000000000000000e+00, 0.000000000000000e+00, 00201 0.000000000000000e+00, 0.000000000000000e+00, 00202 0.000000000000000e+00, 1.500000000000000e+00, 00203 0.000000000000000e+00, 0.000000000000000e+00, 00204 -7.500000000000000e-01, 0.000000000000000e+00, 00205 7.500000000000000e-01, 7.500000000000000e-01, 00206 0.000000000000000e+00, 0.000000000000000e+00 00207 }; 00208 00209 int cur=0; 00210 for (int i=0;i<myBasisValues.dimension(0);i++) { 00211 for (int j=0;j<myBasisValues.dimension(1);j++) { 00212 for (int k=0;k<myBasisValues.dimension(2);k++) { 00213 if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) { 00214 errorFlag++; 00215 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00216 00217 // Output the multi-index of the value where the error is: 00218 *outStream << " At multi-index { "; 00219 *outStream << i << " " << j << " " << k; 00220 *outStream << "} computed value: " << myBasisValues(i,j,k) 00221 << " but correct value: " << fiat_vals[cur] << "\n"; 00222 *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n"; 00223 } 00224 cur++; 00225 } 00226 } 00227 } 00228 } 00229 catch (std::exception err) { 00230 *outStream << err.what() << "\n\n"; 00231 errorFlag = -1000; 00232 } 00233 00234 try { 00235 const int deg = 2; 00236 Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED ); 00237 const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 ); 00238 FieldContainer<double> lattice( np_lattice , 2 ); 00239 PointTools::getLattice<double,FieldContainer<double> >( lattice , 00240 myBasis.getBaseCellTopology() , 00241 deg , 00242 0 , 00243 POINTTYPE_EQUISPACED ); 00244 FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice ); 00245 00246 00247 myBasis.getValues( myBasisValues, lattice , OPERATOR_CURL ); 00248 00249 const double fiat_curls[] = { 00250 -7.000000000000000e+00, 00251 -2.500000000000000e+00, 00252 2.000000000000000e+00, 00253 -2.500000000000000e+00, 00254 2.000000000000000e+00, 00255 2.000000000000000e+00, 00256 2.000000000000000e+00, 00257 -2.500000000000000e+00, 00258 -7.000000000000000e+00, 00259 2.000000000000000e+00, 00260 -2.500000000000000e+00, 00261 2.000000000000000e+00, 00262 2.000000000000000e+00, 00263 -2.500000000000000e+00, 00264 -7.000000000000000e+00, 00265 2.000000000000000e+00, 00266 -2.500000000000000e+00, 00267 2.000000000000000e+00, 00268 2.000000000000000e+00, 00269 2.000000000000000e+00, 00270 2.000000000000000e+00, 00271 -2.500000000000000e+00, 00272 -2.500000000000000e+00, 00273 -7.000000000000000e+00, 00274 2.000000000000000e+00, 00275 2.000000000000000e+00, 00276 2.000000000000000e+00, 00277 -2.500000000000000e+00, 00278 -2.500000000000000e+00, 00279 -7.000000000000000e+00, 00280 -7.000000000000000e+00, 00281 -2.500000000000000e+00, 00282 2.000000000000000e+00, 00283 -2.500000000000000e+00, 00284 2.000000000000000e+00, 00285 2.000000000000000e+00, 00286 9.000000000000000e+00, 00287 4.500000000000000e+00, 00288 0.000000000000000e+00, 00289 0.000000000000000e+00, 00290 -4.500000000000000e+00, 00291 -9.000000000000000e+00, 00292 -9.000000000000000e+00, 00293 0.000000000000000e+00, 00294 9.000000000000000e+00, 00295 -4.500000000000000e+00, 00296 4.500000000000000e+00, 00297 0.000000000000000e+00 00298 }; 00299 00300 int cur=0; 00301 for (int i=0;i<myBasisValues.dimension(0);i++) { 00302 for (int j=0;j<myBasisValues.dimension(1);j++) { 00303 if (std::abs( myBasisValues(i,j) - fiat_curls[cur] ) > INTREPID_TOL ) { 00304 errorFlag++; 00305 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00306 00307 // Output the multi-index of the value where the error is: 00308 *outStream << " At multi-index { "; 00309 *outStream << i << " " << j; 00310 *outStream << "} computed value: " << myBasisValues(i,j) 00311 << " but correct value: " << fiat_curls[cur] << "\n"; 00312 *outStream << "Difference: " << std::abs( myBasisValues(i,j) - fiat_curls[cur] ) << "\n"; 00313 } 00314 cur++; 00315 } 00316 } 00317 } 00318 catch (std::exception err) { 00319 *outStream << err.what() << "\n\n"; 00320 errorFlag = -1000; 00321 } 00322 00323 if (errorFlag != 0) 00324 std::cout << "End Result: TEST FAILED\n"; 00325 else 00326 std::cout << "End Result: TEST PASSED\n"; 00327 00328 // reset format state of std::cout 00329 std::cout.copyfmt(oldFormatState); 00330 00331 return errorFlag; 00332 }
1.7.4