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