|
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_QUAD_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_QUAD_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_QUAD_I1_FEM<double, FieldContainer<double> > quadBasis; 00110 int errorFlag = 0; 00111 00112 // Initialize throw counter for exception testing 00113 int nException = 0; 00114 int throwCounter = 0; 00115 00116 // Array with the 4 vertices of the reference Quadrilateral, its center and 4 more points 00117 FieldContainer<double> quadNodes(9, 2); 00118 quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0; 00119 quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0; 00120 quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0; 00121 quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0; 00122 00123 quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0; 00124 quadNodes(5,0) = 0.0; quadNodes(5,1) = -0.5; 00125 quadNodes(6,0) = 0.0; quadNodes(6,1) = 0.5; 00126 quadNodes(7,0) = -0.5; quadNodes(7,1) = 0.0; 00127 quadNodes(8,0) = 0.5; quadNodes(8,1) = 0.0; 00128 00129 // Generic array for the output values; needs to be properly resized depending on the operator type 00130 FieldContainer<double> vals; 00131 00132 try{ 00133 00134 int spaceDim = quadBasis.getBaseCellTopology().getDimension(); 00135 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(quadBasis.getCardinality(), quadNodes.dimension(0), spaceDim ); 00139 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException ); 00140 00141 // exception #2: CURL cannot be applied to HDIV functions 00142 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, 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( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00148 // exception #4 00149 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00150 // exception #5 00151 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00152 // exception #6 00153 INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException ); 00154 // exception #7 00155 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException ); 00156 00157 #ifdef HAVE_INTREPID_DEBUG 00158 // Exceptions 8- 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( quadBasis.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, spaceDim + 1); 00165 INTREPID_TEST_COMMAND( quadBasis.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, 5); 00169 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00170 00171 // exception #11 output values must be of rank-2 for OPERATOR_DIV 00172 FieldContainer<double> badVals2(4, 5, spaceDim); 00173 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00174 00175 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00176 FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0), spaceDim); 00177 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00178 00179 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions) 00180 FieldContainer<double> badVals4(quadBasis.getCardinality() + 1, quadNodes.dimension(0)); 00181 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00182 00183 // exception #14 incorrect 1st dimension of output array (must equal number of points) 00184 FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, spaceDim); 00185 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00186 00187 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00188 FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0) + 1); 00189 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00190 00191 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00192 FieldContainer<double> badVals7(quadBasis.getCardinality(), quadNodes.dimension(0), spaceDim + 1); 00193 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals7, quadNodes, 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 = quadBasis.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 = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00222 00223 std::vector<int> myTag = quadBasis.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 < quadBasis.getCardinality(); bfOrd++) { 00245 std::vector<int> myTag = quadBasis.getDofTag(bfOrd); 00246 int myBfOrd = quadBasis.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 6x3 correct basis set values at an evaluation point: (P,F,D) layout 00276 double basisValues[] = { 00277 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, -0.500000, 0.500000, 0, 0, \ 00278 0, 0, 0, 0, 0, 0.500000, 0, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, \ 00279 0.500000, -0.500000, 0, 0, -0.250000, 0.250000, 0, 0, 0.250000, \ 00280 -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.125000, -0.250000, 0, \ 00281 0, -0.125000, 0.250000, 0, 0, 0.375000, -0.250000, 0, 0, -0.250000, \ 00282 0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.375000, 0, 0, \ 00283 0.250000, -0.125000, 0 }; 00284 00285 // DIV: each row gives the 6 correct values of the divergence of the 6 basis functions: (P,F) layout 00286 double basisDivs[] = { 00287 0.25, 0.25, 0.25, 0.25, 00288 0.25, 0.25, 0.25, 0.25, 00289 0.25, 0.25, 0.25, 0.25, 00290 0.25, 0.25, 0.25, 0.25, 00291 0.25, 0.25, 0.25, 0.25, 00292 0.25, 0.25, 0.25, 0.25, 00293 0.25, 0.25, 0.25, 0.25, 00294 0.25, 0.25, 0.25, 0.25, 00295 0.25, 0.25, 0.25, 0.25, 00296 }; 00297 00298 try{ 00299 00300 // Dimensions for the output arrays: 00301 int numPoints = quadNodes.dimension(0); 00302 int numFields = quadBasis.getCardinality(); 00303 int spaceDim = quadBasis.getBaseCellTopology().getDimension(); 00304 00305 // Generic array for values and curls that will be properly sized before each call 00306 FieldContainer<double> vals; 00307 00308 // Check VALUE of basis functions: resize vals to rank-3 container: 00309 vals.resize(numFields, numPoints, spaceDim); 00310 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE); 00311 for (int i = 0; i < numFields; i++) { 00312 for (int j = 0; j < numPoints; j++) { 00313 for (int k = 0; k < spaceDim; k++) { 00314 00315 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k 00316 int l = k + i * spaceDim + j * spaceDim * numFields; 00317 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00318 errorFlag++; 00319 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00320 00321 // Output the multi-index of the value where the error is: 00322 *outStream << " At multi-index { "; 00323 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00324 *outStream << "} computed value: " << vals(i,j,k) 00325 << " but reference value: " << basisValues[l] << "\n"; 00326 } 00327 } 00328 } 00329 } 00330 00331 // Check DIV of basis function: resize vals to rank-2 container 00332 vals.resize(numFields, numPoints); 00333 quadBasis.getValues(vals, quadNodes, OPERATOR_DIV); 00334 for (int i = 0; i < numFields; i++) { 00335 for (int j = 0; j < numPoints; j++) { 00336 int l = i + j * numFields; 00337 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) { 00338 errorFlag++; 00339 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00340 00341 // Output the multi-index of the value where the error is: 00342 *outStream << " At multi-index { "; 00343 *outStream << i << " ";*outStream << j << " "; 00344 *outStream << "} computed divergence component: " << vals(i,j) 00345 << " but reference divergence component: " << basisDivs[l] << "\n"; 00346 } 00347 } 00348 } 00349 00350 } 00351 00352 // Catch unexpected errors 00353 catch (std::logic_error err) { 00354 *outStream << err.what() << "\n\n"; 00355 errorFlag = -1000; 00356 }; 00357 00358 if (errorFlag != 0) 00359 std::cout << "End Result: TEST FAILED\n"; 00360 else 00361 std::cout << "End Result: TEST PASSED\n"; 00362 00363 // reset format state of std::cout 00364 std::cout.copyfmt(oldFormatState); 00365 00366 return errorFlag; 00367 }
1.7.6.1