|
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_TET_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_TET_I1_FEM) |\n" \ 00080 << "| |\n" \ 00081 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00082 << "| 3) Exception tests on input arguments and invalid operators |\n" \ 00083 << "| 2) Basis values for VALUE and DIV operators |\n" \ 00084 << "| |\n" \ 00085 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00086 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00087 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00088 << "| |\n" \ 00089 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00090 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00091 << "| |\n" \ 00092 << "===============================================================================\n"\ 00093 << "| TEST 1: Basis creation, exception testing |\n"\ 00094 << "===============================================================================\n"; 00095 00096 // Define basis and error flag 00097 Basis_HDIV_TET_I1_FEM<double, FieldContainer<double> > tetBasis; 00098 int errorFlag = 0; 00099 00100 // Initialize throw counter for exception testing 00101 int nException = 0; 00102 int throwCounter = 0; 00103 00104 // Define array containing the 4 vertices of the reference TET and its 6 edge midpoints. 00105 FieldContainer<double> tetNodes(10, 3); 00106 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0; 00107 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0; 00108 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0; 00109 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0; 00110 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0; 00111 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0; 00112 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0; 00113 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5; 00114 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5; 00115 tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5; 00116 00117 // Generic array for the output values; needs to be properly resized depending on the operator type 00118 FieldContainer<double> vals; 00119 00120 try{ 00121 // exception #1: GRAD cannot be applied to HDIV functions 00122 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00123 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 3 ); 00124 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException ); 00125 00126 // exception #2: CURL cannot be applied to HDIV functions 00127 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException ); 00128 00129 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00130 // getDofTag() to access invalid array elements thereby causing bounds check exception 00131 // exception #3 00132 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00133 // exception #4 00134 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00135 // exception #5 00136 INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00137 // exception #6 00138 INTREPID_TEST_COMMAND( tetBasis.getDofTag(7), throwCounter, nException ); 00139 // exception #7 00140 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException ); 00141 00142 #ifdef HAVE_INTREPID_DEBUG 00143 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays 00144 // exception #8: input points array must be of rank-2 00145 FieldContainer<double> badPoints1(4, 5, 3); 00146 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00147 00148 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00149 FieldContainer<double> badPoints2(4, 2); 00150 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00151 00152 // exception #10 output values must be of rank-3 for OPERATOR_VALUE 00153 FieldContainer<double> badVals1(4, 3); 00154 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00155 00156 // exception #11 output values must be of rank-2 for OPERATOR_DIV 00157 FieldContainer<double> badVals2(4, 3, 1); 00158 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00159 00160 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00161 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0), 3); 00162 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00163 00164 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions) 00165 FieldContainer<double> badVals4(tetBasis.getCardinality() + 1, tetNodes.dimension(0)); 00166 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException ); 00167 00168 // exception #14 incorrect 1st dimension of output array (must equal number of points) 00169 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0) + 1, 3); 00170 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00171 00172 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00173 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0) + 1); 00174 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException ); 00175 00176 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00177 FieldContainer<double> badVals7(tetBasis.getCardinality(), tetNodes.dimension(0), 4); 00178 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00179 #endif 00180 00181 } 00182 catch (std::logic_error err) { 00183 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00184 *outStream << err.what() << '\n'; 00185 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00186 errorFlag = -1000; 00187 }; 00188 00189 // Check if number of thrown exceptions matches the one we expect 00190 if (throwCounter != nException) { 00191 errorFlag++; 00192 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00193 } 00194 00195 *outStream \ 00196 << "\n" 00197 << "===============================================================================\n"\ 00198 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00199 << "===============================================================================\n"; 00200 00201 try{ 00202 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags(); 00203 00204 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00205 for (unsigned i = 0; i < allTags.size(); i++) { 00206 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00207 00208 std::vector<int> myTag = tetBasis.getDofTag(bfOrd); 00209 if( !( (myTag[0] == allTags[i][0]) && 00210 (myTag[1] == allTags[i][1]) && 00211 (myTag[2] == allTags[i][2]) && 00212 (myTag[3] == allTags[i][3]) ) ) { 00213 errorFlag++; 00214 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00215 *outStream << " getDofOrdinal( {" 00216 << allTags[i][0] << ", " 00217 << allTags[i][1] << ", " 00218 << allTags[i][2] << ", " 00219 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00220 *outStream << " getDofTag(" << bfOrd << ") = { " 00221 << myTag[0] << ", " 00222 << myTag[1] << ", " 00223 << myTag[2] << ", " 00224 << myTag[3] << "}\n"; 00225 } 00226 } 00227 00228 // Now do the same but loop over basis functions 00229 for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) { 00230 std::vector<int> myTag = tetBasis.getDofTag(bfOrd); 00231 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00232 if( bfOrd != myBfOrd) { 00233 errorFlag++; 00234 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00235 *outStream << " getDofTag(" << bfOrd << ") = { " 00236 << myTag[0] << ", " 00237 << myTag[1] << ", " 00238 << myTag[2] << ", " 00239 << myTag[3] << "} but getDofOrdinal({" 00240 << myTag[0] << ", " 00241 << myTag[1] << ", " 00242 << myTag[2] << ", " 00243 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00244 } 00245 } 00246 } 00247 catch (std::logic_error err){ 00248 *outStream << err.what() << "\n\n"; 00249 errorFlag = -1000; 00250 }; 00251 00252 *outStream \ 00253 << "\n" 00254 << "===============================================================================\n"\ 00255 << "| TEST 3: correctness of basis function values |\n"\ 00256 << "===============================================================================\n"; 00257 00258 outStream -> precision(20); 00259 00260 // VALUE: Correct basis values in (P,F,D) format: each row gives the 4x3 correct basis set values 00261 // at an evaluation point. Note that getValues returns results as an (F,P,D) array. 00262 double basisValues[] = { 00263 // 4 vertices 00264 0.,-2.0,0., 0.,0.,0., -2.0,0.,0., 0.,0.,-2.0, 00265 2.0,-2.0,0., 2.0,0.,0., 0.,0.,0., 2.0,0.,-2.0, 00266 0.,0.,0., 0.,2.0,0., -2.0,2.0,0., 0,2.0,-2.0, 00267 0.,-2.0,2.0, 0.,0.,2.0, -2.0,0.,2.0, 0.,0.,0., 00268 // 6 edge midpoints 00269 1.0,-2.0,0., 1.0,0.,0., -1.0,0.,0., 1.0,0.,-2.0, 00270 1.0,-1.0,0., 1.0,1.0,0., -1.0,1.0,0., 1.0,1.0,-2.0, 00271 0.,-1.0,0., 0.,1.0,0., -2.0,1.0,0., 0.,1.0,-2.0, 00272 0.,-2.0,1.0, 0.,0.,1.0, -2.0,0.,1.0, 0.,0.,-1.0, 00273 1.0,-2.0,1.0, 1.0,0.,1.0, -1.0,0.,1.0, 1.0,0.,-1.0, 00274 0.,-1.0,1.0, 0.,1.0,1.0, -2.0,1.0,1.0, 0.,1.0,-1.0 00275 // bf0 bf1 bf2 bf3 00276 }; 00277 00278 // DIV: each row gives the 4 correct values of the divergence of the 4 basis functions 00279 double basisDivs[] = { 00280 // 4 vertices 00281 6.0, 6.0, 6.0, 6.0, 00282 6.0, 6.0, 6.0, 6.0, 00283 6.0, 6.0, 6.0, 6.0, 00284 6.0, 6.0, 6.0, 6.0, 00285 // 6 edge midpoints 00286 6.0, 6.0, 6.0, 6.0, 00287 6.0, 6.0, 6.0, 6.0, 00288 6.0, 6.0, 6.0, 6.0, 00289 6.0, 6.0, 6.0, 6.0, 00290 6.0, 6.0, 6.0, 6.0, 00291 6.0, 6.0, 6.0, 6.0 00292 }; 00293 00294 try{ 00295 00296 // Dimensions for the output arrays: 00297 int numFields = tetBasis.getCardinality(); 00298 int numPoints = tetNodes.dimension(0); 00299 int spaceDim = tetBasis.getBaseCellTopology().getDimension(); 00300 00301 // Generic array for values and curls that will be properly sized before each call 00302 FieldContainer<double> vals; 00303 00304 // Check VALUE of basis functions: resize vals to rank-3 container: 00305 vals.resize(numFields, numPoints, spaceDim); 00306 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE); 00307 for (int i = 0; i < numFields; i++) { 00308 for (int j = 0; j < numPoints; j++) { 00309 for (int k = 0; k < spaceDim; k++) { 00310 // basisValues is (P,F,D) array so its multiindex is (j,i,k) and not (i,j,k)! 00311 int l = k + i * spaceDim + j * spaceDim * numFields; 00312 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00313 errorFlag++; 00314 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00315 00316 // Output the multi-index of the value where the error is: 00317 *outStream << " At (Field,Point,Dim) multi-index { "; 00318 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00319 *outStream << "} computed value: " << vals(i,j,k) 00320 << " but reference value: " << basisValues[l] << "\n"; 00321 } 00322 } 00323 } 00324 } 00325 00326 00327 // Check DIV of basis function: resize vals to rank-2 container 00328 vals.resize(numFields, numPoints); 00329 tetBasis.getValues(vals, tetNodes, OPERATOR_DIV); 00330 for (int i = 0; i < numFields; i++) { 00331 for (int j = 0; j < numPoints; j++) { 00332 int l = i + j * numFields; 00333 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) { 00334 errorFlag++; 00335 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00336 00337 // Output the multi-index of the value where the error is: 00338 *outStream << " At multi-index { "; 00339 *outStream << i << " ";*outStream << j << " "; 00340 *outStream << "} computed divergence component: " << vals(i,j) 00341 << " but reference divergence component: " << basisDivs[l] << "\n"; 00342 } 00343 } 00344 } 00345 00346 } 00347 00348 // Catch unexpected errors 00349 catch (std::logic_error err) { 00350 *outStream << err.what() << "\n\n"; 00351 errorFlag = -1000; 00352 }; 00353 00354 if (errorFlag != 0) 00355 std::cout << "End Result: TEST FAILED\n"; 00356 else 00357 std::cout << "End Result: TEST PASSED\n"; 00358 00359 // reset format state of std::cout 00360 std::cout.copyfmt(oldFormatState); 00361 00362 return errorFlag; 00363 }
1.7.4