|
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_HCURL_HEX_In_FEM.hpp" 00037 #include "Intrepid_PointTools.hpp" 00038 #include "Teuchos_oblackholestream.hpp" 00039 #include "Teuchos_RCP.hpp" 00040 #include "Teuchos_GlobalMPISession.hpp" 00041 00042 using namespace std; 00043 using namespace Intrepid; 00044 00045 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \ 00046 { \ 00047 ++nException; \ 00048 try { \ 00049 S ; \ 00050 } \ 00051 catch (std::logic_error err) { \ 00052 ++throwCounter; \ 00053 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \ 00054 *outStream << err.what() << '\n'; \ 00055 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00056 }; \ 00057 } 00058 00059 int main(int argc, char *argv[]) { 00060 00061 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00062 00063 // This little trick lets us print to std::cout only if 00064 // a (dummy) command-line argument is provided. 00065 int iprint = argc - 1; 00066 Teuchos::RCP<std::ostream> outStream; 00067 Teuchos::oblackholestream bhs; // outputs nothing 00068 if (iprint > 0) 00069 outStream = Teuchos::rcp(&std::cout, false); 00070 else 00071 outStream = Teuchos::rcp(&bhs, false); 00072 00073 // Save the format state of the original std::cout. 00074 Teuchos::oblackholestream oldFormatState; 00075 oldFormatState.copyfmt(std::cout); 00076 00077 *outStream \ 00078 << "===============================================================================\n" \ 00079 << "| |\n" \ 00080 << "| Unit Test (Basis_HCURL_HEX_In_FEM) |\n" \ 00081 << "| |\n" \ 00082 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00083 << "| 2) Basis values for VALUE and CURL operators |\n" \ 00084 << "| |\n" \ 00085 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00086 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00087 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00088 << "| |\n" \ 00089 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00090 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00091 << "| |\n" \ 00092 << "===============================================================================\n"\ 00093 << "| TEST 1: Basis creation, exception testing |\n"\ 00094 << "===============================================================================\n"; 00095 00096 // Define basis and error flag 00097 const int deg = 1; 00098 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); 00099 FieldContainer<double> closedPts(PointTools::getLatticeSize(line,deg),1); 00100 FieldContainer<double> openPts(PointTools::getLatticeSize(line,deg+1,1),1); 00101 PointTools::getLattice<double,FieldContainer<double> >(closedPts,line,deg); 00102 PointTools::getLattice<double,FieldContainer<double> >(openPts,line,deg+1,1); 00103 00104 Basis_HCURL_HEX_In_FEM<double, FieldContainer<double> > hexBasis(deg,closedPts,openPts); 00105 int errorFlag = 0; 00106 00107 // Initialize throw counter for exception testing 00108 int nException = 0; 00109 int throwCounter = 0; 00110 00111 // Define array containing the 8 vertices of the reference HEX, its center and 6 face centers 00112 FieldContainer<double> hexNodes(8, 3); 00113 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0; 00114 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0; 00115 hexNodes(2,0) = -1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0; 00116 hexNodes(3,0) = 1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0; 00117 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0; 00118 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0; 00119 hexNodes(6,0) = -1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0; 00120 hexNodes(7,0) = 1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0; 00121 00122 00123 00124 // Generic array for the output values; needs to be properly resized depending on the operator type 00125 FieldContainer<double> vals; 00126 00127 try{ 00128 // exception #1: GRAD cannot be applied to HCURL functions 00129 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00130 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 ); 00131 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException ); 00132 00133 // exception #2: DIV cannot be applied to HCURL functions 00134 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00135 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) ); 00136 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException ); 00137 00138 // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00139 // getDofTag() to access invalid array elements thereby causing bounds check exception 00140 // exception #3 00141 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00142 // exception #4 00143 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00144 // exception #5 00145 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException ); 00146 // exception #6 00147 INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException ); 00148 // exception #7 00149 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException ); 00150 00151 #ifdef HAVE_INTREPID_DEBUG 00152 // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays 00153 // exception #8: input points array must be of rank-2 00154 FieldContainer<double> badPoints1(4, 5, 3); 00155 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00156 00157 // exception #9 dimension 1 in the input point array must equal space dimension of the cell 00158 FieldContainer<double> badPoints2(4, 2); 00159 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00160 00161 // exception #10 output values must be of rank-3 for OPERATOR_VALUE 00162 FieldContainer<double> badVals1(4, 3); 00163 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException ); 00164 00165 // exception #11 output values must be of rank-3 for OPERATOR_CURL 00166 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException ); 00167 00168 // exception #12 incorrect 0th dimension of output array (must equal number of basis functions) 00169 FieldContainer<double> badVals2(hexBasis.getCardinality() + 1, hexNodes.dimension(0), 3); 00170 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00171 00172 // exception #13 incorrect 1st dimension of output array (must equal number of points) 00173 FieldContainer<double> badVals3(hexBasis.getCardinality(), hexNodes.dimension(0) + 1, 3); 00174 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00175 00176 // exception #14: incorrect 2nd dimension of output array (must equal the space dimension) 00177 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0), 4); 00178 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ; 00179 00180 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension) 00181 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ; 00182 00183 // exception #16: D2 cannot be applied to HCURL functions 00184 // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary) 00185 vals.resize(hexBasis.getCardinality(), 00186 hexNodes.dimension(0), 00187 Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension())); 00188 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException ); 00189 00190 #endif 00191 00192 } 00193 catch (std::logic_error err) { 00194 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00195 *outStream << err.what() << '\n'; 00196 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00197 errorFlag = -1000; 00198 }; 00199 00200 // Check if number of thrown exceptions matches the one we expect 00201 // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match. 00202 if (throwCounter != nException) { 00203 errorFlag++; 00204 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00205 } 00206 //#endif 00207 00208 *outStream \ 00209 << "\n" 00210 << "===============================================================================\n"\ 00211 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00212 << "===============================================================================\n"; 00213 00214 try{ 00215 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags(); 00216 00217 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00218 for (unsigned i = 0; i < allTags.size(); i++) { 00219 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00220 00221 for (unsigned j=0;j<4;j++) std::cout << allTags[i][j] << " "; std::cout << std::endl; 00222 00223 std::vector<int> myTag = hexBasis.getDofTag(bfOrd); 00224 if( !( (myTag[0] == allTags[i][0]) && 00225 (myTag[1] == allTags[i][1]) && 00226 (myTag[2] == allTags[i][2]) && 00227 (myTag[3] == allTags[i][3]) ) ) { 00228 errorFlag++; 00229 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00230 *outStream << " getDofOrdinal( {" 00231 << allTags[i][0] << ", " 00232 << allTags[i][1] << ", " 00233 << allTags[i][2] << ", " 00234 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00235 *outStream << " getDofTag(" << bfOrd << ") = { " 00236 << myTag[0] << ", " 00237 << myTag[1] << ", " 00238 << myTag[2] << ", " 00239 << myTag[3] << "}\n"; 00240 } 00241 } 00242 00243 // Now do the same but loop over basis functions 00244 for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) { 00245 std::vector<int> myTag = hexBasis.getDofTag(bfOrd); 00246 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00247 if( bfOrd != myBfOrd) { 00248 errorFlag++; 00249 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00250 *outStream << " getDofTag(" << bfOrd << ") = { " 00251 << myTag[0] << ", " 00252 << myTag[1] << ", " 00253 << myTag[2] << ", " 00254 << myTag[3] << "} but getDofOrdinal({" 00255 << myTag[0] << ", " 00256 << myTag[1] << ", " 00257 << myTag[2] << ", " 00258 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00259 } 00260 } 00261 } 00262 catch (std::logic_error err){ 00263 *outStream << err.what() << "\n\n"; 00264 errorFlag = -1000; 00265 }; 00266 00267 *outStream \ 00268 << "\n" 00269 << "===============================================================================\n"\ 00270 << "| TEST 3: correctness of basis function values |\n"\ 00271 << "===============================================================================\n"; 00272 00273 outStream -> precision(20); 00274 00275 // VALUE: Each row pair gives the 12x3 correct basis set values at an evaluation point: (P,F,D) layout 00276 double basisValues[] = { 00277 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 00278 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 00279 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 00280 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 00281 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 00282 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 00283 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 00284 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 00285 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 00286 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 00287 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 00288 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1 00289 }; 00290 00291 // CURL 00292 00293 double basisCurls[] = { 00294 0,0.5,-0.5, 0,0.5,-0.5, 0,0,-0.5, 0,0,-0.5, 0,0.5,0, 0,0.5,0, 0,0,0, 0,0,0, 00295 0,0,0.5, 0,0,0.5, 0,0.5,0.5, 0,0.5,0.5, 0,0,0, 0,0,0, 0,0.5,0, 0,0.5,0, 00296 0,-0.5,0, 0,-0.5,0, 0,0,0, 0,0,0, 0,-0.5,-0.5, 0,-0.5,-0.5, 0,0,-0.5, 0,0,-0.5, 00297 0,0,0, 0,0,0, 0,-0.5,0, 0,-0.5,0, 0,0,0.5, 0,0,0.5, 0,-0.5,0.5, 0,-0.5,0.5, 00298 // y-component basis functions 00299 // first y-component bf is (0,1/4(1-x)(1-z),0) 00300 // curl is (1/4(1-x),0,-1/4(1-z)) 00301 0.5,0,-0.5, 0,0,-0.5, 0.5,0,-0.5, 0,0,-0.5, 0.5,0,0, 0,0,0, 0.5,0,0, 0,0,0, 00302 // second y-component bf is (0,1/4(1+x)(1-z),0) 00303 // curl is (1/4(1+x),0,1/4(1-z)) 00304 0,0,0.5, 0.5,0,0.5, 0,0,0.5, 0.5,0,0.5, 0,0,0, 0.5,0,0, 0,0,0, 0.5,0,0, 00305 // third y-component bf is (0,1/4(1-x)(1+z),0) 00306 // curl is (-1/4(1-x),0,-1/4(1+z)) 00307 -0.5,0,0, 0,0,0, -0.5,0,0, 0,0,0, -0.5,0,-0.5, 0,0,-0.5, -0.5,0,-0.5, 0,0,-0.5, 00308 // fourth y-component bf is (0,1/4(1+x)(1+z),0) 00309 // curl is (-1/4(1+x),0,1/4(1+z)) 00310 0.0,0,0, -0.5,0,0, 0.0,0,0, -0.5,0,0, 0.0,0,0.5, -0.5,0,0.5, 0.0,0,0.5, -0.5,0,0.5, 00311 // first z-component bf is (0,0,1/4(1-x)(1-y)) 00312 // curl is (-1/4(1-x),1/4(1-y),0) 00313 -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0, -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0, 00314 // second z-component bf is (0,0,1/4(1+x)(1-y)) 00315 // curl is (-1/4(1+x),1/4(1-y),0) 00316 0.0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0, 0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0, 00317 // third z-component bf is (0,0,1/4(1-x)(1+y)) 00318 // curl is (1/4(1-x),1/4(1+y),0) 00319 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0, 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0, 00320 // fourth z-component bf is (0,0,1/4(1+x)(1+y)) 00321 // curl is (1/4(1+x),-1/4(1+y),0) 00322 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0, 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0 00323 }; 00324 00325 00326 try{ 00327 00328 // Dimensions for the output arrays: 00329 int numFields = hexBasis.getCardinality(); 00330 int numPoints = hexNodes.dimension(0); 00331 int spaceDim = hexBasis.getBaseCellTopology().getDimension(); 00332 00333 // Generic array for values and curls that will be properly sized before each call 00334 FieldContainer<double> vals; 00335 00336 // Check VALUE of basis functions: resize vals to rank-3 container: 00337 vals.resize(numFields, numPoints, spaceDim); 00338 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE); 00339 for (int i = 0; i < numFields; i++) { 00340 for (int j = 0; j < numPoints; j++) { 00341 for (int k = 0; k < spaceDim; k++) { 00342 00343 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k 00344 int l = k + i * spaceDim * numPoints + j * spaceDim; 00345 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) { 00346 errorFlag++; 00347 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00348 00349 // Output the multi-index of the value where the error is: 00350 *outStream << " At multi-index { "; 00351 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00352 *outStream << "} computed value: " << vals(i,j,k) 00353 << " but reference value: " << basisValues[l] << "\n"; 00354 } 00355 } 00356 } 00357 } 00358 00359 // Check CURL of basis function: resize vals to rank-3 container 00360 vals.resize(numFields, numPoints, spaceDim); 00361 hexBasis.getValues(vals, hexNodes, OPERATOR_CURL); 00362 for (int i = 0; i < numFields; i++) { 00363 for (int j = 0; j < numPoints; j++) { 00364 for (int k = 0; k < spaceDim; k++) { 00365 // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k 00366 int l = k + i * spaceDim * numPoints + j * spaceDim; 00367 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) { 00368 errorFlag++; 00369 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00370 00371 // Output the multi-index of the value where the error is: 00372 *outStream << " At multi-index { "; 00373 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00374 *outStream << "} computed curl component: " << vals(i,j,k) 00375 << " but reference curl component: " << basisCurls[l] << "\n"; 00376 } 00377 } 00378 } 00379 } 00380 00381 } 00382 00383 // Catch unexpected errors 00384 catch (std::logic_error err) { 00385 *outStream << err.what() << "\n\n"; 00386 errorFlag = -1000; 00387 }; 00388 00389 if (errorFlag != 0) 00390 std::cout << "End Result: TEST FAILED\n"; 00391 else 00392 std::cout << "End Result: TEST PASSED\n"; 00393 00394 // reset format state of std::cout 00395 std::cout.copyfmt(oldFormatState); 00396 00397 return errorFlag; 00398 }
1.7.4