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