|
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 00049 #include "Intrepid_FieldContainer.hpp" 00050 #include "Intrepid_HGRAD_TRI_C2_FEM.hpp" 00051 #include "Teuchos_oblackholestream.hpp" 00052 #include "Teuchos_RCP.hpp" 00053 #include "Teuchos_GlobalMPISession.hpp" 00054 00055 using namespace std; 00056 using namespace Intrepid; 00057 00058 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00059 { \ 00060 ++nException; \ 00061 try { \ 00062 S ; \ 00063 } \ 00064 catch (std::logic_error err) { \ 00065 ++throwCounter; \ 00066 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00067 *outStream << err.what() << '\n'; \ 00068 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 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_TRI_C2_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 (dridzal@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 00109 // Define basis and error flag 00110 Basis_HGRAD_TRI_C2_FEM<double, FieldContainer<double> > triBasis; 00111 int errorFlag = 0; 00112 00113 // Initialize throw counter for exception testing 00114 int nException = 0; 00115 int throwCounter = 0; 00116 00117 // The unisolvent set for the default Lagrangian basis of degree 2 uses the standard equispaced 00118 // lattice on Triangle consisting of the 3 vertices and the 3 edge midpoints. To check basis 00119 // correctness it suffices to evaluate each basis function at the lattice points, defined below, 00120 // but we throw in an additional random point. 00121 FieldContainer<double> triNodes(7, 2); 00122 triNodes(0,0) = 0.0; triNodes(0,1) = 0.0; 00123 triNodes(1,0) = 1.0; triNodes(1,1) = 0.0; 00124 triNodes(2,0) = 0.0; triNodes(2,1) = 1.0; 00125 triNodes(3,0) = 0.5; triNodes(3,1) = 0.0; 00126 triNodes(4,0) = 0.5; triNodes(4,1) = 0.5; 00127 triNodes(5,0) = 0.0; triNodes(5,1) = 0.5; 00128 triNodes(6,0) = 0.1732; triNodes(6,1) = 0.6714; 00129 00130 // Generic array for the output values; needs to be properly resized depending on the operator type 00131 FieldContainer<double> vals; 00132 00133 try{ 00134 // exception #1: DIV cannot be applied to scalar functions 00135 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00136 vals.resize(triBasis.getCardinality(), triNodes.dimension(0) ); 00137 INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException ); 00138 00139 // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00140 // getDofTag() to access invalid array elements thereby causing bounds check exception 00141 // exception #2 00142 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException ); 00143 // exception #3 00144 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00145 // exception #4 00146 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException ); 00147 // exception #5 00148 INTREPID_TEST_COMMAND( triBasis.getDofTag(6), throwCounter, nException ); 00149 // exception #6 00150 INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException ); 00151 00152 #ifdef HAVE_INTREPID_DEBUG 00153 // Exceptions 7-17 test exception handling with incorrectly dimensioned input/output arrays 00154 // exception #7: input points array must be of rank-2 00155 FieldContainer<double> badPoints1(4, 5, 3); 00156 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00157 00158 // exception #8 dimension 1 in the input point array must equal space dimension of the cell 00159 FieldContainer<double> badPoints2(4, triBasis.getBaseCellTopology().getDimension() + 1); 00160 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00161 00162 // exception #9 output values must be of rank-2 for OPERATOR_VALUE 00163 FieldContainer<double> badVals1(4, 3, 1); 00164 INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException ); 00165 00166 // exception #10 output values must be of rank-3 for OPERATOR_GRAD 00167 FieldContainer<double> badVals2(4, 3); 00168 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException ); 00169 00170 // exception #11 output values must be of rank-3 for OPERATOR_CURL 00171 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException ); 00172 00173 // exception #12 output values must be of rank-3 for OPERATOR_D2 00174 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException ); 00175 00176 // exception #13 incorrect 1st dimension of output array (must equal number of basis functions) 00177 FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0)); 00178 INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException ); 00179 00180 // exception #14 incorrect 0th dimension of output array (must equal number of points) 00181 FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1); 00182 INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException ); 00183 00184 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension) 00185 FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), triBasis.getBaseCellTopology().getDimension() + 1); 00186 INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException ); 00187 00188 // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D) 00189 FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40); 00190 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException ); 00191 00192 // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D) 00193 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException ); 00194 #endif 00195 00196 } 00197 catch (std::logic_error err) { 00198 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00199 *outStream << err.what() << '\n'; 00200 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00201 errorFlag = -1000; 00202 }; 00203 00204 // Check if number of thrown exceptions matches the one we expect 00205 if (throwCounter != nException) { 00206 errorFlag++; 00207 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00208 } 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 = triBasis.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 = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00222 00223 std::vector<int> myTag = triBasis.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 < triBasis.getCardinality(); bfOrd++) { 00245 std::vector<int> myTag = triBasis.getDofTag(bfOrd); 00246 int myBfOrd = triBasis.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: Correct basis values in (F,P) format: each row gives the 6 correct basis values at 00276 // the lattice points (3 vertices and 3 edge midpoints) 00277 double basisValues[] = { 00278 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.10710168, 00279 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -0.11320352, 00280 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.23015592, 00281 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.10766112, 00282 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.46514592, 00283 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.41734224 00284 }; 00285 00286 // GRAD and D1: Correct gradients and D1 in (F,P,D) format: each row contains 6x2 values of 00287 // gradients of basis functions. Same array is used to check correctness of CURL. 00288 double basisGrads[] = { 00289 // V0 V1 V2 M0 M1 M2 RP 00290 -3.0, -3.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0,-1.0, 0.3784, 0.3784, 00291 -1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, -0.3072, 0.0, 00292 0.0, -1.0, 0.0,-1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.6856, 00293 4.0, 0.0, -4.0,-4.0, 0.0, 0.0, 0.0, -2.0, -2.0,-2.0, 2.0, 0.0, -0.0712, -0.6928, 00294 0.0, 0.0, 0.0, 4.0, 4.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 2.6856, 0.6928, 00295 0.0, 4.0, 0.0, 0.0, -4.0,-4.0, 0.0, 2.0, -2.0,-2.0, -2.0, 0.0, -2.6856, -2.064 00296 }; 00297 00298 // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D 00299 double basisD2[] = { 00300 // V0 V1 V2 M0 M1 M2 RP 00301 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 00302 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 00303 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 00304 -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, 00305 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 00306 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0 00307 }; 00308 try{ 00309 00310 // Dimensions for the output arrays: 00311 int numFields = triBasis.getCardinality(); 00312 int numPoints = triNodes.dimension(0); 00313 int spaceDim = triBasis.getBaseCellTopology().getDimension(); 00314 00315 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00316 FieldContainer<double> vals; 00317 00318 // Check VALUE of basis functions: resize vals to rank-2 container: 00319 vals.resize(numFields, numPoints); 00320 triBasis.getValues(vals, triNodes, OPERATOR_VALUE); 00321 for (int i = 0; i < numFields; i++) { 00322 for (int j = 0; j < numPoints; j++) { 00323 00324 // Compute offset for (F,P) container 00325 int l = j + i * numPoints; 00326 if (std::abs(vals(i,j) - basisValues[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 value: " << vals(i,j) 00334 << " but reference value: " << basisValues[l] << "\n"; 00335 } 00336 } 00337 } 00338 00339 // Check GRAD of basis function: resize vals to rank-3 container 00340 vals.resize(numFields, numPoints, spaceDim); 00341 triBasis.getValues(vals, triNodes, OPERATOR_GRAD); 00342 for (int i = 0; i < numFields; i++) { 00343 for (int j = 0; j < numPoints; j++) { 00344 for (int k = 0; k < spaceDim; k++) { 00345 // basisGrads is (F,P,D), compute offset: 00346 int l = k + j * spaceDim + i * spaceDim * numPoints; 00347 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00348 errorFlag++; 00349 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00350 00351 // Output the multi-index of the value where the error is: 00352 *outStream << " At multi-index { "; 00353 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00354 *outStream << "} computed grad component: " << vals(i,j,k) 00355 << " but reference grad component: " << basisGrads[l] << "\n"; 00356 } 00357 } 00358 } 00359 } 00360 00361 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00362 triBasis.getValues(vals, triNodes, OPERATOR_D1); 00363 for (int i = 0; i < numFields; i++) { 00364 for (int j = 0; j < numPoints; j++) { 00365 for (int k = 0; k < spaceDim; k++) { 00366 // basisGrads is (F,P,D), compute offset: 00367 int l = k + j * spaceDim + i * spaceDim * numPoints; 00368 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00369 errorFlag++; 00370 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00371 00372 // Output the multi-index of the value where the error is: 00373 *outStream << " At multi-index { "; 00374 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00375 *outStream << "} computed D1 component: " << vals(i,j,k) 00376 << " but reference D1 component: " << basisGrads[l] << "\n"; 00377 } 00378 } 00379 } 00380 } 00381 00382 // Check CURL of basis function: resize vals just for illustration! 00383 vals.resize(numFields, numPoints, spaceDim); 00384 triBasis.getValues(vals, triNodes, OPERATOR_CURL); 00385 for (int i = 0; i < numFields; i++) { 00386 for (int j = 0; j < numPoints; j++) { 00387 // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x) 00388 int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints; // position of y-derivative 00389 int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints; // position of x-derivative 00390 00391 double curl_value_0 = basisGrads[curl_0]; 00392 double curl_value_1 =-basisGrads[curl_1]; 00393 if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) { 00394 errorFlag++; 00395 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00396 // Output the multi-index of the value where the error is: 00397 *outStream << " At multi-index { "; 00398 *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " "; 00399 *outStream << "} computed curl component: " << vals(i,j,0) 00400 << " but reference curl component: " << curl_value_0 << "\n"; 00401 } 00402 if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) { 00403 errorFlag++; 00404 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00405 // Output the multi-index of the value where the error is: 00406 *outStream << " At multi-index { "; 00407 *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " "; 00408 *outStream << "} computed curl component: " << vals(i,j,1) 00409 << " but reference curl component: " << curl_value_1 << "\n"; 00410 } 00411 } 00412 } 00413 00414 // Check D2 of basis function 00415 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00416 vals.resize(numFields, numPoints, D2cardinality); 00417 triBasis.getValues(vals, triNodes, OPERATOR_D2); 00418 for (int i = 0; i < numFields; i++) { 00419 for (int j = 0; j < numPoints; j++) { 00420 for (int k = 0; k < D2cardinality; k++) { 00421 00422 // basisD2 is (F,P,Dk), compute offset: 00423 int l = k + j * D2cardinality + i * D2cardinality * numPoints; 00424 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00425 errorFlag++; 00426 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00427 00428 // Output the multi-index of the value where the error is: 00429 *outStream << " At multi-index { "; 00430 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00431 *outStream << "} computed D2 component: " << vals(i,j,k) 00432 << " but reference D2 component: " << basisD2[l] << "\n"; 00433 } 00434 } 00435 } 00436 } 00437 00438 // Check all higher derivatives - must be zero. 00439 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) { 00440 00441 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00442 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00443 vals.resize(numFields, numPoints, DkCardin); 00444 00445 triBasis.getValues(vals, triNodes, op); 00446 for (int i = 0; i < vals.size(); i++) { 00447 if (std::abs(vals[i]) > INTREPID_TOL) { 00448 errorFlag++; 00449 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00450 00451 // Get the multi-index of the value where the error is and the operator order 00452 std::vector<int> myIndex; 00453 vals.getMultiIndex(myIndex,i); 00454 int ord = Intrepid::getOperatorOrder(op); 00455 *outStream << " At multi-index { "; 00456 for(int j = 0; j < vals.rank(); j++) { 00457 *outStream << myIndex[j] << " "; 00458 } 00459 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00460 << " but reference D" << ord << " component: 0 \n"; 00461 } 00462 } 00463 } 00464 } 00465 00466 // Catch unexpected errors 00467 catch (std::logic_error err) { 00468 *outStream << err.what() << "\n\n"; 00469 errorFlag = -1000; 00470 }; 00471 00472 if (errorFlag != 0) 00473 std::cout << "End Result: TEST FAILED\n"; 00474 else 00475 std::cout << "End Result: TEST PASSED\n"; 00476 00477 // reset format state of std::cout 00478 std::cout.copyfmt(oldFormatState); 00479 00480 return errorFlag; 00481 }
1.7.6.1