|
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 00036 #include "Intrepid_FieldContainer.hpp" 00037 #include "Intrepid_HGRAD_LINE_Cn_FEM_JACOBI.hpp" 00038 #include "Intrepid_DefaultCubatureFactory.hpp" 00039 #include "Intrepid_ArrayTools.hpp" 00040 #include "Intrepid_FunctionSpaceTools.hpp" 00041 #include "Teuchos_oblackholestream.hpp" 00042 #include "Teuchos_RCP.hpp" 00043 #include "Teuchos_GlobalMPISession.hpp" 00044 00045 using namespace std; 00046 using namespace Intrepid; 00047 00048 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00049 { \ 00050 ++nException; \ 00051 try { \ 00052 S ; \ 00053 } \ 00054 catch (std::logic_error err) { \ 00055 ++throwCounter; \ 00056 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00057 *outStream << err.what() << '\n'; \ 00058 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00059 }; \ 00060 } 00061 00062 00063 int main(int argc, char *argv[]) { 00064 00065 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00066 00067 // This little trick lets us print to std::cout only if 00068 // a (dummy) command-line argument is provided. 00069 int iprint = argc - 1; 00070 Teuchos::RCP<std::ostream> outStream; 00071 Teuchos::oblackholestream bhs; // outputs nothing 00072 if (iprint > 0) 00073 outStream = Teuchos::rcp(&std::cout, false); 00074 else 00075 outStream = Teuchos::rcp(&bhs, false); 00076 00077 // Save the format state of the original std::cout. 00078 Teuchos::oblackholestream oldFormatState; 00079 oldFormatState.copyfmt(std::cout); 00080 00081 *outStream \ 00082 << "===============================================================================\n" \ 00083 << "| |\n" \ 00084 << "| Unit Test (Basis_HGRAD_LINE_Cn_FEM_JACOBI) |\n" \ 00085 << "| |\n" \ 00086 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00087 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \ 00088 << "| |\n" \ 00089 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00090 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00091 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00092 << "| |\n" \ 00093 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00094 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00095 << "| |\n" \ 00096 << "===============================================================================\n"\ 00097 << "| TEST 1: Basis creation, exception testing |\n"\ 00098 << "===============================================================================\n"; 00099 00100 00101 // Define basis and error flag 00102 double alpha = 0.0, beta = 0.0; 00103 Basis_HGRAD_LINE_Cn_FEM_JACOBI<double, FieldContainer<double> > lineBasis(5, alpha, beta); 00104 int errorFlag = 0; 00105 00106 // Initialize throw counter for exception testing 00107 int nException = 0; 00108 int throwCounter = 0; 00109 00110 // Define array containing vertices of the reference Line and a few other points 00111 int numIntervals = 100; 00112 FieldContainer<double> lineNodes(numIntervals+1, 1); 00113 for (int i=0; i<numIntervals+1; i++) { 00114 lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals; 00115 } 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 // Exceptions 1-5: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00122 // getDofTag() to access invalid array elements thereby causing bounds check exception 00123 // exception #1 00124 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException ); 00125 // exception #2 00126 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00127 // exception #3 00128 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,7), throwCounter, nException ); 00129 // not an exception 00130 INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,5), throwCounter, nException ); --nException; 00131 // exception #4 00132 INTREPID_TEST_COMMAND( lineBasis.getDofTag(6), throwCounter, nException ); 00133 // exception #5 00134 INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException ); 00135 // not an exception 00136 INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException; 00137 #ifdef HAVE_INTREPID_DEBUG 00138 // Exceptions 6-16 test exception handling with incorrectly dimensioned input/output arrays 00139 // exception #6: input points array must be of rank-2 00140 FieldContainer<double> badPoints1(4, 5, 3); 00141 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00142 00143 // exception #7: dimension 1 in the input point array must equal space dimension of the cell 00144 FieldContainer<double> badPoints2(4, 3); 00145 INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00146 00147 // exception #8: output values must be of rank-2 for OPERATOR_VALUE 00148 FieldContainer<double> badVals1(4, 3, 1); 00149 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException ); 00150 00151 // exception #9: output values must be of rank-3 for OPERATOR_GRAD 00152 FieldContainer<double> badVals2(4, 3); 00153 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException ); 00154 00155 // exception #10: output values must be of rank-3 for OPERATOR_CURL 00156 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_CURL), throwCounter, nException ); 00157 00158 // exception #11: output values must be of rank-2 for OPERATOR_DIV 00159 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_DIV), throwCounter, nException ); 00160 00161 // exception #12: output values must be of rank-2 for OPERATOR_D1 00162 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException ); 00163 00164 // exception #13: incorrect 0th dimension of output array (must equal number of basis functions) 00165 FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0)); 00166 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException ); 00167 00168 // exception #14: incorrect 1st dimension of output array (must equal number of points) 00169 FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1); 00170 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException ); 00171 00172 // exception #15: incorrect 2nd dimension of output array (must equal spatial dimension) 00173 FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2); 00174 INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException ); 00175 00176 // not an exception 00177 FieldContainer<double> goodVals2(lineBasis.getCardinality(), lineNodes.dimension(0)); 00178 INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --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! Incorrect number of exceptions." << "\n"; 00193 } 00194 00195 00196 *outStream \ 00197 << "\n" 00198 << "===============================================================================\n"\ 00199 << "| TEST 3: orthogonality of basis functions |\n"\ 00200 << "===============================================================================\n"; 00201 00202 outStream -> precision(20); 00203 00204 try { 00205 00206 // Check orthogonality property for Legendre polynomials. 00207 int maxorder = 10; 00208 00209 DefaultCubatureFactory<double> cubFactory; // create factory 00210 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); // create cell topology 00211 for (int ordi=0; ordi < maxorder; ordi++) { 00212 //create left basis 00213 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisLeft = 00214 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(ordi) ); 00215 00216 for (int ordj=0; ordj < maxorder; ordj++) { 00217 00218 //create right basis 00219 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasisRight = 00220 Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(ordj) ); 00221 00222 // get cubature points and weights 00223 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, ordi+ordj); 00224 int numPoints = lineCub->getNumPoints(); 00225 FieldContainer<double> cubPoints (numPoints, lineCub->getDimension()); 00226 FieldContainer<double> cubWeights(numPoints); 00227 FieldContainer<double> cubWeightsC(1, numPoints); 00228 lineCub->getCubature(cubPoints, cubWeights); 00229 // "reshape" weights 00230 for (int i=0; i<numPoints; i++) { cubWeightsC(0,i) = cubWeights(i); } 00231 00232 00233 // get basis values 00234 int numFieldsLeft = lineBasisLeft ->getCardinality(); 00235 int numFieldsRight = lineBasisRight->getCardinality(); 00236 FieldContainer<double> valsLeft(numFieldsLeft,numPoints), 00237 valsRight(numFieldsRight,numPoints); 00238 lineBasisLeft ->getValues(valsLeft, cubPoints, OPERATOR_VALUE); 00239 lineBasisRight->getValues(valsRight, cubPoints, OPERATOR_VALUE); 00240 00241 // reshape by cloning and integrate 00242 FieldContainer<double> valsLeftC(1, numFieldsLeft,numPoints), 00243 valsRightC(1, numFieldsRight,numPoints), 00244 massMatrix(1, numFieldsLeft, numFieldsRight); 00245 ArrayTools::cloneFields<double>(valsLeftC, valsLeft); 00246 ArrayTools::cloneFields<double>(valsRightC, valsRight); 00247 ArrayTools::scalarMultiplyDataField<double>(valsRightC, cubWeightsC, valsRightC); 00248 FunctionSpaceTools::integrate<double>(massMatrix, valsLeftC, valsRightC, COMP_CPP); 00249 00250 // check orthogonality property 00251 for (int i=0; i<numFieldsLeft; i++) { 00252 for (int j=0; j<numFieldsRight; j++) { 00253 00254 if (i==j) { 00255 if ( std::abs(massMatrix(0,i,j)-(double)(2.0/(2.0*j+1.0))) > INTREPID_TOL ) { 00256 *outStream << "Incorrect ii (\"diagonal\") value for i=" << i << ", j=" << j << ": " 00257 << massMatrix(0,i,j) << " != " << "2/(2*" << j << "+1)\n\n"; 00258 errorFlag++; 00259 } 00260 } 00261 else { 00262 if ( std::abs(massMatrix(0,i,j)) > INTREPID_TOL ) { 00263 *outStream << "Incorrect ij (\"off-diagonal\") value for i=" << i << ", j=" << j << ": " 00264 << massMatrix(0,i,j) << " != " << "0\n\n"; 00265 errorFlag++; 00266 } 00267 } 00268 } 00269 } 00270 00271 } 00272 } 00273 00274 } 00275 // Catch unexpected errors 00276 catch (std::logic_error err) { 00277 *outStream << err.what() << "\n\n"; 00278 errorFlag = -1000; 00279 }; 00280 00281 *outStream \ 00282 << "\n" 00283 << "===============================================================================\n"\ 00284 << "| TEST 4: correctness of basis function derivatives |\n"\ 00285 << "===============================================================================\n"; 00286 00287 outStream -> precision(20); 00288 00289 // function values stored by bf, then pt 00290 double basisValues[] = { 00291 1.000000000000000, 1.000000000000000, 1.000000000000000, \ 00292 1.000000000000000, -1.000000000000000, -0.3333333333333333, \ 00293 0.3333333333333333, 1.000000000000000, 1.000000000000000, \ 00294 -0.3333333333333333, -0.3333333333333333, 1.000000000000000, \ 00295 -1.000000000000000, 0.4074074074074074, -0.4074074074074074, \ 00296 1.000000000000000}; 00297 00298 double basisD1Values[] = 00299 {0, 0, 0, 0, 1.000000000000000, 1.000000000000000, 1.000000000000000, \ 00300 1.000000000000000, -3.000000000000000, -1.000000000000000, \ 00301 1.000000000000000, 3.000000000000000, 6.000000000000000, \ 00302 -0.6666666666666667, -0.6666666666666667, 6.000000000000000}; 00303 00304 double basisD2Values[] = 00305 {0, 0, 0, 0, 0, 0, 0, 0, 3.000000000000000, 3.000000000000000, \ 00306 3.000000000000000, 3.000000000000000, -15.00000000000000, \ 00307 -5.000000000000000, 5.000000000000000, 15.00000000000000}; 00308 00309 double basisD3Values[] = 00310 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.00000000000000, \ 00311 15.00000000000000, 15.00000000000000, 15.00000000000000}; 00312 00313 00314 00315 try { 00316 Basis_HGRAD_LINE_Cn_FEM_JACOBI<double, FieldContainer<double> > lineBasis3(3, alpha, beta); 00317 int numIntervals = 3; 00318 FieldContainer<double> lineNodes3(numIntervals+1, 1); 00319 FieldContainer<double> vals; 00320 for (int i=0; i<numIntervals+1; i++) { 00321 lineNodes3(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals; 00322 } 00323 int numFields = lineBasis3.getCardinality(); 00324 int numPoints = lineNodes3.dimension(0); 00325 00326 // test basis values 00327 vals.resize(numFields, numPoints); 00328 lineBasis3.getValues(vals,lineNodes3,OPERATOR_VALUE); 00329 for (int i = 0; i < numFields; i++) { 00330 for (int j = 0; j < numPoints; j++) { 00331 00332 // Compute offset for (F,P) container 00333 int l = j + i * numPoints; 00334 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00335 errorFlag++; 00336 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00337 00338 // Output the multi-index of the value where the error is: 00339 *outStream << " At multi-index { "; 00340 *outStream << i << " ";*outStream << j << " "; 00341 *outStream << "} computed value: " << vals(i,j) 00342 << " but reference value: " << basisValues[l] << "\n"; 00343 } 00344 } 00345 } 00346 00347 // test basis derivatives 00348 vals.resize(numFields, numPoints,1); 00349 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D1); 00350 for (int i = 0; i < numFields; i++) { 00351 for (int j = 0; j < numPoints; j++) { 00352 00353 // Compute offset for (F,P) container 00354 int l = j + i * numPoints; 00355 if (std::abs(vals(i,j,0) - basisD1Values[l]) > INTREPID_TOL) { 00356 errorFlag++; 00357 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00358 00359 // Output the multi-index of the value where the error is: 00360 *outStream << " At multi-index { "; 00361 *outStream << i << " ";*outStream << j << " "; 00362 *outStream << "} computed value: " << vals(i,j,0) 00363 << " but reference value: " << basisD1Values[l] << "\n"; 00364 } 00365 } 00366 } 00367 00368 vals.resize(numFields, numPoints,1); 00369 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D2); 00370 for (int i = 0; i < numFields; i++) { 00371 for (int j = 0; j < numPoints; j++) { 00372 00373 // Compute offset for (F,P) container 00374 int l = j + i * numPoints; 00375 if (std::abs(vals(i,j,0) - basisD2Values[l]) > INTREPID_TOL) { 00376 errorFlag++; 00377 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00378 00379 // Output the multi-index of the value where the error is: 00380 *outStream << " At multi-index { "; 00381 *outStream << i << " ";*outStream << j << " "; 00382 *outStream << "} computed value: " << vals(i,j,0) 00383 << " but reference value: " << basisD2Values[l] << "\n"; 00384 } 00385 } 00386 } 00387 00388 vals.resize(numFields, numPoints,1); 00389 lineBasis3.getValues(vals,lineNodes3,OPERATOR_D3); 00390 for (int i = 0; i < numFields; i++) { 00391 for (int j = 0; j < numPoints; j++) { 00392 00393 // Compute offset for (F,P) container 00394 int l = j + i * numPoints; 00395 if (std::abs(vals(i,j,0) - basisD3Values[l]) > INTREPID_TOL) { 00396 errorFlag++; 00397 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00398 00399 // Output the multi-index of the value where the error is: 00400 *outStream << " At multi-index { "; 00401 *outStream << i << " ";*outStream << j << " "; 00402 *outStream << "} computed value: " << vals(i,j,0) 00403 << " but reference value: " << basisD3Values[l] << "\n"; 00404 } 00405 } 00406 } 00407 } 00408 // Catch unexpected errors 00409 catch (std::logic_error err) { 00410 *outStream << err.what() << "\n\n"; 00411 errorFlag = -1000; 00412 }; 00413 00414 00415 if (errorFlag != 0) 00416 std::cout << "End Result: TEST FAILED\n"; 00417 else 00418 std::cout << "End Result: TEST PASSED\n"; 00419 00420 // reset format state of std::cout 00421 std::cout.copyfmt(oldFormatState); 00422 00423 return errorFlag; 00424 }
1.7.4