|
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_QUAD_Cn_FEM.hpp" 00050 #include "Teuchos_oblackholestream.hpp" 00051 #include "Teuchos_RCP.hpp" 00052 #include "Teuchos_GlobalMPISession.hpp" 00053 #include "Intrepid_PointTools.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_HGRAD_QUAD_Cn_FEM) |\n" \ 00094 << "| |\n" \ 00095 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \ 00096 << "| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \ 00097 << "| |\n" \ 00098 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00099 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \ 00100 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00101 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00102 << "| |\n" \ 00103 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00104 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00105 << "| |\n" \ 00106 << "===============================================================================\n"\ 00107 << "| TEST 1: Basis creation, exception testing |\n"\ 00108 << "===============================================================================\n"; 00109 00110 // Define basis and error flag 00111 // get points for basis 00112 const int deg=2; 00113 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); 00114 FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1); 00115 PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg); 00116 00117 Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadBasis(deg,deg,pts,pts); 00118 int errorFlag = 0; 00119 00120 // Initialize throw counter for exception testing 00121 int nException = 0; 00122 int throwCounter = 0; 00123 00124 // Array with the 4 vertices, 4 edge midpoints, center of the reference QUAD and a random point. 00125 FieldContainer<double> quadNodes(10, 2); 00126 quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0; 00127 quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0; 00128 quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0; 00129 quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0; 00130 // edge midpoints 00131 quadNodes(4,0) = 0.0; quadNodes(4,1) = -1.0; 00132 quadNodes(5,0) = 1.0; quadNodes(5,1) = 0.0; 00133 quadNodes(6,0) = 0.0; quadNodes(6,1) = 1.0; 00134 quadNodes(7,0) = -1.0; quadNodes(7,1) = 0.0; 00135 // center & random point 00136 quadNodes(8,0) = 0.0; quadNodes(8,1) = 0.0; 00137 quadNodes(9,0) =1./3.; quadNodes(9,1) =-3./5.; 00138 00139 // Generic array for the output values; needs to be properly resized depending on the operator type 00140 FieldContainer<double> vals; 00141 00142 try{ 00143 // exception #1: DIV cannot be applied to scalar functions 00144 // resize vals to rank-2 container with dimensions (num. points, num. basis functions) 00145 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0)); 00146 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException ); 00147 00148 // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 00149 // getDofTag() to access invalid array elements thereby causing bounds check exception 00150 // exception #2 00151 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException ); 00152 // exception #3 00153 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException ); 00154 // exception #4 00155 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException ); 00156 // exception #5 00157 INTREPID_TEST_COMMAND( quadBasis.getDofTag(10), throwCounter, nException ); 00158 // exception #6 00159 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException ); 00160 00161 #ifdef HAVE_INTREPID_DEBUG 00162 // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays 00163 // exception #7: input points array must be of rank-2 00164 FieldContainer<double> badPoints1(4, 5, 3); 00165 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException ); 00166 00167 // exception #8 dimension 1 in the input point array must equal space dimension of the cell 00168 FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1); 00169 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException ); 00170 00171 // exception #9 output values must be of rank-2 for OPERATOR_VALUE 00172 FieldContainer<double> badVals1(4, 3, 1); 00173 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00174 00175 // exception #10 output values must be of rank-3 for OPERATOR_GRAD 00176 FieldContainer<double> badVals2(4, 3); 00177 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException ); 00178 00179 // exception #11 output values must be of rank-3 for OPERATOR_CURL 00180 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException ); 00181 00182 // exception #12 output values must be of rank-3 for OPERATOR_D2 00183 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException ); 00184 00185 // exception #13 incorrect 0th dimension of output array (must equal number of basis functions) 00186 FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0)); 00187 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00188 00189 // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes) 00190 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1); 00191 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ); 00192 00193 // exception #15: incorrect 2nd dimension of output array (must equal the space dimension) 00194 FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() + 1); 00195 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException ); 00196 00197 // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D) 00198 FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40); 00199 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException ); 00200 00201 // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D) 00202 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException ); 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 if (throwCounter != nException) { 00215 errorFlag++; 00216 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00217 } 00218 00219 *outStream \ 00220 << "\n" 00221 << "===============================================================================\n"\ 00222 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\ 00223 << "===============================================================================\n"; 00224 00225 try{ 00226 std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags(); 00227 00228 // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again 00229 for (unsigned i = 0; i < allTags.size(); i++) { 00230 int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]); 00231 00232 std::vector<int> myTag = quadBasis.getDofTag(bfOrd); 00233 if( !( (myTag[0] == allTags[i][0]) && 00234 (myTag[1] == allTags[i][1]) && 00235 (myTag[2] == allTags[i][2]) && 00236 (myTag[3] == allTags[i][3]) ) ) { 00237 errorFlag++; 00238 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00239 *outStream << " getDofOrdinal( {" 00240 << allTags[i][0] << ", " 00241 << allTags[i][1] << ", " 00242 << allTags[i][2] << ", " 00243 << allTags[i][3] << "}) = " << bfOrd <<" but \n"; 00244 *outStream << " getDofTag(" << bfOrd << ") = { " 00245 << myTag[0] << ", " 00246 << myTag[1] << ", " 00247 << myTag[2] << ", " 00248 << myTag[3] << "}\n"; 00249 } 00250 } 00251 00252 // Now do the same but loop over basis functions 00253 for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) { 00254 std::vector<int> myTag = quadBasis.getDofTag(bfOrd); 00255 int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]); 00256 if( bfOrd != myBfOrd) { 00257 errorFlag++; 00258 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00259 *outStream << " getDofTag(" << bfOrd << ") = { " 00260 << myTag[0] << ", " 00261 << myTag[1] << ", " 00262 << myTag[2] << ", " 00263 << myTag[3] << "} but getDofOrdinal({" 00264 << myTag[0] << ", " 00265 << myTag[1] << ", " 00266 << myTag[2] << ", " 00267 << myTag[3] << "} ) = " << myBfOrd << "\n"; 00268 } 00269 } 00270 } 00271 catch (std::logic_error err){ 00272 *outStream << err.what() << "\n\n"; 00273 errorFlag = -1000; 00274 }; 00275 00276 *outStream \ 00277 << "\n" 00278 << "===============================================================================\n"\ 00279 << "| TEST 3: correctness of basis function values |\n"\ 00280 << "===============================================================================\n"; 00281 00282 outStream -> precision(20); 00283 00284 // VALUE: Correct basis values in (F,P) format: 00285 double basisValues[] = { 00286 1, 0, 0, 0, 0, 0, 0, 0, 0, -0.05333333333333334, \ 00287 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.4266666666666667, \ 00288 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.1066666666666667, \ 00289 0, 0, 0, 0, 0, 0, 0, 1, 0, -0.07111111111111112 , \ 00290 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5688888888888890, \ 00291 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.1422222222222222 ,\ 00292 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.01333333333333333, \ 00293 0, 0, 0, 0, 0, 0, 1, 0, 0, -0.1066666666666667, \ 00294 0, 0, 1, 0, 0, 0, 0, 0, 0, -0.02666666666666666 }; 00295 00296 // FIXME HERE: needs to be reordered. 00297 00298 // GRAD and D1: Correct gradients and D1 in (F,P,D) format 00299 // 9 basis functions, each evaluated at 10 points, with two 00300 // components at each point. 00301 // that looks like 10 per to me. 00302 double basisGrads[] = { 00303 // 00304 -1.500000000000000, -1.500000000000000, 0.5000000000000000, 0, 0, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \ 00305 0, 0, 0, 0, 0, -0.5000000000000000, 0, 0, -0.08000000000000002, 0.1222222222222222, \ 00306 // 00307 2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, 0, 0, -1.500000000000000, \ 00308 0, 0, 0, 0.5000000000000000, 0, 0, 0, -0.5000000000000000, -0.3199999999999999, -0.9777777777777779, \ 00309 // 00310 -0.5000000000000000, 0, 1.500000000000000, -1.500000000000000, 0, 0.5000000000000000, 0, 0, 0.5000000000000000, 0, \ 00311 0, -0.5000000000000000, 0, 0, 0, 0, 0, 0, 0.3999999999999999, -0.2444444444444444, \ 00312 // 00313 0, 2.0, 0, 0, 0, 0, 0, -2.000000000000000, 0, 0, \ 00314 0.5000000000000000, 0, 0, 0, -1.50, 0, -0.50, 0, -0.1066666666666667, -0.1333333333333333, \ 00315 // 00316 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.0,\ 00317 -2.00, 0, 0, -2.0, 2.0, 0, 0, 0, -0.4266666666666667, 1.066666666666667 , \ 00318 // 00319 0, 0, 0, 2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, \ 00320 1.5, 0, 0, 0, -0.5, 0, 0.5000000000000000, 0, 0.5333333333333334, 0.2666666666666666 , \ 00321 // 00322 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, -1.500000000000000, 1.500000000000000, 0, 0, \ 00323 0, 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, 0.02000000000000000, 0.01111111111111112 , \ 00324 // 00325 0, 0, 0, 0, -2.0, 0, 2.0, 0, 0, -0.50, \ 00326 0, 0, 0, 1.5, 0, 0, 0, 0.5000000000000000, 0.07999999999999997, -0.08888888888888888, \ 00327 // 00328 0, 0, 0, -0.5000000000000000, 1.500000000000000, 1.500000000000000, -0.5000000000000000, 0, 0, 0, \ 00329 0, 0.5000000000000000, 0.5000000000000000, 0, 0, 0, 0, 0, -0.09999999999999998, -0.02222222222222221 \ 00330 // 00331 }; 00332 00333 // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D 00334 // 10 quad points and 3 values per point, so 00335 // each bf consists of 30 values. 00336 double basisD2[] = { 00337 1.0, 2.25, 1.0, 1.0, -0.75, 0, 0, 0.25, 0, 0, -0.75, 1.0, 1.0, 0.75, 0, 0, -0.25, 0, 0, -0.25, 0, 0, 0.75, 1.0, 0, 0.25, 0, 0.48, 0.1833333333333334, -0.1111111111111111, 00338 // 00339 -2.0, -3.0, 0, -2.0, 3.0, 0, 0, -1.0, 0, \ 00340 0, 1.0, 0, -2.0, 0, 1.0, 0, \ 00341 1.0, 0, 0, 0, 1.0, 0, -1.0, \ 00342 0, 0, 0, 1.0, -0.96, 0.7333333333333332, \ 00343 0.8888888888888890, \ 00344 // 00345 1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, 0.75, 1.0, 0, -0.25, 0, \ 00346 1.0, -0.75, 0, 0, -0.75, 1.0, 0, 0.25, 0, 0, 0.25, \ 00347 0, 0, -0.25, 0, 0.48, -0.9166666666666666, 0.2222222222222222, 00348 // 00349 0, -3.0, -2.0, 0, 1.0, 0, 0, -1.0, 0, 0, 3.0, \ 00350 -2.0, 0, -1.0, 0, 1.0, 0, 0, 0, 1.0, 0, 1.0, 0, -2.0, \ 00351 1.0, 0, 0, 0.6400000000000001, -0.2000000000000001, 0.2222222222222222, \ 00352 // 00353 0, 4.0, 0, 0, -4.0, 0, 0, 4.0, 0, 0, -4.0, 0, 0, 0, \ 00354 -2.0, -2.0, 0, 0, 0, 0, -2.0, -2.0, 0, 0, -2.0, 0, \ 00355 -2.0, -1.280000000000000, -0.7999999999999998, -1.777777777777778 , 00356 // 00357 0, -1.0, 0, 0, 3.0, -2.0, 0, -3.0, -2.0, 0, \ 00358 1.0, 0, 0, 1.0, 0, 1.0, 0, -2.0, 0, -1.0, 0, 1.0, 0, \ 00359 0, 1.0, 0, 0, 0.6400000000000001, 1.0, -0.4444444444444444, \ 00360 // 00361 0, 0.75, 1.0, 0, -0.25, 0, 1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, \ 00362 0.25, 0, 0, 0.25, 0, 1.0, -0.75, 0, 0, -0.75, 1.0, 0, \ 00363 -0.25, 0, -0.12, 0.01666666666666666, -0.1111111111111111, \ 00364 // 00365 0, -1.0, 0, 0, 1.0, 0, -2.0, -3.0, 0, -2.0, 3.0, 0, 0, 0, 1.0, 0, -1.0, \ 00366 0, -2.0, 0, 1.0, 0, 1.0, 0, \ 00367 0, 0, 1.0, 0.24, 0.06666666666666665, 0.8888888888888890, \ 00368 // 00369 0, 0.25, 0, 0, -0.75, 1.0, 1.0, 2.25, 1.0, 1.0, \ 00370 -0.75, 0, 0, -0.25, 0, 0, 0.75, 1.0, 1.0, \ 00371 0.75, 0, 0, -0.25, 0, 0, 0.25, 0, -0.12, -0.08333333333333331, 0.2222222222222222 \ 00372 }; 00373 00374 //D3: Correct multiset of second order partials in (F,P,Dk) format. D3 cardinality = 4 for 2D 00375 double basisD3[] = { 00376 0, -1.5, -1.5, 0, 0, -1.5, 0.5, 0, 0, 0.5, 00377 0.5, 0, 0, 0.5, -1.5, 0, 0, -1.5, -0.5, 0, 00378 0, -0.5, 0.5, 0, 0, 0.5, -0.5, 0, 0, -0.5, 00379 -1.5, 0, 0, -0.5, -0.5, 0, 0, -1.1, -0.1666666666666667, 0, 00380 // 00381 0, 3.0, 2.0, 0, 0, 3.0, -2.0, 0, 0, -1.0, 00382 -2.0, 0, 0, -1.0, 2.0, 0, 0, 3.0, 0, 0, 00383 0, 1.0, -2.0, 0, 0, -1.0, 0, 0, 0, 1.0, 00384 2.0, 0, 0, 1.0, 0, 0, 0, 2.2, -0.6666666666666665, 0, 00385 // 00386 0, -1.5, -0.5, 0, 0, -1.5, 1.5, 0, 0, 0.5, 00387 1.5, 0, 0, 0.5, -0.5, 0, 0, -1.5, 0.5, 0, 00388 0, -0.5, 1.5, 0, 0, 0.5, 0.5, 0, 0, -0.5, 00389 -0.5, 0, 0, -0.5, 0.5, 0, 0, -1.1, 0.8333333333333333, 0, 00390 // 00391 0, 2.0, 3.0, 0, 0, 2.0, -1.0, 0, 0, -2.0, 00392 -1.0, 0, 0, -2.0, 3.0, 0, 0, 2.0, 1.0, 0, 00393 0, 0, -1.0, 0, 0, -2.0, 1.0, 0, 0, 0, 00394 3.0, 0, 0, 0, 1.0, 0, 0, 1.2, 0.3333333333333334, 0, 00395 // 00396 0, -4.0, -4.0, 0, 0, -4.0, 4.0, 0, 0, 4.0, 00397 4.0, 0, 0, 4.0, -4.0, 0, 0, -4.0, 0, 0, 00398 0, 0, 4.0, 0, 0, 4.0, 0, 0, 0, 0, 00399 -4.0, 0, 0, 0, 0, 0, 0, -2.40, 1.333333333333333, 0, 00400 // 00401 0, 2.0, 1.0, 0, 0, 2.0, -3.0, 0, 0, -2.0, 00402 -3.0, 0, 0, -2.0, 1.0, 0, 0, 2.0, -1.0, 0, 00403 0, 0, -3.0, 0, 0, -2.0, -1.0, 0, 0, 0, 00404 1.0, 0, 0, 0, -1.0, 0, 0, 1.2, -1.666666666666667, 0 , 00405 // 00406 0, -0.5, -1.5, 0, 0, -0.5, 0.5, 0, 0, 1.5, 00407 0.5, 0, 0, 1.5, -1.5, 0, 0, -0.5, -0.5, 0, 00408 0, 0.5, 0.5, 0, 0, 1.5, -0.5, 0, 0, 0.5, 00409 -1.5, 0, 0, 0.5, -0.5, 0, 0, -0.09999999999999998, -0.1666666666666667, 0, 00410 // 00411 0, 1.0, 2.0, 0, 0, 1.0, -2.0, 0, 0, -3.0, 00412 -2.0, 0, 0, -3.0, 2.0, 0, 0, 1.0, 0, 0, 00413 0, -1.0, -2.0, 0, 0, -3.0, 0, 0, 0, -1.0, 00414 2.0, 0, 0, -1.0, 0, 0, 0, 0.2, -0.6666666666666665, 0, 00415 // 00416 0, -0.5, -0.5, 0, 0, -0.5, 1.5, 0, 0, 1.5, 00417 1.5, 0, 0, 1.5, -0.5, 0, 0, -0.5, 0.5, 0, 00418 0, 0.5, 1.5, 0, 0, 1.5, 0.5, 0, 0, 0.5, 00419 -0.5, 0, 0, 0.5, 0.5, 0, 0, -0.09999999999999998, 0.8333333333333333, 0 00420 }; 00421 //D4: Correct multiset of second order partials in (F,P,Dk) format. D4 cardinality = 5 for 2D 00422 double basisD4[] = { 00423 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00424 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00425 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00426 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00427 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00428 // 00429 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00430 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00431 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00432 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00433 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00434 // 00435 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00436 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00437 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00438 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00439 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00440 // 00441 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00442 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00443 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00444 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00445 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00446 // 00447 0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 00448 0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 00449 0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 00450 0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 00451 0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 00452 // 00453 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00454 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00455 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00456 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00457 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00458 // 00459 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00460 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00461 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00462 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00463 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00464 // 00465 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00466 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00467 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00468 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00469 0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 00470 // 00471 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00472 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00473 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00474 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 00475 0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0 00476 }; 00477 00478 try{ 00479 00480 // Dimensions for the output arrays: 00481 int numFields = quadBasis.getCardinality(); 00482 int numPoints = quadNodes.dimension(0); 00483 int spaceDim = quadBasis.getBaseCellTopology().getDimension(); 00484 00485 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00486 FieldContainer<double> vals; 00487 00488 // Check VALUE of basis functions: resize vals to rank-2 container: 00489 vals.resize(numFields, numPoints); 00490 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE); 00491 for (int i = 0; i < numFields; i++) { 00492 for (int j = 0; j < numPoints; j++) { 00493 00494 // Compute offset for (F,P) container 00495 int l = j + i * numPoints; 00496 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) { 00497 errorFlag++; 00498 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00499 00500 // Output the multi-index of the value where the error is: 00501 *outStream << " At multi-index { "; 00502 *outStream << i << " ";*outStream << j << " "; 00503 *outStream << "} computed value: " << vals(i,j) 00504 << " but reference value: " << basisValues[l] << "\n"; 00505 } 00506 } 00507 } 00508 00509 // Check GRAD of basis function: resize vals to rank-3 container 00510 vals.resize(numFields, numPoints, spaceDim); 00511 quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD); 00512 for (int i = 0; i < numFields; i++) { 00513 for (int j = 0; j < numPoints; j++) { 00514 for (int k = 0; k < spaceDim; k++) { 00515 00516 // basisGrads is (F,P,D), compute offset: 00517 int l = k + j * spaceDim + i * spaceDim * numPoints; 00518 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00519 errorFlag++; 00520 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00521 00522 // Output the multi-index of the value where the error is: 00523 *outStream << " At multi-index { "; 00524 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00525 *outStream << "} computed grad component: " << vals(i,j,k) 00526 << " but reference grad component: " << basisGrads[l] << "\n"; 00527 } 00528 } 00529 } 00530 } 00531 00532 00533 // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD) 00534 quadBasis.getValues(vals, quadNodes, OPERATOR_D1); 00535 for (int i = 0; i < numFields; i++) { 00536 for (int j = 0; j < numPoints; j++) { 00537 for (int k = 0; k < spaceDim; k++) { 00538 00539 // basisGrads is (F,P,D), compute offset: 00540 int l = k + j * spaceDim + i * spaceDim * numPoints; 00541 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) { 00542 errorFlag++; 00543 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00544 00545 // Output the multi-index of the value where the error is: 00546 *outStream << " At multi-index { "; 00547 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00548 *outStream << "} computed D1 component: " << vals(i,j,k) 00549 << " but reference D1 component: " << basisGrads[l] << "\n"; 00550 } 00551 } 00552 } 00553 } 00554 00555 00556 // Check CURL of basis function: resize vals just for illustration! 00557 vals.resize(numFields, numPoints, spaceDim); 00558 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL); 00559 for (int i = 0; i < numFields; i++) { 00560 for (int j = 0; j < numPoints; j++) { 00561 // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x) 00562 int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints; // position of y-derivative 00563 int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints; // position of x-derivative 00564 00565 double curl_value_0 = basisGrads[curl_0]; 00566 double curl_value_1 =-basisGrads[curl_1]; 00567 if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) { 00568 errorFlag++; 00569 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00570 // Output the multi-index of the value where the error is: 00571 *outStream << " At multi-index { "; 00572 *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " "; 00573 *outStream << "} computed curl component: " << vals(i,j,0) 00574 << " but reference curl component: " << curl_value_0 << "\n"; 00575 } 00576 if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) { 00577 errorFlag++; 00578 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00579 // Output the multi-index of the value where the error is: 00580 *outStream << " At multi-index { "; 00581 *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " "; 00582 *outStream << "} computed curl component: " << vals(i,j,1) 00583 << " but reference curl component: " << curl_value_1 << "\n"; 00584 } 00585 } 00586 } 00587 00588 // Check D2 of basis function 00589 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim); 00590 vals.resize(numFields, numPoints, D2cardinality); 00591 quadBasis.getValues(vals, quadNodes, OPERATOR_D2); 00592 for (int i = 0; i < numFields; i++) { 00593 for (int j = 0; j < numPoints; j++) { 00594 for (int k = 0; k < D2cardinality; k++) { 00595 00596 // basisD2 is (F,P,Dk), compute offset: 00597 int l = k + j * D2cardinality + i * D2cardinality * numPoints; 00598 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) { 00599 errorFlag++; 00600 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00601 00602 // Output the multi-index of the value where the error is: 00603 *outStream << " At multi-index { "; 00604 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00605 *outStream << "} computed D2 component: " << vals(i,j,k) 00606 << " but reference D2 component: " << basisD2[l] << "\n"; 00607 } 00608 } 00609 } 00610 } 00611 00612 00613 // Check D3 of basis function 00614 int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim); 00615 vals.resize(numFields, numPoints, D3cardinality); 00616 quadBasis.getValues(vals, quadNodes, OPERATOR_D3); 00617 for (int i = 0; i < numFields; i++) { 00618 for (int j = 0; j < numPoints; j++) { 00619 for (int k = 0; k < D3cardinality; k++) { 00620 00621 // basisD3 is (F,P,Dk), compute offset: 00622 int l = k + j * D3cardinality + i * D3cardinality * numPoints; 00623 if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) { 00624 errorFlag++; 00625 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00626 00627 // Output the multi-index of the value where the error is: 00628 *outStream << " At multi-index { "; 00629 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00630 *outStream << "} computed D3 component: " << vals(i,j,k) 00631 << " but reference D3 component: " << basisD2[l] << "\n"; 00632 } 00633 } 00634 } 00635 } 00636 00637 // Check D4 of basis function 00638 int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim); 00639 vals.resize(numFields, numPoints, D4cardinality); 00640 quadBasis.getValues(vals, quadNodes, OPERATOR_D4); 00641 for (int i = 0; i < numFields; i++) { 00642 for (int j = 0; j < numPoints; j++) { 00643 for (int k = 0; k < D4cardinality; k++) { 00644 00645 // basisD4 is (F,P,Dk), compute offset: 00646 int l = k + j * D4cardinality + i * D4cardinality * numPoints; 00647 if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) { 00648 errorFlag++; 00649 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00650 00651 // Output the multi-index of the value where the error is: 00652 *outStream << " At multi-index { "; 00653 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00654 *outStream << "} computed D4 component: " << vals(i,j,k) 00655 << " but reference D4 component: " << basisD4[l] << "\n"; 00656 } 00657 } 00658 } 00659 } 00660 00661 00662 // Check all higher derivatives - must be zero. 00663 for(EOperator op = OPERATOR_D5; op <= OPERATOR_D6; op++) { 00664 00665 // The last dimension is the number of kth derivatives and needs to be resized for every Dk 00666 int DkCardin = Intrepid::getDkCardinality(op, spaceDim); 00667 vals.resize(numFields, numPoints, DkCardin); 00668 00669 quadBasis.getValues(vals, quadNodes, op); 00670 for (int i = 0; i < vals.size(); i++) { 00671 if (std::abs(vals[i]) > INTREPID_TOL) { 00672 errorFlag++; 00673 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00674 00675 // Get the multi-index of the value where the error is and the operator order 00676 std::vector<int> myIndex; 00677 vals.getMultiIndex(myIndex,i); 00678 int ord = Intrepid::getOperatorOrder(op); 00679 *outStream << " At multi-index { "; 00680 for(int j = 0; j < vals.rank(); j++) { 00681 *outStream << myIndex[j] << " "; 00682 } 00683 *outStream << "} computed D"<< ord <<" component: " << vals[i] 00684 << " but reference D" << ord << " component: 0 \n"; 00685 } 00686 } 00687 } 00688 } 00689 00690 // Catch unexpected errors 00691 catch (std::logic_error err) { 00692 *outStream << err.what() << "\n\n"; 00693 errorFlag = -1000; 00694 }; 00695 00696 if (errorFlag != 0) 00697 std::cout << "End Result: TEST FAILED\n"; 00698 else 00699 std::cout << "End Result: TEST PASSED\n"; 00700 00701 // reset format state of std::cout 00702 std::cout.copyfmt(oldFormatState); 00703 00704 return errorFlag; 00705 }
1.7.6.1