|
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_HEX_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_HEX_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_HEX_C1_FEM<double, FieldContainer<double> > hexBasis; 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 8 vertices of the reference HEX, its center and 6 face centers 00104 FieldContainer<double> hexNodes(15, 3); 00105 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0; 00106 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0; 00107 hexNodes(2,0) = 1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0; 00108 hexNodes(3,0) = -1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0; 00109 00110 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0; 00111 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0; 00112 hexNodes(6,0) = 1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0; 00113 hexNodes(7,0) = -1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0; 00114 00115 hexNodes(8,0) = 0.0; hexNodes(8,1) = 0.0; hexNodes(8,2) = 0.0; 00116 00117 hexNodes(9,0) = 1.0; hexNodes(9,1) = 0.0; hexNodes(9,2) = 0.0; 00118 hexNodes(10,0)= -1.0; hexNodes(10,1)= 0.0; hexNodes(10,2)= 0.0; 00119 00120 hexNodes(11,0)= 0.0; hexNodes(11,1)= 1.0; hexNodes(11,2)= 0.0; 00121 hexNodes(12,0)= 0.0; hexNodes(12,1)= -1.0; hexNodes(12,2)= 0.0; 00122 00123 hexNodes(13,0)= 0.0; hexNodes(13,1)= 0.0; hexNodes(13,2)= 1.0; 00124 hexNodes(14,0)= 0.0; hexNodes(14,1)= 0.0; hexNodes(14,2)= -1.0; 00125 00126 00127 // Generic array for the output values; needs to be properly resized depending on the operator type 00128 FieldContainer<double> vals; 00129 00130 try{ 00131 // exception #1: CURL cannot be applied to scalar functions in 3D 00132 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00133 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 ); 00134 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException ); 00135 00136 // exception #2: DIV cannot be applied to scalar functions in 3D 00137 // resize vals to rank-2 container with dimensions (num. basis functions, num. points) 00138 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) ); 00139 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException ); 00140 00141 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00142 // getDofTag() to access invalid array elements thereby causing bounds check exception 00143 // exception #3 00144 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00145 // exception #4 00146 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00147 // exception #5 00148 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00149 // exception #6 00150 INTREPID_TEST_COMMAND( hexBasis.getDofTag(8), throwCounter, nException ); 00151 // exception #7 00152 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException ); 00153 00154 #ifdef HAVE_INTREPID_DEBUG 00155 // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays 00156 // exception #8: input points array must be of rank-2 00157 FieldContainer<double> badPoints1(4, 5, 3); 00158 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00159 00160 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00161 FieldContainer<double> badPoints2(4, 2); 00162 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00163 00164 // exception #10 output values must be of rank-2 for OPERATOR_VALUE 00165 FieldContainer<double> badVals1(4, 3, 1); 00166 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException ); 00167 00168 // exception #11 output values must be of rank-3 for OPERATOR_GRAD 00169 FieldContainer<double> badVals2(4, 3); 00170 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException ); 00171 00172 // exception #12 output values must be of rank-3 for OPERATOR_D1 00173 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException ); 00174 00175 // exception #13 output values must be of rank-3 for OPERATOR_D2 00176 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException ); 00177 00178 // exception #14 incorrect 0th dimension of output array (must equal number of basis functions) 00179 FieldContainer<double> badVals3(hexBasis.getCardinality() + 1, hexNodes.dimension(0)); 00180 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ); 00181 00182 // exception #15 incorrect 1st dimension of output array (must equal number of points) 00183 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0) + 1); 00184 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ); 00185 00186 // exception #16: incorrect 2nd dimension of output array (must equal the space dimension) 00187 FieldContainer<double> badVals5(hexBasis.getCardinality(), hexNodes.dimension(0), 4); 00188 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException ); 00189 00190 // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D) 00191 FieldContainer<double> badVals6(hexBasis.getCardinality(), hexNodes.dimension(0), 40); 00192 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException ); 00193 00194 // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D) 00195 FieldContainer<double> badVals7(hexBasis.getCardinality(), hexNodes.dimension(0), 50); 00196 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException ); 00197 #endif 00198 00199 } 00200 catch (std::logic_error err) { 00201 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00202 *outStream << err.what() << '\n'; 00203 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00204 errorFlag = -1000; 00205 }; 00206 00207 // Check if number of thrown exceptions matches the one we expect 00208 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match. 00209 if (throwCounter != nException) { 00210 errorFlag++; 00211 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00212 } 00213 00214 *outStream \ 00215 << "\n" 00216 << "===============================================================================\n"\ 00217 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00218 << "===============================================================================\n"; 00219 00220 try{ 00221 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags(); 00222 00223 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00224 for (unsigned i = 0; i < allTags.size(); i++) { 00225 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00226 00227 std::vector<int> myTag = hexBasis.getDofTag(bfOrd); 00228 if( !( (myTag[0] == allTags[i][0]) && 00229 (myTag[1] == allTags[i][1]) && 00230 (myTag[2] == allTags[i][2]) && 00231 (myTag[3] == allTags[i][3]) ) ) { 00232 errorFlag++; 00233 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00234 *outStream << " getDofOrdinal( {" 00235 << allTags[i][0] << ", " 00236 << allTags[i][1] << ", " 00237 << allTags[i][2] << ", " 00238 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00239 *outStream << " getDofTag(" << bfOrd << ") = { " 00240 << myTag[0] << ", " 00241 << myTag[1] << ", " 00242 << myTag[2] << ", " 00243 << myTag[3] << "}\n"; 00244 } 00245 } 00246 00247 // Now do the same but loop over basis functions 00248 for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) { 00249 std::vector<int> myTag = hexBasis.getDofTag(bfOrd); 00250 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00251 if( bfOrd != myBfOrd) { 00252 errorFlag++; 00253 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00254 *outStream << " getDofTag(" << bfOrd << ") = { " 00255 << myTag[0] << ", " 00256 << myTag[1] << ", " 00257 << myTag[2] << ", " 00258 << myTag[3] << "} but getDofOrdinal({" 00259 << myTag[0] << ", " 00260 << myTag[1] << ", " 00261 << myTag[2] << ", " 00262 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00263 } 00264 } 00265 } 00266 catch (std::logic_error err){ 00267 *outStream << err.what() << "\n\n"; 00268 errorFlag = -1000; 00269 }; 00270 00271 *outStream \ 00272 << "\n" 00273 << "===============================================================================\n"\ 00274 << "| TEST 3: correctness of basis function values |\n"\ 00275 << "===============================================================================\n"; 00276 00277 outStream -> precision(20); 00278 00279 // VALUE: Each row gives the 8 correct basis set values at an evaluation point 00280 double basisValues[] = { 00281 // bottom 4 vertices 00282 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00283 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00284 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00285 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 00286 // top 4 vertices 00287 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 00288 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00289 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00290 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 00291 // center {0, 0, 0} 00292 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 00293 // faces { 1, 0, 0} and {-1, 0, 0} 00294 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 00295 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 00296 // faces { 0, 1, 0} and { 0,-1, 0} 00297 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 00298 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 00299 // faces {0, 0, 1} and {0, 0, -1} 00300 0.0, 0.0, 0.0, 0.0, 0.25, 0.25, 0.25, 0.25, 00301 0.25, 0.25, 0.25, 0.25, 0.0, 0.0, 0.0, 0.0, 00302 }; 00303 00304 // GRAD and D1: each row gives the 3x8 correct values of the gradients of the 8 basis functions 00305 double basisGrads[] = { 00306 // points 0-3 00307 -0.5,-0.5,-0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00308 -0.5, 0.0, 0.0, 0.5,-0.5,-0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00309 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5,-0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 00310 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 00311 // points 4-7 00312 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5,-0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 00313 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5, 0.0, 0.0, 0.5,-0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 00314 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5, 0.5, -0.5, 0.0, 0.0, 00315 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5, 0.5, 00316 // point 8 00317 -0.125,-0.125,-0.125, 0.125,-0.125,-0.125, 0.125, 0.125,-0.125, \ 00318 -0.125, 0.125,-0.125, -0.125,-0.125, 0.125, 0.125,-0.125, 0.125, \ 00319 0.125, 0.125, 0.125, -0.125, 0.125, 0.125, 00320 // point 9 00321 -0.125, 0.0, 0.0, 0.125,-0.25, -0.25, 0.125, 0.25, -0.25, -0.125, 0.0, 0.0, \ 00322 -0.125, 0.0, 0.0, 0.125,-0.25, 0.25, 0.125, 0.25, 0.25, -0.125, 0.0, 0.0, 00323 // point 10 00324 -0.125,-0.25, -0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, -0.25,\ 00325 -0.125,-0.25, 0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, 0.25, 00326 // point 11 00327 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125,-0.25, -0.25, 0.125,-0.25,\ 00328 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125, 0.25, -0.25, 0.125, 0.25, 00329 // point 12 00330 -0.25, -0.125,-0.25, 0.25, -0.125,-0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, \ 00331 -0.25, -0.125, 0.25, 0.25, -0.125, 0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 00332 // point 13 00333 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, \ 00334 -0.25, -0.25, 0.125, 0.25, -0.25, 0.125, 0.25, 0.25, 0.125, -0.25, 0.25, 0.125, 00335 // point 14 00336 -0.25, -0.25, -0.125, 0.25, -0.25, -0.125, 0.25, 0.25, -0.125, -0.25, 0.25, -0.125, \ 00337 0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125 00338 }; 00339 00340 //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K) 00341 double basisD2[] = { 00342 // point 0 00343 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0, 0, 0.25, 0., 0, \ 00344 0., 0, 0, -0.25, 0., 0, -0.25, 0, 0, 0., -0.25, 0, -0.25, 0, 0, 0., \ 00345 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0., \ 00346 // point 1 00347 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0, 0, 0.25, 0., 0, \ 00348 -0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, 0, 0., \ 00349 0.25, 0, -0.25, 0, 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0., \ 00350 // Point 2 00351 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0, 0, 0.25, -0.25, 0, \ 00352 -0.25, 0, 0, -0.25, 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., \ 00353 0, -0.25, 0, 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0., \ 00354 // Point 3 00355 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0.25, -0.25, 0, \ 00356 0., 0, 0, -0.25, 0.25, 0, -0.25, 0, 0, 0., 0., 0, -0.25, 0, 0, 0., \ 00357 0., 0, 0., 0, 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0.,\ 00358 // Point 4 00359 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, \ 00360 0, 0., 0., 0, -0.25, 0, 0, 0.25, -0.25, 0, -0.25, 0, 0, -0.25, 0.25, \ 00361 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0., \ 00362 // Point 5 00363 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0, 0, 0., 0., 0, -0.25, \ 00364 0, 0, 0., 0., 0, 0., 0, 0, 0.25, -0.25, 0, 0., 0, 0, -0.25, 0.25, 0, \ 00365 -0.25, 0, 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0., \ 00366 // Point 6 00367 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0, 0, 0., -0.25, 0, -0.25, \ 00368 0, 0, 0., 0.25, 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, \ 00369 -0.25, 0, 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0., \ 00370 // Point 7 00371 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, \ 00372 0, 0., 0.25, 0, -0.25, 0, 0, 0.25, 0., 0, -0.25, 0, 0, -0.25, 0., 0, \ 00373 0., 0, 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0., \ 00374 // Point 8 00375 0, 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0, 0, \ 00376 0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \ 00377 0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \ 00378 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0., \ 00379 // Point 9 00380 0, 0.125, 0.125, 0, 0., 0, 0, -0.125, -0.125, 0, 0.25, 0, 0, 0.125, \ 00381 -0.125, 0, -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, -0.125, 0, \ 00382 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, 0.125, 0, 0.25, 0, 0, \ 00383 -0.125, -0.125, 0, 0., 0., \ 00384 // Point 10 00385 0, 0.125, 0.125, 0, 0.25, 0, 0, -0.125, -0.125, 0, 0., 0, 0, 0.125, \ 00386 -0.125, 0, 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, -0.125, 0, \ 00387 -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, 0.125, 0, 0., 0, 0, \ 00388 -0.125, -0.125, 0, 0.25, 0., \ 00389 // Point 11 00390 0, 0.125, 0., 0, 0.125, 0, 0, -0.125, 0., 0, 0.125, 0, 0, 0.125, \ 00391 -0.25, 0, -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, \ 00392 -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, 0.25, 0, 0.125, 0, \ 00393 0, -0.125, -0.25, 0, 0.125, 0., \ 00394 // Point 12 00395 0, 0.125, 0.25, 0, 0.125, 0, 0, -0.125, -0.25, 0, 0.125, 0, 0, 0.125, \ 00396 0., 0, -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, -0.25, 0, \ 00397 -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, 0.125, 0, \ 00398 0, -0.125, 0., 0, 0.125, 0., \ 00399 // Point 13 00400 0, 0., 0.125, 0, 0.125, 0, 0, 0., -0.125, 0, 0.125, 0, 0, 0., -0.125, \ 00401 0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0.25, -0.125, 0, -0.125, \ 00402 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0.25, 0.125, 0, 0.125, 0, 0, \ 00403 -0.25, -0.125, 0, 0.125, 0., \ 00404 // Point 14 00405 0, 0.25, 0.125, 0, 0.125, 0, 0, -0.25, -0.125, 0, 0.125, 0, 0, 0.25, \ 00406 -0.125, 0, -0.125, 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0., -0.125, \ 00407 0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0., 0.125, 0, 0.125, 0, \ 00408 0, 0., -0.125, 0, 0.125, 0. 00409 }; 00410 00411 try{ 00412 00413 // Dimensions for the output arrays: 00414 int numFields = hexBasis.getCardinality(); 00415 int numPoints = hexNodes.dimension(0); 00416 int spaceDim = hexBasis.getBaseCellTopology().getDimension(); 00417 int D2Cardin = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00418 00419 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00420 FieldContainer<double> vals; 00421 00422 // Check VALUE of basis functions: resize vals to rank-2 container: 00423 vals.resize(numFields, numPoints); 00424 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE); 00425 for (int i = 0; i < numFields; i++) { 00426 for (int j = 0; j < numPoints; j++) { 00427 int l = i + j * numFields; 00428 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00429 errorFlag++; 00430 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00431 00432 // Output the multi-index of the value where the error is: 00433 *outStream << " At multi-index { "; 00434 *outStream << i << " ";*outStream << j << " "; 00435 *outStream << "} computed value: " << vals(i,j) 00436 << " but reference value: " << basisValues[l] << "\n"; 00437 } 00438 } 00439 } 00440 00441 // Check GRAD of basis function: resize vals to rank-3 container 00442 vals.resize(numFields, numPoints, spaceDim); 00443 hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD); 00444 for (int i = 0; i < numFields; i++) { 00445 for (int j = 0; j < numPoints; j++) { 00446 for (int k = 0; k < spaceDim; k++) { 00447 int l = k + i * spaceDim + j * spaceDim * numFields; 00448 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00449 errorFlag++; 00450 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00451 00452 // Output the multi-index of the value where the error is: 00453 *outStream << " At multi-index { "; 00454 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00455 *outStream << "} computed grad component: " << vals(i,j,k) 00456 << " but reference grad component: " << basisGrads[l] << "\n"; 00457 } 00458 } 00459 } 00460 } 00461 00462 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00463 hexBasis.getValues(vals, hexNodes, OPERATOR_D1); 00464 for (int i = 0; i < numFields; i++) { 00465 for (int j = 0; j < numPoints; j++) { 00466 for (int k = 0; k < spaceDim; k++) { 00467 int l = k + i * spaceDim + j * spaceDim * numFields; 00468 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00469 errorFlag++; 00470 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00471 00472 // Output the multi-index of the value where the error is: 00473 *outStream << " At multi-index { "; 00474 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00475 *outStream << "} computed D1 component: " << vals(i,j,k) 00476 << " but reference D1 component: " << basisGrads[l] << "\n"; 00477 } 00478 } 00479 } 00480 } 00481 00482 00483 // Check D2 of basis function 00484 vals.resize(numFields, numPoints, D2Cardin); 00485 hexBasis.getValues(vals, hexNodes, OPERATOR_D2); 00486 for (int i = 0; i < numFields; i++) { 00487 for (int j = 0; j < numPoints; j++) { 00488 for (int k = 0; k < D2Cardin; k++) { 00489 int l = k + i * D2Cardin + j * D2Cardin * numFields; 00490 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00491 errorFlag++; 00492 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00493 00494 // Output the multi-index of the value where the error is: 00495 *outStream << " At multi-index { "; 00496 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00497 *outStream << "} computed D2 component: " << vals(i,j,k) 00498 << " but reference D2 component: " << basisD2[l] << "\n"; 00499 } 00500 } 00501 } 00502 } 00503 00504 // Check all higher derivatives - must be zero. 00505 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) { 00506 00507 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00508 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00509 vals.resize(numFields, numPoints, DkCardin); 00510 00511 hexBasis.getValues(vals, hexNodes, op); 00512 for (int i = 0; i < vals.size(); i++) { 00513 if (std::abs(vals[i]) > INTREPID_TOL) { 00514 errorFlag++; 00515 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00516 00517 // Get the multi-index of the value where the error is and the operator order 00518 std::vector<int> myIndex; 00519 vals.getMultiIndex(myIndex,i); 00520 int ord = Intrepid::getOperatorOrder(op); 00521 *outStream << " At multi-index { "; 00522 for(int j = 0; j < vals.rank(); j++) { 00523 *outStream << myIndex[j] << " "; 00524 } 00525 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00526 << " but reference D" << ord << " component: 0 \n"; 00527 } 00528 } 00529 } 00530 } 00531 00532 // Catch unexpected errors 00533 catch (std::logic_error err) { 00534 *outStream << err.what() << "\n\n"; 00535 errorFlag = -1000; 00536 }; 00537 00538 *outStream \ 00539 << "\n" 00540 << "===============================================================================\n"\ 00541 << "| TEST 4: correctness of DoF locations |\n"\ 00542 << "===============================================================================\n"; 00543 00544 try{ 00545 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis = 00546 Teuchos::rcp(new Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> >); 00547 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface = 00548 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis); 00549 00550 FieldContainer<double> cvals; 00551 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality()); 00552 00553 // Check exceptions. 00554 #ifdef HAVE_INTREPID_DEBUG 00555 cvals.resize(1,2,3); 00556 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00557 cvals.resize(3,2); 00558 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00559 cvals.resize(8,2); 00560 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); 00561 #endif 00562 cvals.resize(8,3); 00563 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--; 00564 // Check if number of thrown exceptions matches the one we expect 00565 if (throwCounter != nException) { 00566 errorFlag++; 00567 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00568 } 00569 00570 // Check mathematical correctness. 00571 basis->getValues(bvals, cvals, OPERATOR_VALUE); 00572 char buffer[120]; 00573 for (int i=0; i<bvals.dimension(0); i++) { 00574 for (int j=0; j<bvals.dimension(1); j++) { 00575 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) { 00576 errorFlag++; 00577 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 0.0); 00578 *outStream << buffer; 00579 } 00580 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) { 00581 errorFlag++; 00582 sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 1.0); 00583 *outStream << buffer; 00584 } 00585 } 00586 } 00587 00588 } 00589 catch (std::logic_error err){ 00590 *outStream << err.what() << "\n\n"; 00591 errorFlag = -1000; 00592 }; 00593 00594 if (errorFlag != 0) 00595 std::cout << "End Result: TEST FAILED\n"; 00596 else 00597 std::cout << "End Result: TEST PASSED\n"; 00598 00599 // reset format state of std::cout 00600 std::cout.copyfmt(oldFormatState); 00601 00602 return errorFlag; 00603 }
1.7.4