|
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_PYR_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_PYR_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_PYR_C1_FEM<double, FieldContainer<double> > pyrBasis; 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 6 vertices of the reference PYRMID and some other points. 00117 FieldContainer<double> pyrNodes(10, 3); 00118 pyrNodes(0,0) = -1.0; pyrNodes(0,1) = -1.0; pyrNodes(0,2) = 0; 00119 pyrNodes(1,0) = 1.0; pyrNodes(1,1) = -1.0; pyrNodes(1,2) = 0; 00120 pyrNodes(2,0) = 1.0; pyrNodes(2,1) = 1.0; pyrNodes(2,2) = 0; 00121 pyrNodes(3,0) = -1.0; pyrNodes(3,1) = 1.0; pyrNodes(3,2) = 0; 00122 pyrNodes(4,0) = 0.0; pyrNodes(4,1) = 0.0; pyrNodes(4,2) = 1.0; 00123 00124 pyrNodes(5,0) = 0.25; pyrNodes(5,1) = 0.5; pyrNodes(5,2) = 0.2; 00125 pyrNodes(6,0) = -0.7 ; pyrNodes(6,1) = 0.3; pyrNodes(6,2) = 0.3; 00126 pyrNodes(7,0) = 0.; pyrNodes(7,1) = -0.05; pyrNodes(7,2) = 0.95; 00127 pyrNodes(8,0) = -0.15; pyrNodes(8,1) = -0.2; pyrNodes(8,2) = 0.75; 00128 pyrNodes(9,0) = -0.4; pyrNodes(9,1) = 0.9; pyrNodes(9,2) = 0.0; 00129 00130 00131 // Generic array for the output values; needs to be properly resized depending on the operator type 00132 FieldContainer<double> vals; 00133 00134 try{ 00135 // exception #1: CURL cannot be applied to scalar functions 00136 // resize vals to rank-3 container with dimensions (num. points, num. basis functions) 00137 vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0), 3 ); 00138 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_CURL), throwCounter, nException ); 00139 00140 // exception #2: DIV cannot be applied to scalar functions 00141 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00142 vals.resize(pyrBasis.getCardinality(), pyrNodes.dimension(0) ); 00143 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, pyrNodes, OPERATOR_DIV), throwCounter, nException ); 00144 00145 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00146 // getDofTag() to access invalid array elements thereby causing bounds check exception 00147 // exception #3 00148 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00149 // exception #4 00150 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00151 // exception #5 00152 INTREPID_TEST_COMMAND( pyrBasis.getDofOrdinal(0,6,0), throwCounter, nException ); 00153 // exception #6 00154 INTREPID_TEST_COMMAND( pyrBasis.getDofTag(7), throwCounter, nException ); 00155 // exception #7 00156 INTREPID_TEST_COMMAND( pyrBasis.getDofTag(-1), throwCounter, nException ); 00157 00158 #ifdef HAVE_INTREPID_DEBUG 00159 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays 00160 // exception #8: input points array must be of rank-2 00161 FieldContainer<double> badPoints1(4, 5, 3); 00162 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00163 00164 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00165 FieldContainer<double> badPoints2(4, 2); 00166 INTREPID_TEST_COMMAND( pyrBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00167 00168 // exception #10 output values must be of rank-2 for OPERATOR_VALUE 00169 FieldContainer<double> badVals1(4, 3, 1); 00170 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals1, pyrNodes, OPERATOR_VALUE), throwCounter, nException ); 00171 00172 // exception #11 output values must be of rank-3 for OPERATOR_GRAD 00173 FieldContainer<double> badVals2(4, 3); 00174 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_GRAD), throwCounter, nException ); 00175 00176 // exception #12 output values must be of rank-3 for OPERATOR_D1 00177 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_D1), throwCounter, nException ); 00178 00179 // exception #13 output values must be of rank-3 for OPERATOR_D2 00180 //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals2, pyrNodes, OPERATOR_D2), throwCounter, nException ); 00181 00182 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) 00183 FieldContainer<double> badVals3(pyrBasis.getCardinality() + 1, pyrNodes.dimension(0)); 00184 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals3, pyrNodes, OPERATOR_VALUE), throwCounter, nException ); 00185 00186 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00187 FieldContainer<double> badVals4(pyrBasis.getCardinality(), pyrNodes.dimension(0) + 1); 00188 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals4, pyrNodes, OPERATOR_VALUE), throwCounter, nException ); 00189 00190 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00191 FieldContainer<double> badVals5(pyrBasis.getCardinality(), pyrNodes.dimension(0), 4); 00192 INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals5, pyrNodes, OPERATOR_GRAD), throwCounter, nException ); 00193 00194 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) 00195 //FieldContainer<double> badVals6(pyrBasis.getCardinality(), pyrNodes.dimension(0), 40); 00196 //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals6, pyrNodes, OPERATOR_D2), throwCounter, nException ); 00197 00198 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) 00199 //INTREPID_TEST_COMMAND( pyrBasis.getValues(badVals6, pyrNodes, OPERATOR_D3), throwCounter, nException ); 00200 #endif 00201 00202 } 00203 catch (std::logic_error err) { 00204 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00205 *outStream << err.what() << '\n'; 00206 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00207 errorFlag = -1000; 00208 }; 00209 00210 // Check if number of thrown exceptions matches the one we expect - 18 00211 if (throwCounter != nException) { 00212 errorFlag++; 00213 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00214 } 00215 00216 *outStream \ 00217 << "\n" 00218 << "===============================================================================\n"\ 00219 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00220 << "===============================================================================\n"; 00221 00222 try{ 00223 std::vector<std::vector<int> > allTags = pyrBasis.getAllDofTags(); 00224 00225 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00226 for (unsigned i = 0; i < allTags.size(); i++) { 00227 int bfOrd = pyrBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00228 00229 std::vector<int> myTag = pyrBasis.getDofTag(bfOrd); 00230 if( !( (myTag[0] == allTags[i][0]) && 00231 (myTag[1] == allTags[i][1]) && 00232 (myTag[2] == allTags[i][2]) && 00233 (myTag[3] == allTags[i][3]) ) ) { 00234 errorFlag++; 00235 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00236 *outStream << " getDofOrdinal( {" 00237 << allTags[i][0] << ", " 00238 << allTags[i][1] << ", " 00239 << allTags[i][2] << ", " 00240 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00241 *outStream << " getDofTag(" << bfOrd << ") = { " 00242 << myTag[0] << ", " 00243 << myTag[1] << ", " 00244 << myTag[2] << ", " 00245 << myTag[3] << "}\n"; 00246 } 00247 } 00248 00249 // Now do the same but loop over basis functions 00250 for( int bfOrd = 0; bfOrd < pyrBasis.getCardinality(); bfOrd++) { 00251 std::vector<int> myTag = pyrBasis.getDofTag(bfOrd); 00252 int myBfOrd = pyrBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00253 if( bfOrd != myBfOrd) { 00254 errorFlag++; 00255 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00256 *outStream << " getDofTag(" << bfOrd << ") = { " 00257 << myTag[0] << ", " 00258 << myTag[1] << ", " 00259 << myTag[2] << ", " 00260 << myTag[3] << "} but getDofOrdinal({" 00261 << myTag[0] << ", " 00262 << myTag[1] << ", " 00263 << myTag[2] << ", " 00264 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00265 } 00266 } 00267 } 00268 catch (std::logic_error err){ 00269 *outStream << err.what() << "\n\n"; 00270 errorFlag = -1000; 00271 }; 00272 00273 *outStream \ 00274 << "\n" 00275 << "===============================================================================\n"\ 00276 << "| TEST 3: correctness of basis function values |\n"\ 00277 << "===============================================================================\n"; 00278 00279 outStream -> precision(20); 00280 00281 // VALUE: Each row gives the 4 correct basis set values at an evaluation point 00282 double basisValues[] = { 00283 1.0, 0.0, 0.0, 0.0, 0.0, 00284 0.0, 1.0, 0.0, 0.0, 0.0, 00285 0.0, 0.0, 1.0, 0.0, 0.0, 00286 0.0, 0.0, 0.0, 1.0, 0.0, 00287 0.0, 0.0, 0.0, 0.0, 1.0, 00288 // 00289 0.0515625, 0.0984375, 0.4265625, 0.2234375, 0.2, 00290 0.2, 0, 0., 0.5, 0.3, 00291 0.025, 0.025, 0., 0, 0.95, 00292 0.18, 0.045, 0.005, 0.02, 0.75, 00293 0.035, 0.015, 0.285, 0.665, 0., 00294 }; 00295 00296 // GRAD and D1: each row gives the 3 x 5 correct values of the gradients of the 5 basis functions 00297 double basisGrads[] = { 00298 -0.5, -0.5, 0.0, 0.5, 0.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 1.0, \ 00299 -0.5, 0.0, -0.5, 0.5, -0.5, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, \ 00300 0.0, 0.0, 0.0, 0.0, -0.5, -0.5, 0.5, 0.5, 0.0, -0.5, 0.0, -0.5, 0.0, 0.0, 1.0, \ 00301 0.0, -0.5, -0.5, 0.0, 0.0, 0.0, 0.5, 0.0, -0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 1.0, \ 00302 -0.25,-0.25,-0.25, 0.25,-0.25,-0.25, 0.25, 0.25,-0.25,-0.25, 0.25,-0.25, 0.0, 0.0, 1.0, \ 00303 -0.09375, -0.171875, -0.201171875, 0.09375, -0.328125, -0.298828125, 0.40625, 0.328125, -0.201171875, -0.40625, 0.171875, -0.298828125, 0.0, 0.0, 1.0, \ 00304 -0.1428571428571429, -0.5, -0.3571428571428571, 0.1428571428571429, 0.0, -0.1428571428571429, 0.3571428571428571, 0.0, -0.3571428571428571, -0.3571428571428571, 0.5, -0.1428571428571429, 0.0, 0.0, 1.0, \ 00305 -0.5, -0.25, -0.25, 0.5, -0.25, -0.25, 0.0, 0.25, -0.25, 0., 0.25, -0.25, 0.0, 0.0, 1.0, \ 00306 -0.45, -0.4, -0.13, 0.45, -0.1, -0.37, 0.05, 0.1, -0.13, -0.05, 0.4, -0.37, 0.0, 0.0, 1.0, \ 00307 -0.025, -0.35, -0.34, 0.025, -0.15, -0.16, 0.475, 0.15, -0.34, -0.475, 0.35, -0.16, 0.0, 0.0, 1.0 00308 }; 00309 00310 00311 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K) 00312 const double eps = std::numeric_limits<double>::epsilon( ); 00313 double basisD2[] = { 00314 0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0.25,-0.25, 0,-0.25, 0.5, 0,-0.25, 0.25, 0, 0.25,-0.5, 0, 0, 0, 0, 0, 0, \ 00315 0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0.25,-0.25, 0, 0.25,-0.5, 0,-0.25, 0.25, 0,-0.25, 0.5, 0, 0, 0, 0, 0, 0, \ 00316 0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0.25, 0.25, 0, 0.25, 0.5, 0,-0.25,-0.25, 0,-0.25,-0.5, 0, 0, 0, 0, 0, 0, \ 00317 0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0.25, 0.25, 0,-0.25,-0.5, 0,-0.25,-0.25, 0, 0.25, 0.5, 0, 0, 0, 0, 0, 0, \ 00318 0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0.25/eps, 0, 0, 0, 0, 0,-0.25/eps, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, \ 00319 0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0.3125, 0.1953125, 0, 0.09765625, 0.1220703125, 0,-0.3125,-0.1953125, 0,-0.09765625,-0.1220703125, 0, 0, 0, 0, 0, 0, \ 00320 0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428572,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0.3571428571428571, 0.1530612244897959, 0,-0.3571428571428571,-0.306122448979592, 0,-0.3571428571428571,-0.1530612244897959, 0, 0.3571428571428571, 0.306122448979592, 0, 0, 0, 0, 0, 0, \ 00321 0, 5,-5, 0, 0, 0, 0,-5, 5, 0, 0, 0, 0, 5,-5, 0, 0, 0, 0, -5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, \ 00322 0, 1,-0.8, 0,-0.6, 0.96, 0, -1, 0.8, 0, 0.6,-0.96, 0,1,-0.8, 0, -0.6, 0.96, 0, -1, 0.8, 0, 0.6, -0.96, 0, 0, 0, 0, 0, 0, \ 00323 0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225, 0, 0.1, 0.18, 0, 0.25, 0.225, 0,-0.1,-0.18, 0,-0.25,-0.225,0,0.1,0.18, 0, 0, 0, 0, 0, 0 00324 }; 00325 try{ 00326 00327 // Dimensions for the output arrays: 00328 int numFields = pyrBasis.getCardinality(); 00329 int numPoints = pyrNodes.dimension(0); 00330 int spaceDim = pyrBasis.getBaseCellTopology().getDimension(); 00331 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00332 00333 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00334 FieldContainer<double> vals; 00335 00336 // Check VALUE of basis functions: resize vals to rank-2 container: 00337 vals.resize(numFields, numPoints); 00338 pyrBasis.getValues(vals, pyrNodes, OPERATOR_VALUE); 00339 for (int i = 0; i < numFields; i++) { 00340 for (int j = 0; j < numPoints; j++) { 00341 int l = i + j * numFields; 00342 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00343 errorFlag++; 00344 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00345 00346 // Output the multi-index of the value where the error is: 00347 *outStream << " At multi-index { "; 00348 *outStream << i << " ";*outStream << j << " "; 00349 *outStream << "} computed value: " << vals(i,j) 00350 << " but reference value: " << basisValues[l] << "\n"; 00351 } 00352 } 00353 } 00354 00355 // Check GRAD of basis function: resize vals to rank-3 container 00356 vals.resize(numFields, numPoints, spaceDim); 00357 pyrBasis.getValues(vals, pyrNodes, OPERATOR_GRAD); 00358 for (int i = 0; i < numFields; i++) { 00359 for (int j = 0; j < numPoints; j++) { 00360 for (int k = 0; k < spaceDim; k++) { 00361 int l = k + i * spaceDim + j * spaceDim * numFields; 00362 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00363 errorFlag++; 00364 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00365 00366 // Output the multi-index of the value where the error is: 00367 *outStream << " At multi-index { "; 00368 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00369 *outStream << "} computed grad component: " << vals(i,j,k) 00370 << " but reference grad component: " << basisGrads[l] << "\n"; 00371 } 00372 } 00373 } 00374 } 00375 00376 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00377 pyrBasis.getValues(vals, pyrNodes, OPERATOR_D1); 00378 for (int i = 0; i < numFields; i++) { 00379 for (int j = 0; j < numPoints; j++) { 00380 for (int k = 0; k < spaceDim; k++) { 00381 int l = k + i * spaceDim + j * spaceDim * numFields; 00382 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00383 errorFlag++; 00384 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00385 00386 // Output the multi-index of the value where the error is: 00387 *outStream << " At multi-index { "; 00388 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00389 *outStream << "} computed D1 component: " << vals(i,j,k) 00390 << " but reference D1 component: " << basisGrads[l] << "\n"; 00391 } 00392 } 00393 } 00394 } 00395 00396 // Check D2 of basis function 00397 vals.resize(numFields, numPoints, D2Cardin); 00398 pyrBasis.getValues(vals, pyrNodes, OPERATOR_D2); 00399 for (int i = 0; i < numFields; i++) { 00400 for (int j = 0; j < numPoints; j++) { 00401 if(j == 4) continue; //derivatives are singular when z = 1. 00402 for (int k = 0; k < D2Cardin; k++) { 00403 int l = k + i * D2Cardin + j * D2Cardin * numFields; 00404 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00405 errorFlag++; 00406 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00407 00408 // Output the multi-index of the value where the error is: 00409 *outStream << " At multi-index { "; 00410 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00411 *outStream << "} computed D2 component: " << vals(i,j,k) 00412 << " but reference D2 component: " << basisD2[l] << "\n"; 00413 } 00414 } 00415 } 00416 } 00417 } 00418 00419 00420 // Catch unexpected errors 00421 catch (std::logic_error err) { 00422 *outStream << err.what() << "\n\n"; 00423 errorFlag = -1000; 00424 }; 00425 00426 if (errorFlag != 0) 00427 std::cout << "End Result: TEST FAILED\n"; 00428 else 00429 std::cout << "End Result: TEST PASSED\n"; 00430 00431 // reset format state of std::cout 00432 std::cout.copyfmt(oldFormatState); 00433 00434 return errorFlag; 00435 }
1.7.6.1