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