|
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 00048 #include "Intrepid_FieldContainer.hpp" 00049 #include "Intrepid_HDIV_WEDGE_I1_FEM.hpp" 00050 #include "Teuchos_oblackholestream.hpp" 00051 #include "Teuchos_RCP.hpp" 00052 #include "Teuchos_GlobalMPISession.hpp" 00053 00054 using namespace std; 00055 using namespace Intrepid; 00056 00057 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00058 { \ 00059 ++nException; \ 00060 try { \ 00061 S ; \ 00062 } \ 00063 catch (std::logic_error err) { \ 00064 ++throwCounter; \ 00065 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00066 *outStream << err.what() << '\n'; \ 00067 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00068 }; \ 00069 } 00070 00071 int main(int argc, char *argv[]) { 00072 00073 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00074 00075 // This little trick lets us print to std::cout only if 00076 // a (dummy) command-line argument is provided. 00077 int iprint = argc - 1; 00078 Teuchos::RCP<std::ostream> outStream; 00079 Teuchos::oblackholestream bhs; // outputs nothing 00080 if (iprint > 0) 00081 outStream = Teuchos::rcp(&std::cout, false); 00082 else 00083 outStream = Teuchos::rcp(&bhs, false); 00084 00085 // Save the format state of the original std::cout. 00086 Teuchos::oblackholestream oldFormatState; 00087 oldFormatState.copyfmt(std::cout); 00088 00089 *outStream \ 00090 << "===============================================================================\n" \ 00091 << "| |\n" \ 00092 << "| Unit Test (Basis_HDIV_WEDGE_I1_FEM) |\n" \ 00093 << "| |\n" \ 00094 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00095 << "| 2) Basis values for VALUE and DIV operators |\n" \ 00096 << "| |\n" \ 00097 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00098 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00099 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00100 << "| |\n" \ 00101 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00102 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00103 << "| |\n" \ 00104 << "===============================================================================\n"\ 00105 << "| TEST 1: Basis creation, exception testing |\n"\ 00106 << "===============================================================================\n"; 00107 00108 // Define basis and error flag 00109 Basis_HDIV_WEDGE_I1_FEM<double, FieldContainer<double> > wedgeBasis; 00110 int errorFlag = 0; 00111 00112 // Initialize throw counter for exception testing 00113 int nException = 0; 00114 int throwCounter = 0; 00115 00116 // Define array containing the 6 vertices of the reference WEDGE and 6 other points. 00117 FieldContainer<double> wedgeNodes(12, 3); 00118 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0; 00119 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0; 00120 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0; 00121 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0; 00122 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0; 00123 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0; 00124 00125 wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0; 00126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0; 00127 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.25; wedgeNodes(8,2) = 1.0; 00128 wedgeNodes(9,0) = 0.25; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75; 00129 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.5; wedgeNodes(10,2)= -0.25; 00130 wedgeNodes(11,0)= 0.5; wedgeNodes(11,1)= 0.5; wedgeNodes(11,2)= 0.0; 00131 00132 // Generic array for the output values; needs to be properly resized depending on the operator type 00133 FieldContainer<double> vals; 00134 00135 try{ 00136 // exception #1: GRAD cannot be applied to HDIV functions 00137 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00138 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 ); 00139 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException ); 00140 00141 // exception #2: CURL cannot be applied to HDIV functions 00142 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_CURL), throwCounter, nException ); 00143 00144 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00145 // getDofTag() to access invalid array elements thereby causing bounds check exception 00146 // exception #3 00147 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00148 // exception #4 00149 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00150 // exception #5 00151 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00152 // exception #6 00153 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(11), throwCounter, nException ); 00154 // exception #7 00155 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException ); 00156 00157 #ifdef HAVE_INTREPID_DEBUG 00158 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays 00159 // exception #8: input points array must be of rank-2 00160 FieldContainer<double> badPoints1(4, 5, 3); 00161 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00162 00163 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00164 FieldContainer<double> badPoints2(4, 2); 00165 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00166 00167 // exception #10 output values must be of rank-3 for OPERATOR_VALUE 00168 FieldContainer<double> badVals1(4, 3); 00169 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00170 00171 // exception #11 output values must be of rank-2 for OPERATOR_DIV 00172 FieldContainer<double> badVals2(4, 3, 1); 00173 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00174 00175 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00176 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0), 3); 00177 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00178 00179 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions) 00180 FieldContainer<double> badVals4(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0)); 00181 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00182 00183 // exception #14 incorrect 1st dimension of output array (must equal number of points) 00184 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1, 3); 00185 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00186 00187 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00188 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1); 00189 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00190 00191 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00192 FieldContainer<double> badVals7(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4); 00193 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals7, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00194 #endif 00195 00196 } 00197 catch (std::logic_error err) { 00198 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00199 *outStream << err.what() << '\n'; 00200 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00201 errorFlag = -1000; 00202 }; 00203 00204 // Check if number of thrown exceptions matches the one we expect 00205 if (throwCounter != nException) { 00206 errorFlag++; 00207 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00208 } 00209 00210 *outStream \ 00211 << "\n" 00212 << "===============================================================================\n"\ 00213 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00214 << "===============================================================================\n"; 00215 00216 try{ 00217 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags(); 00218 00219 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00220 for (unsigned i = 0; i < allTags.size(); i++) { 00221 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00222 00223 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00224 if( !( (myTag[0] == allTags[i][0]) && 00225 (myTag[1] == allTags[i][1]) && 00226 (myTag[2] == allTags[i][2]) && 00227 (myTag[3] == allTags[i][3]) ) ) { 00228 errorFlag++; 00229 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00230 *outStream << " getDofOrdinal( {" 00231 << allTags[i][0] << ", " 00232 << allTags[i][1] << ", " 00233 << allTags[i][2] << ", " 00234 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00235 *outStream << " getDofTag(" << bfOrd << ") = { " 00236 << myTag[0] << ", " 00237 << myTag[1] << ", " 00238 << myTag[2] << ", " 00239 << myTag[3] << "}\n"; 00240 } 00241 } 00242 00243 // Now do the same but loop over basis functions 00244 for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) { 00245 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00246 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00247 if( bfOrd != myBfOrd) { 00248 errorFlag++; 00249 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00250 *outStream << " getDofTag(" << bfOrd << ") = { " 00251 << myTag[0] << ", " 00252 << myTag[1] << ", " 00253 << myTag[2] << ", " 00254 << myTag[3] << "} but getDofOrdinal({" 00255 << myTag[0] << ", " 00256 << myTag[1] << ", " 00257 << myTag[2] << ", " 00258 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00259 } 00260 } 00261 } 00262 catch (std::logic_error err){ 00263 *outStream << err.what() << "\n\n"; 00264 errorFlag = -1000; 00265 }; 00266 00267 *outStream \ 00268 << "\n" 00269 << "===============================================================================\n"\ 00270 << "| TEST 3: correctness of basis function values |\n"\ 00271 << "===============================================================================\n"; 00272 00273 outStream -> precision(20); 00274 00275 // VALUE: Each row pair gives the 5x3 correct basis set values at an evaluation point 00276 double basisValues[] = { 00277 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -2.00000, 0, 0, 0, \ 00278 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, -2.00000, 0, \ 00279 0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, 0, 0, \ 00280 -2.00000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, \ 00281 0, 0, 0, 2.00000, 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, \ 00282 0, 0, 0, 0, 2.00000, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, \ 00283 0, 0, 0, 0, 0, 2.00000, 0.125000, -0.250000, 0, 0.125000, 0.250000, \ 00284 0, -0.375000, 0.250000, 0, 0, 0, -2.00000, 0, 0, 0, 0.250000, \ 00285 -0.375000, 0, 0.250000, 0.125000, 0, -0.250000, 0.125000, 0, 0, 0, \ 00286 -1.00000, 0, 0, 1.00000, 0.125000, -0.375000, 0, 0.125000, 0.125000, \ 00287 0, -0.375000, 0.125000, 0, 0, 0, 0, 0, 0, 2.00000, 0.125000, \ 00288 -0.500000, 0, 0.125000, 0, 0, -0.375000, 0, 0, 0, 0, -0.250000, 0, 0, \ 00289 1.75000, 0, -0.250000, 0, 0, 0.250000, 0, -0.500000, 0.250000, 0, 0, \ 00290 0, -1.25000, 0, 0, 0.750000, 0.250000, -0.250000, 0, 0.250000, \ 00291 0.250000, 0, -0.250000, 0.250000, 0, 0, 0, -1.00000, 0, 0, 1.00000}; 00292 00293 // DIV: each row pair gives the 5 correct values of the divergence of the 5 basis functions 00294 double basisDivs[] = { 00295 // 6 vertices 00296 1.0, 1.0, 1.0, 1.0, 1.0, 00297 1.0, 1.0, 1.0, 1.0, 1.0, 00298 1.0, 1.0, 1.0, 1.0, 1.0, 00299 1.0, 1.0, 1.0, 1.0, 1.0, 00300 1.0, 1.0, 1.0, 1.0, 1.0, 00301 1.0, 1.0, 1.0, 1.0, 1.0, 00302 // 6 other points 00303 1.0, 1.0, 1.0, 1.0, 1.0, 00304 1.0, 1.0, 1.0, 1.0, 1.0, 00305 1.0, 1.0, 1.0, 1.0, 1.0, 00306 1.0, 1.0, 1.0, 1.0, 1.0, 00307 1.0, 1.0, 1.0, 1.0, 1.0, 00308 1.0, 1.0, 1.0, 1.0, 1.0 00309 }; 00310 00311 try{ 00312 00313 // Dimensions for the output arrays: 00314 int numFields = wedgeBasis.getCardinality(); 00315 int numPoints = wedgeNodes.dimension(0); 00316 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension(); 00317 00318 // Generic array for values and curls that will be properly sized before each call 00319 FieldContainer<double> vals; 00320 00321 // Check VALUE of basis functions: resize vals to rank-3 container: 00322 vals.resize(numFields, numPoints, spaceDim); 00323 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE); 00324 for (int i = 0; i < numFields; i++) { 00325 for (int j = 0; j < numPoints; j++) { 00326 for (int k = 0; k < spaceDim; k++) { 00327 int l = k + i * spaceDim + j * spaceDim * numFields; 00328 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00329 errorFlag++; 00330 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00331 00332 // Output the multi-index of the value where the error is: 00333 *outStream << " At multi-index { "; 00334 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00335 *outStream << "} computed value: " << vals(i,j,k) 00336 << " but reference value: " << basisValues[l] << "\n"; 00337 } 00338 } 00339 } 00340 } 00341 00342 00343 // Check DIV of basis function: resize vals to rank-2 container 00344 vals.resize(numFields, numPoints); 00345 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV); 00346 for (int i = 0; i < numFields; i++) { 00347 for (int j = 0; j < numPoints; j++) { 00348 int l = i + j * numFields; 00349 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) { 00350 errorFlag++; 00351 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00352 00353 // Output the multi-index of the value where the error is: 00354 *outStream << " At multi-index { "; 00355 *outStream << i << " ";*outStream << j << " "; 00356 *outStream << "} computed divergence component: " << vals(i,j) 00357 << " but reference divergence component: " << basisDivs[l] << "\n"; 00358 } 00359 } 00360 } 00361 00362 } 00363 00364 // Catch unexpected errors 00365 catch (std::logic_error err) { 00366 *outStream << err.what() << "\n\n"; 00367 errorFlag = -1000; 00368 }; 00369 00370 if (errorFlag != 0) 00371 std::cout << "End Result: TEST FAILED\n"; 00372 else 00373 std::cout << "End Result: TEST PASSED\n"; 00374 00375 // reset format state of std::cout 00376 std::cout.copyfmt(oldFormatState); 00377 00378 return errorFlag; 00379 }
1.7.6.1