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