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