|
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_HCURL_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_HCURL_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 CURL 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_HCURL_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 // exception #1: GRAD cannot be applied to HCURL functions 00134 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00135 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 ); 00136 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException ); 00137 00138 // exception #2: DIV cannot be applied to HCURL functions 00139 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00140 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) ); 00141 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00142 00143 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00144 // getDofTag() to access invalid array elements thereby causing bounds check exception 00145 // exception #3 00146 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00147 // exception #4 00148 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00149 // exception #5 00150 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00151 // exception #6 00152 INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException ); 00153 // exception #7 00154 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException ); 00155 00156 #ifdef HAVE_INTREPID_DEBUG 00157 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays 00158 // exception #8: input points array must be of rank-2 00159 FieldContainer<double> badPoints1(4, 5, 3); 00160 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00161 00162 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00163 FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1); 00164 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00165 00166 // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D 00167 FieldContainer<double> badVals1(4, 3); 00168 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00169 00170 FieldContainer<double> badCurls1(4,3,2); 00171 // exception #11 output values must be of rank-2 for OPERATOR_CURL 00172 INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException ); 00173 00174 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00175 FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension()); 00176 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00177 00178 // exception #13 incorrect 1st dimension of output array (must equal number of points) 00179 FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() ); 00180 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00181 00182 // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension) 00183 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1); 00184 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00185 00186 // exception #15: D2 cannot be applied to HCURL functions 00187 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00188 vals.resize(quadBasis.getCardinality(), 00189 quadNodes.dimension(0), 00190 Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension())); 00191 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), 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 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match. 00204 if (throwCounter != nException) { 00205 errorFlag++; 00206 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00207 } 00208 //#endif 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 4x2 correct basis set values at an evaluation point: (P,F,D) layout 00276 double basisValues[] = { 00277 0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \ 00278 0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \ 00279 -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \ 00280 0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \ 00281 0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \ 00282 0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \ 00283 -0.250000, 0, 0, -0.125000 00284 }; 00285 00286 // CURL: correct values in (F,P) format 00287 double basisCurls[] = { 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 0.25, 0.25, 0.25, 0.25 00297 }; 00298 00299 00300 try{ 00301 00302 // Dimensions for the output arrays: 00303 int numFields = quadBasis.getCardinality(); 00304 int numPoints = quadNodes.dimension(0); 00305 int spaceDim = quadBasis.getBaseCellTopology().getDimension(); 00306 00307 // Generic array for values and curls that will be properly sized before each call 00308 FieldContainer<double> vals; 00309 00310 // Check VALUE of basis functions: resize vals to rank-3 container: 00311 vals.resize(numFields, numPoints, spaceDim); 00312 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE); 00313 for (int i = 0; i < numFields; i++) { 00314 for (int j = 0; j < numPoints; j++) { 00315 for (int k = 0; k < spaceDim; k++) { 00316 00317 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k 00318 int l = k + i * spaceDim + j * spaceDim * numFields; 00319 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00320 errorFlag++; 00321 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00322 00323 // Output the multi-index of the value where the error is: 00324 *outStream << " At multi-index { "; 00325 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00326 *outStream << "} computed value: " << vals(i,j,k) 00327 << " but reference value: " << basisValues[l] << "\n"; 00328 } 00329 } 00330 } 00331 } 00332 00333 // Check CURL of basis function: resize vals to rank-2 container 00334 vals.resize(numFields, numPoints); 00335 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL); 00336 for (int i = 0; i < numFields; i++) { 00337 for (int j = 0; j < numPoints; j++) { 00338 int l = i + j * numFields; 00339 if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) { 00340 errorFlag++; 00341 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00342 00343 // Output the multi-index of the value where the error is: 00344 *outStream << " At multi-index { "; 00345 *outStream << i << " ";*outStream << j << " "; 00346 *outStream << "} computed curl component: " << vals(i,j) 00347 << " but reference curl component: " << basisCurls[l] << "\n"; 00348 } 00349 } 00350 } 00351 } 00352 00353 // Catch unexpected errors 00354 catch (std::logic_error err) { 00355 *outStream << err.what() << "\n\n"; 00356 errorFlag = -1000; 00357 }; 00358 00359 *outStream \ 00360 << "\n" 00361 << "===============================================================================\n"\ 00362 << "| TEST 4: correctness of DoF locations |\n"\ 00363 << "===============================================================================\n"; 00364 00365 try{ 00366 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis = 00367 Teuchos::rcp(new Basis_HCURL_QUAD_I1_FEM<double, FieldContainer<double> >); 00368 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface = 00369 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis); 00370 00371 int spaceDim = 2; 00372 FieldContainer<double> cvals; 00373 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),2); // last dimension is spatial dim 00374 00375 // Check exceptions. 00376 #ifdef HAVE_INTREPID_DEBUG 00377 cvals.resize(1,2,3); 00378 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00379 cvals.resize(3,2); 00380 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00381 cvals.resize(4,3); 00382 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00383 #endif 00384 cvals.resize(4,spaceDim); 00385 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--; 00386 // Check if number of thrown exceptions matches the one we expect 00387 if (throwCounter != nException) { 00388 errorFlag++; 00389 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00390 } 00391 00392 // Check mathematical correctness 00393 FieldContainer<double> tangents(basis->getCardinality(),spaceDim); // tangents at each point basis point 00394 tangents(0,0) = 2.0; tangents(0,1) = 0.0; 00395 tangents(1,0) = 0.0; tangents(1,1) = 2.0; 00396 tangents(2,0) = -2.0; tangents(2,1) = 0.0; 00397 tangents(3,0) = 0.0; tangents(3,1) = -2.0; 00398 00399 basis->getValues(bvals, cvals, OPERATOR_VALUE); 00400 char buffer[120]; 00401 for (int i=0; i<bvals.dimension(0); i++) { 00402 for (int j=0; j<bvals.dimension(1); j++) { 00403 00404 double tangent = 0.0; 00405 for(int d=0;d<spaceDim;d++) 00406 tangent += bvals(i,j,d)*tangents(j,d); 00407 00408 if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) { 00409 errorFlag++; 00410 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), tangent, 0.0); 00411 *outStream << buffer; 00412 } 00413 else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) { 00414 errorFlag++; 00415 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), tangent, 1.0); 00416 *outStream << buffer; 00417 } 00418 } 00419 } 00420 00421 } 00422 catch (std::logic_error err){ 00423 *outStream << err.what() << "\n\n"; 00424 errorFlag = -1000; 00425 }; 00426 00427 00428 if (errorFlag != 0) 00429 std::cout << "End Result: TEST FAILED\n"; 00430 else 00431 std::cout << "End Result: TEST PASSED\n"; 00432 00433 // reset format state of std::cout 00434 std::cout.copyfmt(oldFormatState); 00435 00436 return errorFlag; 00437 }
1.7.6.1