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