|
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_HCURL_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_HCURL_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 CURL 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_HCURL_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 // exception #1: GRAD cannot be applied to HCURL functions 00121 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00122 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 ); 00123 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException ); 00124 00125 // exception #2: DIV cannot be applied to HCURL functions 00126 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00127 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) ); 00128 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00129 00130 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00131 // getDofTag() to access invalid array elements thereby causing bounds check exception 00132 // exception #3 00133 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00134 // exception #4 00135 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00136 // exception #5 00137 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00138 // exception #6 00139 INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException ); 00140 // exception #7 00141 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException ); 00142 00143 #ifdef HAVE_INTREPID_DEBUG 00144 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays 00145 // exception #8: input points array must be of rank-2 00146 FieldContainer<double> badPoints1(4, 5, 3); 00147 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00148 00149 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00150 FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1); 00151 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00152 00153 // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D 00154 FieldContainer<double> badVals1(4, 3); 00155 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00156 00157 FieldContainer<double> badCurls1(4,3,2); 00158 // exception #11 output values must be of rank-2 for OPERATOR_CURL 00159 INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException ); 00160 00161 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00162 FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension()); 00163 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00164 00165 // exception #13 incorrect 1st dimension of output array (must equal number of points) 00166 FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() ); 00167 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00168 00169 // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension) 00170 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1); 00171 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00172 00173 // exception #15: D2 cannot be applied to HCURL functions 00174 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00175 vals.resize(quadBasis.getCardinality(), 00176 quadNodes.dimension(0), 00177 Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension())); 00178 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), 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 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match. 00191 if (throwCounter != nException) { 00192 errorFlag++; 00193 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00194 } 00195 //#endif 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 4x2 correct basis set values at an evaluation point: (P,F,D) layout 00263 double basisValues[] = { 00264 0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \ 00265 0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \ 00266 -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \ 00267 0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \ 00268 0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \ 00269 0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \ 00270 -0.250000, 0, 0, -0.125000 00271 }; 00272 00273 // CURL: correct values in (F,P) format 00274 double basisCurls[] = { 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 0.25, 0.25, 0.25, 0.25 00284 }; 00285 00286 00287 try{ 00288 00289 // Dimensions for the output arrays: 00290 int numFields = quadBasis.getCardinality(); 00291 int numPoints = quadNodes.dimension(0); 00292 int spaceDim = quadBasis.getBaseCellTopology().getDimension(); 00293 00294 // Generic array for values and curls that will be properly sized before each call 00295 FieldContainer<double> vals; 00296 00297 // Check VALUE of basis functions: resize vals to rank-3 container: 00298 vals.resize(numFields, numPoints, spaceDim); 00299 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE); 00300 for (int i = 0; i < numFields; i++) { 00301 for (int j = 0; j < numPoints; j++) { 00302 for (int k = 0; k < spaceDim; k++) { 00303 00304 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k 00305 int l = k + i * spaceDim + j * spaceDim * numFields; 00306 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00307 errorFlag++; 00308 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00309 00310 // Output the multi-index of the value where the error is: 00311 *outStream << " At multi-index { "; 00312 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00313 *outStream << "} computed value: " << vals(i,j,k) 00314 << " but reference value: " << basisValues[l] << "\n"; 00315 } 00316 } 00317 } 00318 } 00319 00320 // Check CURL of basis function: resize vals to rank-2 container 00321 vals.resize(numFields, numPoints); 00322 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL); 00323 for (int i = 0; i < numFields; i++) { 00324 for (int j = 0; j < numPoints; j++) { 00325 int l = i + j * numFields; 00326 if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) { 00327 errorFlag++; 00328 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00329 00330 // Output the multi-index of the value where the error is: 00331 *outStream << " At multi-index { "; 00332 *outStream << i << " ";*outStream << j << " "; 00333 *outStream << "} computed curl component: " << vals(i,j) 00334 << " but reference curl component: " << basisCurls[l] << "\n"; 00335 } 00336 } 00337 } 00338 } 00339 00340 // Catch unexpected errors 00341 catch (std::logic_error err) { 00342 *outStream << err.what() << "\n\n"; 00343 errorFlag = -1000; 00344 }; 00345 00346 if (errorFlag != 0) 00347 std::cout << "End Result: TEST FAILED\n"; 00348 else 00349 std::cout << "End Result: TEST PASSED\n"; 00350 00351 // reset format state of std::cout 00352 std::cout.copyfmt(oldFormatState); 00353 00354 return errorFlag; 00355 }
1.7.4