|
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 00035 #include "Intrepid_FieldContainer.hpp" 00036 #include "Intrepid_HGRAD_WEDGE_C1_FEM.hpp" 00037 #include "Teuchos_oblackholestream.hpp" 00038 #include "Teuchos_RCP.hpp" 00039 #include "Teuchos_GlobalMPISession.hpp" 00040 00041 using namespace std; 00042 using namespace Intrepid; 00043 00044 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00045 { \ 00046 ++nException; \ 00047 try { \ 00048 S ; \ 00049 } \ 00050 catch (std::logic_error err) { \ 00051 ++throwCounter; \ 00052 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00053 *outStream << err.what() << '\n'; \ 00054 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00055 }; \ 00056 } 00057 00058 int main(int argc, char *argv[]) { 00059 00060 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00061 00062 // This little trick lets us print to std::cout only if 00063 // a (dummy) command-line argument is provided. 00064 int iprint = argc - 1; 00065 Teuchos::RCP<std::ostream> outStream; 00066 Teuchos::oblackholestream bhs; // outputs nothing 00067 if (iprint > 0) 00068 outStream = Teuchos::rcp(&std::cout, false); 00069 else 00070 outStream = Teuchos::rcp(&bhs, false); 00071 00072 // Save the format state of the original std::cout. 00073 Teuchos::oblackholestream oldFormatState; 00074 oldFormatState.copyfmt(std::cout); 00075 00076 *outStream \ 00077 << "===============================================================================\n" \ 00078 << "| |\n" \ 00079 << "| Unit Test (Basis_HGRAD_WEDGE_C1_FEM) |\n" \ 00080 << "| |\n" \ 00081 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00082 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \ 00083 << "| |\n" \ 00084 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00085 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00086 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00087 << "| |\n" \ 00088 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00089 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00090 << "| |\n" \ 00091 << "===============================================================================\n"\ 00092 << "| TEST 1: Basis creation, exception testing |\n"\ 00093 << "===============================================================================\n"; 00094 00095 // Define basis and error flag 00096 Basis_HGRAD_WEDGE_C1_FEM<double, FieldContainer<double> > wedgeBasis; 00097 int errorFlag = 0; 00098 00099 // Initialize throw counter for exception testing 00100 int nException = 0; 00101 int throwCounter = 0; 00102 00103 // Define array containing the 6 vertices of the reference TRIPRISM and some other points. 00104 FieldContainer<double> wedgeNodes(12, 3); 00105 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0; 00106 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0; 00107 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0; 00108 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0; 00109 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0; 00110 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0; 00111 00112 wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0; 00113 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0; 00114 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.30; wedgeNodes(8,2) = 1.0; 00115 wedgeNodes(9,0) = 0.3; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75; 00116 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.44; wedgeNodes(10,2)= -0.23; 00117 wedgeNodes(11,0)= 0.4; wedgeNodes(11,1)= 0.6; wedgeNodes(11,2)= 0.0; 00118 00119 00120 // Generic array for the output values; needs to be properly resized depending on the operator type 00121 FieldContainer<double> vals; 00122 00123 try{ 00124 // exception #1: CURL cannot be applied to scalar functions 00125 // resize vals to rank-3 container with dimensions (num. points, num. basis functions) 00126 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 ); 00127 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00128 00129 // exception #2: DIV cannot be applied to scalar functions 00130 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00131 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) ); 00132 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException ); 00133 00134 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00135 // getDofTag() to access invalid array elements thereby causing bounds check exception 00136 // exception #3 00137 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00138 // exception #4 00139 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00140 // exception #5 00141 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,6,0), throwCounter, nException ); 00142 // exception #6 00143 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(7), throwCounter, nException ); 00144 // exception #7 00145 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException ); 00146 00147 #ifdef HAVE_INTREPID_DEBUG 00148 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays 00149 // exception #8: input points array must be of rank-2 00150 FieldContainer<double> badPoints1(4, 5, 3); 00151 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00152 00153 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00154 FieldContainer<double> badPoints2(4, 2); 00155 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00156 00157 // exception #10 output values must be of rank-2 for OPERATOR_VALUE 00158 FieldContainer<double> badVals1(4, 3, 1); 00159 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00160 00161 // exception #11 output values must be of rank-3 for OPERATOR_GRAD 00162 FieldContainer<double> badVals2(4, 3); 00163 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException ); 00164 00165 // exception #12 output values must be of rank-3 for OPERATOR_D1 00166 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException ); 00167 00168 // exception #13 output values must be of rank-3 for OPERATOR_D2 00169 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException ); 00170 00171 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) 00172 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0)); 00173 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00174 00175 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00176 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1); 00177 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException ); 00178 00179 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00180 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4); 00181 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException ); 00182 00183 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) 00184 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40); 00185 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException ); 00186 00187 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) 00188 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException ); 00189 #endif 00190 00191 } 00192 catch (std::logic_error err) { 00193 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00194 *outStream << err.what() << '\n'; 00195 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00196 errorFlag = -1000; 00197 }; 00198 00199 // Check if number of thrown exceptions matches the one we expect - 18 00200 if (throwCounter != nException) { 00201 errorFlag++; 00202 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00203 } 00204 00205 *outStream \ 00206 << "\n" 00207 << "===============================================================================\n"\ 00208 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00209 << "===============================================================================\n"; 00210 00211 try{ 00212 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags(); 00213 00214 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00215 for (unsigned i = 0; i < allTags.size(); i++) { 00216 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00217 00218 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00219 if( !( (myTag[0] == allTags[i][0]) && 00220 (myTag[1] == allTags[i][1]) && 00221 (myTag[2] == allTags[i][2]) && 00222 (myTag[3] == allTags[i][3]) ) ) { 00223 errorFlag++; 00224 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00225 *outStream << " getDofOrdinal( {" 00226 << allTags[i][0] << ", " 00227 << allTags[i][1] << ", " 00228 << allTags[i][2] << ", " 00229 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00230 *outStream << " getDofTag(" << bfOrd << ") = { " 00231 << myTag[0] << ", " 00232 << myTag[1] << ", " 00233 << myTag[2] << ", " 00234 << myTag[3] << "}\n"; 00235 } 00236 } 00237 00238 // Now do the same but loop over basis functions 00239 for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) { 00240 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd); 00241 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00242 if( bfOrd != myBfOrd) { 00243 errorFlag++; 00244 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00245 *outStream << " getDofTag(" << bfOrd << ") = { " 00246 << myTag[0] << ", " 00247 << myTag[1] << ", " 00248 << myTag[2] << ", " 00249 << myTag[3] << "} but getDofOrdinal({" 00250 << myTag[0] << ", " 00251 << myTag[1] << ", " 00252 << myTag[2] << ", " 00253 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00254 } 00255 } 00256 } 00257 catch (std::logic_error err){ 00258 *outStream << err.what() << "\n\n"; 00259 errorFlag = -1000; 00260 }; 00261 00262 *outStream \ 00263 << "\n" 00264 << "===============================================================================\n"\ 00265 << "| TEST 3: correctness of basis function values |\n"\ 00266 << "===============================================================================\n"; 00267 00268 outStream -> precision(20); 00269 00270 // VALUE: Each row gives the 4 correct basis set values at an evaluation point 00271 double basisValues[] = { 00272 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00273 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 00274 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 00275 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00276 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00277 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 00278 // 00279 0.25, 0.25, 0.5, 0., 0., 0., 00280 0.125, 0.25, 0.125, 0.125, 0.25, 0.125, 00281 0., 0., 0., 0.45, 0.25, 0.3, 00282 0.0875, 0.0375, 0, 0.6125, 0.2625, 0, 00283 0.3444, 0, 0.2706, 0.2156, 0, 0.1694, 00284 0., 0.2, 0.3, 0., 0.2, 0.3 00285 }; 00286 00287 // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions 00288 double basisGrads[] = { 00289 -1., -1., -0.5, 1., 0, 0, 0, 1., 0, 0., 0., 0.5, 0., 0, 0, 0, 0., 0, \ 00290 -1., -1., 0., 1., 0, -0.5, 0, 1., 0, 0., 0., 0., 0., 0, 0.5, 0, 0., \ 00291 0, -1., -1., 0., 1., 0, 0, 0, 1., -0.5, 0., 0., 0., 0., 0, 0, 0, 0., \ 00292 0.5, 0., 0., -0.5, 0., 0, 0, 0, 0., 0, -1., -1., 0.5, 1., 0, 0, 0, \ 00293 1., 0, 0., 0., 0., 0., 0, -0.5, 0, 0., 0, -1., -1., 0., 1., 0, 0.5, \ 00294 0, 1., 0, 0., 0., 0., 0., 0, 0, 0, 0., -0.5, -1., -1., 0., 1., 0, 0, \ 00295 0, 1., 0.5, -1., -1., -0.125, 1., 0, -0.125, 0, 1., -0.25, 0., 0., \ 00296 0.125, 0., 0, 0.125, 0, 0., 0.25, -0.5, -0.5, -0.125, 0.5, 0, -0.25, \ 00297 0, 0.5, -0.125, -0.5, -0.5, 0.125, 0.5, 0, 0.25, 0, 0.5, 0.125, 0., \ 00298 0., -0.225, 0., 0, -0.125, 0, 0., -0.15, -1., -1., 0.225, 1., 0, \ 00299 0.125, 0, 1., 0.15, -0.125, -0.125, -0.35, 0.125, 0, -0.15, 0, 0.125, \ 00300 0, -0.875, -0.875, 0.35, 0.875, 0, 0.15, 0, 0.875, 0, -0.615, -0.615, \ 00301 -0.28, 0.615, 0, 0, 0, 0.615, -0.22, -0.385, -0.385, 0.28, 0.385, 0, \ 00302 0, 0, 0.385, 0.22, -0.5, -0.5, 0., 0.5, 0, -0.2, 0, 0.5, -0.3, -0.5, \ 00303 -0.5, 0., 0.5, 0, 0.2, 0, 0.5, 0.3 00304 }; 00305 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K) 00306 double basisD2[] = { 00307 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \ 00308 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \ 00309 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \ 00310 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \ 00311 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \ 00312 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \ 00313 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \ 00314 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, \ 00315 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, \ 00316 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, \ 00317 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, \ 00318 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, \ 00319 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, \ 00320 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, \ 00321 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, \ 00322 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \ 00323 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \ 00324 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \ 00325 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \ 00326 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \ 00327 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \ 00328 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \ 00329 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0 00330 }; 00331 try{ 00332 00333 // Dimensions for the output arrays: 00334 int numFields = wedgeBasis.getCardinality(); 00335 int numPoints = wedgeNodes.dimension(0); 00336 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension(); 00337 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00338 00339 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00340 FieldContainer<double> vals; 00341 00342 // Check VALUE of basis functions: resize vals to rank-2 container: 00343 vals.resize(numFields, numPoints); 00344 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE); 00345 for (int i = 0; i < numFields; i++) { 00346 for (int j = 0; j < numPoints; j++) { 00347 int l = i + j * numFields; 00348 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00349 errorFlag++; 00350 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00351 00352 // Output the multi-index of the value where the error is: 00353 *outStream << " At multi-index { "; 00354 *outStream << i << " ";*outStream << j << " "; 00355 *outStream << "} computed value: " << vals(i,j) 00356 << " but reference value: " << basisValues[l] << "\n"; 00357 } 00358 } 00359 } 00360 00361 // Check GRAD of basis function: resize vals to rank-3 container 00362 vals.resize(numFields, numPoints, spaceDim); 00363 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD); 00364 for (int i = 0; i < numFields; i++) { 00365 for (int j = 0; j < numPoints; j++) { 00366 for (int k = 0; k < spaceDim; k++) { 00367 int l = k + i * spaceDim + j * spaceDim * numFields; 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 grad component: " << vals(i,j,k) 00376 << " but reference grad component: " << basisGrads[l] << "\n"; 00377 } 00378 } 00379 } 00380 } 00381 00382 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00383 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1); 00384 for (int i = 0; i < numFields; i++) { 00385 for (int j = 0; j < numPoints; j++) { 00386 for (int k = 0; k < spaceDim; k++) { 00387 int l = k + i * spaceDim + j * spaceDim * numFields; 00388 if (std::abs(vals(i,j,k) - basisGrads[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 << " ";*outStream << k << " "; 00395 *outStream << "} computed D1 component: " << vals(i,j,k) 00396 << " but reference D1 component: " << basisGrads[l] << "\n"; 00397 } 00398 } 00399 } 00400 } 00401 00402 // Check D2 of basis function 00403 vals.resize(numFields, numPoints, D2Cardin); 00404 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2); 00405 for (int i = 0; i < numFields; i++) { 00406 for (int j = 0; j < numPoints; j++) { 00407 for (int k = 0; k < D2Cardin; k++) { 00408 int l = k + i * D2Cardin + j * D2Cardin * numFields; 00409 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00410 errorFlag++; 00411 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00412 00413 // Output the multi-index of the value where the error is: 00414 *outStream << " At multi-index { "; 00415 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00416 *outStream << "} computed D2 component: " << vals(i,j,k) 00417 << " but reference D2 component: " << basisD2[l] << "\n"; 00418 } 00419 } 00420 } 00421 } 00422 00423 // Check all higher derivatives - must be zero. 00424 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) { 00425 00426 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00427 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00428 vals.resize(numFields, numPoints, DkCardin); 00429 00430 wedgeBasis.getValues(vals, wedgeNodes, op); 00431 for (int i = 0; i < vals.size(); i++) { 00432 if (std::abs(vals[i]) > INTREPID_TOL) { 00433 errorFlag++; 00434 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00435 00436 // Get the multi-index of the value where the error is and the operator order 00437 std::vector<int> myIndex; 00438 vals.getMultiIndex(myIndex,i); 00439 int ord = Intrepid::getOperatorOrder(op); 00440 *outStream << " At multi-index { "; 00441 for(int j = 0; j < vals.rank(); j++) { 00442 *outStream << myIndex[j] << " "; 00443 } 00444 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00445 << " but reference D" << ord << " component: 0 \n"; 00446 } 00447 } 00448 } 00449 } 00450 00451 // Catch unexpected errors 00452 catch (std::logic_error err) { 00453 *outStream << err.what() << "\n\n"; 00454 errorFlag = -1000; 00455 }; 00456 00457 if (errorFlag != 0) 00458 std::cout << "End Result: TEST FAILED\n"; 00459 else 00460 std::cout << "End Result: TEST PASSED\n"; 00461 00462 // reset format state of std::cout 00463 std::cout.copyfmt(oldFormatState); 00464 00465 return errorFlag; 00466 }
1.7.4