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