|
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_HGRAD_TET_C1_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_HGRAD_TET_C1_FEM) |\n" \ 00080 << "| |\n" \ 00081 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00082 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk 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_HGRAD_TET_C1_FEM<double, FieldContainer<double> > tetBasis; 00097 int errorFlag = 0; 00098 00099 // Initialize throw counter for exception testing 00100 int nException = 0; 00101 int throwCounter = 0; 00102 00103 // Define array containing the 4 vertices of the reference TET and its center. 00104 FieldContainer<double> tetNodes(8, 3); 00105 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0; 00106 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0; 00107 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0; 00108 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0; 00109 tetNodes(4,0) = 0.25; tetNodes(4,1) = 0.5; tetNodes(4,2) = 0.1; 00110 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0; 00111 tetNodes(6,0) = 0.5; tetNodes(6,1) = 0.0; tetNodes(6,2) = 0.5; 00112 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.5; tetNodes(7,2) = 0.5; 00113 00114 00115 // Generic array for the output values; needs to be properly resized depending on the operator type 00116 FieldContainer<double> vals; 00117 00118 try{ 00119 // exception #1: CURL cannot be applied to scalar functions 00120 // resize vals to rank-3 container with dimensions (num. points, num. basis functions, arbitrary) 00121 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 ); 00122 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException ); 00123 00124 // exception #2: DIV cannot be applied to scalar functions 00125 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00126 vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) ); 00127 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), 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,0), throwCounter, nException ); 00137 // exception #6 00138 INTREPID_TEST_COMMAND( tetBasis.getDofTag(5), throwCounter, nException ); 00139 // exception #7 00140 INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException ); 00141 00142 #ifdef HAVE_INTREPID_DEBUG 00143 // Exceptions 8-18 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, tetBasis.getBaseCellTopology().getDimension() - 1); 00150 INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00151 00152 // exception #10 output values must be of rank-2 for OPERATOR_VALUE 00153 FieldContainer<double> badVals1(4, 3, 1); 00154 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00155 00156 // exception #11 output values must be of rank-3 for OPERATOR_GRAD 00157 FieldContainer<double> badVals2(4, 3); 00158 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException ); 00159 00160 // exception #12 output values must be of rank-3 for OPERATOR_D1 00161 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException ); 00162 00163 // exception #13 output values must be of rank-3 for OPERATOR_D2 00164 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException ); 00165 00166 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) 00167 FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0)); 00168 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00169 00170 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00171 FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1); 00172 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException ); 00173 00174 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00175 FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1); 00176 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException ); 00177 00178 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D) 00179 FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40); 00180 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException ); 00181 00182 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D) 00183 INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException ); 00184 #endif 00185 00186 } 00187 catch (std::logic_error err) { 00188 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00189 *outStream << err.what() << '\n'; 00190 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00191 errorFlag = -1000; 00192 }; 00193 00194 // Check if number of thrown exceptions matches the one we expect 00195 if (throwCounter != nException) { 00196 errorFlag++; 00197 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00198 } 00199 00200 *outStream \ 00201 << "\n" 00202 << "===============================================================================\n"\ 00203 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00204 << "===============================================================================\n"; 00205 00206 try{ 00207 std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags(); 00208 00209 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00210 for (unsigned i = 0; i < allTags.size(); i++) { 00211 int bfOrd = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00212 00213 std::vector<int> myTag = tetBasis.getDofTag(bfOrd); 00214 if( !( (myTag[0] == allTags[i][0]) && 00215 (myTag[1] == allTags[i][1]) && 00216 (myTag[2] == allTags[i][2]) && 00217 (myTag[3] == allTags[i][3]) ) ) { 00218 errorFlag++; 00219 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00220 *outStream << " getDofOrdinal( {" 00221 << allTags[i][0] << ", " 00222 << allTags[i][1] << ", " 00223 << allTags[i][2] << ", " 00224 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00225 *outStream << " getDofTag(" << bfOrd << ") = { " 00226 << myTag[0] << ", " 00227 << myTag[1] << ", " 00228 << myTag[2] << ", " 00229 << myTag[3] << "}\n"; 00230 } 00231 } 00232 00233 // Now do the same but loop over basis functions 00234 for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) { 00235 std::vector<int> myTag = tetBasis.getDofTag(bfOrd); 00236 int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00237 if( bfOrd != myBfOrd) { 00238 errorFlag++; 00239 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00240 *outStream << " getDofTag(" << bfOrd << ") = { " 00241 << myTag[0] << ", " 00242 << myTag[1] << ", " 00243 << myTag[2] << ", " 00244 << myTag[3] << "} but getDofOrdinal({" 00245 << myTag[0] << ", " 00246 << myTag[1] << ", " 00247 << myTag[2] << ", " 00248 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00249 } 00250 } 00251 } 00252 catch (std::logic_error err){ 00253 *outStream << err.what() << "\n\n"; 00254 errorFlag = -1000; 00255 }; 00256 00257 *outStream \ 00258 << "\n" 00259 << "===============================================================================\n"\ 00260 << "| TEST 3: correctness of basis function values |\n"\ 00261 << "===============================================================================\n"; 00262 00263 outStream -> precision(20); 00264 00265 // VALUE: Each row gives the 4 correct basis set values at an evaluation point 00266 double basisValues[] = { 00267 1.0, 0.0, 0.0, 0.0, 00268 0.0, 1.0, 0.0, 0.0, 00269 0.0, 0.0, 1.0, 0.0, 00270 0.0, 0.0, 0.0, 1.0, 00271 0.15,0.25,0.5, 0.1, 00272 0.0, 0.5, 0.5, 0.0, 00273 0.0, 0.5, 0.0, 0.5, 00274 0.0, 0.0, 0.5, 0.5 00275 }; 00276 00277 // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions 00278 double basisGrads[] = { 00279 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00280 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00281 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00282 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00283 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00284 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00285 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00286 -1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 00287 }; 00288 00289 try{ 00290 00291 // Dimensions for the output arrays: 00292 int numFields = tetBasis.getCardinality(); 00293 int numPoints = tetNodes.dimension(0); 00294 int spaceDim = tetBasis.getBaseCellTopology().getDimension(); 00295 00296 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00297 FieldContainer<double> vals; 00298 00299 // Check VALUE of basis functions: resize vals to rank-2 container: 00300 vals.resize(numFields, numPoints); 00301 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE); 00302 for (int i = 0; i < numFields; i++) { 00303 for (int j = 0; j < numPoints; j++) { 00304 int l = i + j * numFields; 00305 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00306 errorFlag++; 00307 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00308 00309 // Output the multi-index of the value where the error is: 00310 *outStream << " At multi-index { "; 00311 *outStream << i << " ";*outStream << j << " "; 00312 *outStream << "} computed value: " << vals(i,j) 00313 << " but reference value: " << basisValues[l] << "\n"; 00314 } 00315 } 00316 } 00317 00318 // Check GRAD of basis function: resize vals to rank-3 container 00319 vals.resize(numFields, numPoints, spaceDim); 00320 tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD); 00321 for (int i = 0; i < numFields; i++) { 00322 for (int j = 0; j < numPoints; j++) { 00323 for (int k = 0; k < spaceDim; k++) { 00324 int l = k + i * spaceDim + j * spaceDim * numFields; 00325 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00326 errorFlag++; 00327 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00328 00329 // Output the multi-index of the value where the error is: 00330 *outStream << " At multi-index { "; 00331 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00332 *outStream << "} computed grad component: " << vals(i,j,k) 00333 << " but reference grad component: " << basisGrads[l] << "\n"; 00334 } 00335 } 00336 } 00337 } 00338 00339 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00340 tetBasis.getValues(vals, tetNodes, OPERATOR_D1); 00341 for (int i = 0; i < numFields; i++) { 00342 for (int j = 0; j < numPoints; j++) { 00343 for (int k = 0; k < spaceDim; k++) { 00344 int l = k + i * spaceDim + j * spaceDim * numFields; 00345 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00346 errorFlag++; 00347 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00348 00349 // Output the multi-index of the value where the error is: 00350 *outStream << " At multi-index { "; 00351 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00352 *outStream << "} computed D1 component: " << vals(i,j,k) 00353 << " but reference D1 component: " << basisGrads[l] << "\n"; 00354 } 00355 } 00356 } 00357 } 00358 00359 // Check all higher derivatives - must be zero. 00360 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) { 00361 00362 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00363 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00364 vals.resize(numFields, numPoints, DkCardin); 00365 00366 tetBasis.getValues(vals, tetNodes, op); 00367 for (int i = 0; i < vals.size(); i++) { 00368 if (std::abs(vals[i]) > INTREPID_TOL) { 00369 errorFlag++; 00370 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00371 00372 // Get the multi-index of the value where the error is and the operator order 00373 std::vector<int> myIndex; 00374 vals.getMultiIndex(myIndex,i); 00375 int ord = Intrepid::getOperatorOrder(op); 00376 *outStream << " At multi-index { "; 00377 for(int j = 0; j < vals.rank(); j++) { 00378 *outStream << myIndex[j] << " "; 00379 } 00380 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00381 << " but reference D" << ord << " component: 0 \n"; 00382 } 00383 } 00384 } 00385 } 00386 00387 // Catch unexpected errors 00388 catch (std::logic_error err) { 00389 *outStream << err.what() << "\n\n"; 00390 errorFlag = -1000; 00391 }; 00392 00393 if (errorFlag != 0) 00394 std::cout << "End Result: TEST FAILED\n"; 00395 else 00396 std::cout << "End Result: TEST PASSED\n"; 00397 00398 // reset format state of std::cout 00399 std::cout.copyfmt(oldFormatState); 00400 00401 return errorFlag; 00402 }
1.7.4