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