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