|
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_WEDGE_I2_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_WEDGE_I2_FEM) |\n" \ 00093 << "| |\n" \ 00094 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00095 << "| 2) Basis values for VALUE, GRAD, 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_WEDGE_I2_FEM<double, FieldContainer<double> > wedgeBasis; 00110 int errorFlag = 0; 00111 00112 // Initialize throw counter for exception testing 00113 int nException = 0; 00114 int throwCounter = 0; 00115 00116 // Nodes of Wedge<15>: vertices, edge midpoints 00117 FieldContainer<double> wedgeNodes(15, 3); 00118 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0; 00119 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0; 00120 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0; 00121 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0; 00122 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0; 00123 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0; 00124 00125 wedgeNodes(6,0) = 0.5; wedgeNodes(6,1) = 0.0; wedgeNodes(6,2) = -1.0; 00126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.5; wedgeNodes(7,2) = -1.0; 00127 wedgeNodes(8,0) = 0.0; wedgeNodes(8,1) = 0.5; wedgeNodes(8,2) = -1.0; 00128 wedgeNodes(9,0) = 0.0; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.0; 00129 wedgeNodes(10,0)= 1.0; wedgeNodes(10,1)= 0.0; wedgeNodes(10,2)= 0.0; 00130 wedgeNodes(11,0)= 0.0; wedgeNodes(11,1)= 1.0; wedgeNodes(11,2)= 0.0; 00131 wedgeNodes(12,0)= 0.5; wedgeNodes(12,1)= 0.0; wedgeNodes(12,2)= 1.0; 00132 wedgeNodes(13,0)= 0.5; wedgeNodes(13,1)= 0.5; wedgeNodes(13,2)= 1.0; 00133 wedgeNodes(14,0)= 0.0; wedgeNodes(14,1)= 0.5; wedgeNodes(14,2)= 1.0; 00134 00135 00136 // Generic array for the output values; needs to be properly resized depending on the operator type 00137 FieldContainer<double> vals; 00138 00139 try{ 00140 // exception #1: CURL cannot be applied to scalar functions 00141 // resize vals to rank-3 container with dimensions (num. points, num. basis functions) 00142 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 ); 00143 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00144 00145 // exception #2: DIV cannot be applied to scalar functions 00146 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00147 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) ); 00148 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00149 00150 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00151 // getDofTag() to access invalid array elements thereby causing bounds check exception 00152 // exception #3 00153 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00154 // exception #4 00155 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00156 // exception #5 00157 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,9,0), throwCounter, nException ); 00158 // exception #6 00159 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(15), throwCounter, nException ); 00160 // exception #7 00161 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException ); 00162 00163 #ifdef HAVE_INTREPID_DEBUG 00164 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays 00165 // exception #8: input points array must be of rank-2 00166 FieldContainer<double> badPoints1(4, 5, 3); 00167 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00168 00169 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00170 FieldContainer<double> badPoints2(4, wedgeBasis.getBaseCellTopology().getDimension() + 1); 00171 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00172 00173 // exception #10 output values must be of rank-2 for OPERATOR_VALUE 00174 FieldContainer<double> badVals1(4, 3, 1); 00175 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00176 00177 // exception #11 output values must be of rank-3 for OPERATOR_GRAD 00178 FieldContainer<double> badVals2(4, 3); 00179 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException ); 00180 00181 // exception #12 output values must be of rank-3 for OPERATOR_D1 00182 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException ); 00183 00184 // exception #13 output values must be of rank-3 for OPERATOR_D2 00185 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException ); 00186 00187 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) 00188 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0)); 00189 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00190 00191 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00192 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1); 00193 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00194 00195 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00196 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), wedgeBasis.getBaseCellTopology().getDimension() - 1); 00197 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException ); 00198 00199 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) 00200 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40); 00201 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException ); 00202 00203 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) 00204 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException ); 00205 #endif 00206 00207 } 00208 catch (std::logic_error err) { 00209 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00210 *outStream << err.what() << '\n'; 00211 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00212 errorFlag = -1000; 00213 }; 00214 00215 // Check if number of thrown exceptions matches the one we expect - 18 00216 if (throwCounter != nException) { 00217 errorFlag++; 00218 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00219 } 00220 00221 *outStream \ 00222 << "\n" 00223 << "===============================================================================\n"\ 00224 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00225 << "===============================================================================\n"; 00226 00227 try{ 00228 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags(); 00229 00230 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00231 for (unsigned i = 0; i < allTags.size(); i++) { 00232 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00233 00234 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00235 if( !( (myTag[0] == allTags[i][0]) && 00236 (myTag[1] == allTags[i][1]) && 00237 (myTag[2] == allTags[i][2]) && 00238 (myTag[3] == allTags[i][3]) ) ) { 00239 errorFlag++; 00240 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00241 *outStream << " getDofOrdinal( {" 00242 << allTags[i][0] << ", " 00243 << allTags[i][1] << ", " 00244 << allTags[i][2] << ", " 00245 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00246 *outStream << " getDofTag(" << bfOrd << ") = { " 00247 << myTag[0] << ", " 00248 << myTag[1] << ", " 00249 << myTag[2] << ", " 00250 << myTag[3] << "}\n"; 00251 } 00252 } 00253 00254 // Now do the same but loop over basis functions 00255 for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) { 00256 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00257 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00258 if( bfOrd != myBfOrd) { 00259 errorFlag++; 00260 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00261 *outStream << " getDofTag(" << bfOrd << ") = { " 00262 << myTag[0] << ", " 00263 << myTag[1] << ", " 00264 << myTag[2] << ", " 00265 << myTag[3] << "} but getDofOrdinal({" 00266 << myTag[0] << ", " 00267 << myTag[1] << ", " 00268 << myTag[2] << ", " 00269 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00270 } 00271 } 00272 } 00273 catch (std::logic_error err){ 00274 *outStream << err.what() << "\n\n"; 00275 errorFlag = -1000; 00276 }; 00277 00278 *outStream \ 00279 << "\n" 00280 << "===============================================================================\n"\ 00281 << "| TEST 3: correctness of basis function values |\n"\ 00282 << "===============================================================================\n"; 00283 00284 outStream -> precision(20); 00285 00286 // VALUE: correct basis function values in (F,P) format 00287 double basisValues[] = { 00288 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00289 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00290 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00291 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00292 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00293 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00294 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00295 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00296 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00297 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,\ 00298 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,\ 00299 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,\ 00300 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,\ 00301 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,\ 00302 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}; 00303 00304 // GRAD, D1, D2, and D3 test values are stored in files due to their large size 00305 std::string fileName; 00306 std::ifstream dataFile; 00307 00308 // GRAD and D1 values are stored in (F,P,D) format in a data file. Read file and do the test 00309 std::vector<double> basisGrads; // Flat array for the gradient values. 00310 00311 fileName = "./testdata/WEDGE_I2_GradVals.dat"; 00312 dataFile.open(fileName.c_str()); 00313 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error, 00314 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open GRAD values data file, test aborted."); 00315 while (!dataFile.eof() ){ 00316 double temp; 00317 string line; // string for one line of input file 00318 std::getline(dataFile, line); // get next line from file 00319 stringstream data_line(line); // convert to stringstream 00320 while(data_line >> temp){ // extract value from line 00321 basisGrads.push_back(temp); // push into vector 00322 } 00323 } 00324 // It turns out that just closing and then opening the ifstream variable does not reset it 00325 // and subsequent open() command fails. One fix is to explicitely clear the ifstream, or 00326 // scope the variables. 00327 dataFile.close(); 00328 dataFile.clear(); 00329 00330 00331 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (F,P,D2cardinality) 00332 std::vector<double> basisD2; 00333 00334 fileName = "./testdata/WEDGE_I2_D2Vals.dat"; 00335 dataFile.open(fileName.c_str()); 00336 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error, 00337 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D2 values data file, test aborted."); 00338 00339 while (!dataFile.eof() ){ 00340 double temp; 00341 string line; // string for one line of input file 00342 std::getline(dataFile, line); // get next line from file 00343 stringstream data_line(line); // convert to stringstream 00344 while(data_line >> temp){ // extract value from line 00345 basisD2.push_back(temp); // push into vector 00346 } 00347 } 00348 dataFile.close(); 00349 dataFile.clear(); 00350 00351 00352 //D3: flat array with the values of D3 applied to basis functions. Multi-index is (F,P,D3cardinality) 00353 std::vector<double> basisD3; 00354 00355 fileName = "./testdata/WEDGE_I2_D3Vals.dat"; 00356 dataFile.open(fileName.c_str()); 00357 TEUCHOS_TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error, 00358 ">>> ERROR (HGRAD_WEDGE_I2/test01): could not open D3 values data file, test aborted."); 00359 00360 while (!dataFile.eof() ){ 00361 double temp; 00362 string line; // string for one line of input file 00363 std::getline(dataFile, line); // get next line from file 00364 stringstream data_line(line); // convert to stringstream 00365 while(data_line >> temp){ // extract value from line 00366 basisD3.push_back(temp); // push into vector 00367 } 00368 } 00369 dataFile.close(); 00370 dataFile.clear(); 00371 00372 try{ 00373 00374 // Dimensions for the output arrays: 00375 int numFields = wedgeBasis.getCardinality(); 00376 int numPoints = wedgeNodes.dimension(0); 00377 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension(); 00378 00379 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00380 FieldContainer<double> vals; 00381 00382 // Check VALUE of basis functions: resize vals to rank-2 container: 00383 vals.resize(numFields, numPoints); 00384 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE); 00385 for (int i = 0; i < numFields; i++) { 00386 for (int j = 0; j < numPoints; j++) { 00387 int l = i + j * numFields; 00388 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00389 errorFlag++; 00390 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00391 00392 // Output the multi-index of the value where the error is: 00393 *outStream << " At multi-index { "; 00394 *outStream << i << " ";*outStream << j << " "; 00395 *outStream << "} computed value: " << vals(i,j) 00396 << " but reference value: " << basisValues[l] << "\n"; 00397 } 00398 } 00399 } 00400 00401 00402 00403 // Check GRAD of basis function: resize vals to rank-3 container 00404 vals.resize(numFields, numPoints, spaceDim); 00405 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD); 00406 for (int i = 0; i < numFields; i++) { 00407 for (int j = 0; j < numPoints; j++) { 00408 for (int k = 0; k < spaceDim; k++) { 00409 00410 // basisGrads is (F,P,D), compute offset: 00411 int l = k + j * spaceDim + i * spaceDim * numPoints; 00412 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00413 errorFlag++; 00414 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00415 00416 // Output the multi-index of the value where the error is: 00417 *outStream << " At multi-index { "; 00418 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00419 *outStream << "} computed grad component: " << vals(i,j,k) 00420 << " but reference grad component: " << basisGrads[l] << "\n"; 00421 } 00422 } 00423 } 00424 } 00425 00426 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00427 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1); 00428 for (int i = 0; i < numFields; i++) { 00429 for (int j = 0; j < numPoints; j++) { 00430 for (int k = 0; k < spaceDim; k++) { 00431 00432 // basisGrads is (F,P,D), compute offset: 00433 int l = k + j * spaceDim + i * spaceDim * numPoints; 00434 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00435 errorFlag++; 00436 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00437 00438 // Output the multi-index of the value where the error is: 00439 *outStream << " At multi-index { "; 00440 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00441 *outStream << "} computed D1 component: " << vals(i,j,k) 00442 << " but reference D1 component: " << basisGrads[l] << "\n"; 00443 } 00444 } 00445 } 00446 } 00447 00448 00449 // Check D2 of basis function 00450 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00451 vals.resize(numFields, numPoints, D2cardinality); 00452 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2); 00453 for (int i = 0; i < numFields; i++) { 00454 for (int j = 0; j < numPoints; j++) { 00455 for (int k = 0; k < D2cardinality; k++) { 00456 00457 // basisD2 is (F,P,Dk), compute offset: 00458 int l = k + j * D2cardinality + i * D2cardinality * numPoints; 00459 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00460 errorFlag++; 00461 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00462 00463 // Output the multi-index of the value where the error is: 00464 *outStream << " At multi-index { "; 00465 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00466 *outStream << "} computed D2 component: " << vals(i,j,k) 00467 << " but reference D2 component: " << basisD2[l] << "\n"; 00468 } 00469 } 00470 } 00471 } 00472 00473 00474 // Check D3 of basis function 00475 int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim); 00476 vals.resize(numFields, numPoints, D3cardinality); 00477 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D3); 00478 00479 for (int i = 0; i < numFields; i++) { 00480 for (int j = 0; j < numPoints; j++) { 00481 for (int k = 0; k < D3cardinality; k++) { 00482 00483 // basisD3 is (F,P,Dk), compute offset: 00484 int l = k + j * D3cardinality + i * D3cardinality * numPoints; 00485 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) { 00486 errorFlag++; 00487 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00488 00489 // Output the multi-index of the value where the error is: 00490 *outStream << " At multi-index { "; 00491 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00492 *outStream << "} computed D3 component: " << vals(i,j,k) 00493 << " but reference D3 component: " << basisD3[l] << "\n"; 00494 } 00495 } 00496 } 00497 } 00498 00499 00500 // Check all higher derivatives - must be zero. 00501 for(EOperator op = OPERATOR_D4; op < OPERATOR_MAX; op++) { 00502 00503 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00504 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00505 vals.resize(numFields, numPoints, DkCardin); 00506 00507 wedgeBasis.getValues(vals, wedgeNodes, op); 00508 for (int i = 0; i < vals.size(); i++) { 00509 if (std::abs(vals[i]) > INTREPID_TOL) { 00510 errorFlag++; 00511 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00512 00513 // Get the multi-index of the value where the error is and the operator order 00514 std::vector<int> myIndex; 00515 vals.getMultiIndex(myIndex,i); 00516 int ord = Intrepid::getOperatorOrder(op); 00517 *outStream << " At multi-index { "; 00518 for(int j = 0; j < vals.rank(); j++) { 00519 *outStream << myIndex[j] << " "; 00520 } 00521 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00522 << " but reference D" << ord << " component: 0 \n"; 00523 } 00524 } 00525 } 00526 } 00527 00528 // Catch unexpected errors 00529 catch (std::logic_error err) { 00530 *outStream << err.what() << "\n\n"; 00531 errorFlag = -1000; 00532 }; 00533 00534 if (errorFlag != 0) 00535 std::cout << "End Result: TEST FAILED\n"; 00536 else 00537 std::cout << "End Result: TEST PASSED\n"; 00538 00539 // reset format state of std::cout 00540 std::cout.copyfmt(oldFormatState); 00541 00542 return errorFlag; 00543 }
1.7.6.1