|
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_TET_COMP12_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_TET_COMP12_FEM) |\n" \ 00093 << "| |\n" \ 00094 << "| 1) Evaluation of Basis Function Values |\n" \ 00095 << "| |\n" \ 00096 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00097 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00098 << "| Kara Peterson (kjpeter@sandia.gov), |\n" \ 00099 << "| Jake Ostien (jtostie@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, values |\n"\ 00106 << "===============================================================================\n"; 00107 00108 // Define basis and error flag 00109 Basis_HGRAD_TET_COMP12_FEM<double, FieldContainer<double> > tetBasis; 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 10 vertices of the reference TET 00117 FieldContainer<double> tetNodes(10, 3); 00118 tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0; 00119 tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0; 00120 tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0; 00121 tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0; 00122 tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0; 00123 tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0; 00124 tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0; 00125 tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5; 00126 tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5; 00127 tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5; 00128 // Define array containing 5 integration points 00129 FieldContainer<double> tetPoints(9, 3); 00130 // from the 5 point integration 00131 tetPoints(0,0) = 0.25; tetPoints(0,1) = 0.25; tetPoints(0,2) = 0.25; 00132 tetPoints(1,0) = 0.5; tetPoints(1,1) = (1./6.); tetPoints(1,2) = (1./6.); 00133 tetPoints(2,0) = (1./6.); tetPoints(2,1) = 0.5; tetPoints(2,2) = (1./6.); 00134 tetPoints(3,0) = (1./6.); tetPoints(3,1) = (1./6.); tetPoints(3,2) = 0.5; 00135 tetPoints(4,0) = (1./6.); tetPoints(4,1) = (1./6.); tetPoints(4,2) = (1./6.); 00136 // from the 4 point integration 00137 tetPoints(5,0) = 0.1381966011250105151795413165634361882280; 00138 tetPoints(5,1) = 0.1381966011250105151795413165634361882280; 00139 tetPoints(5,2) = 0.1381966011250105151795413165634361882280; 00140 00141 tetPoints(6,0) = 0.5854101966249684544613760503096914353161; 00142 tetPoints(6,1) = 0.1381966011250105151795413165634361882280; 00143 tetPoints(6,2) = 0.1381966011250105151795413165634361882280; 00144 00145 tetPoints(7,0) = 0.1381966011250105151795413165634361882280; 00146 tetPoints(7,1) = 0.5854101966249684544613760503096914353161; 00147 tetPoints(7,2) = 0.1381966011250105151795413165634361882280; 00148 00149 tetPoints(8,0) = 0.1381966011250105151795413165634361882280; 00150 tetPoints(8,1) = 0.1381966011250105151795413165634361882280; 00151 tetPoints(8,2) = 0.5854101966249684544613760503096914353161; 00152 00153 // output precision 00154 outStream -> precision(20); 00155 00156 // VALUE: Each row gives the 10 correct basis set values at an evaluation point 00157 double nodalBasisValues[] = { 00158 // first 4 vertices 00159 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00160 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00161 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00162 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00163 // second 6 vertices 00164 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00165 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 00166 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 00167 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00168 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00169 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 00170 }; 00171 double pointBasisValues[] = { 00172 // pt 0 {0.25, 0.25, 0.25} 00173 0.0, 0.0, 0.0, 0.0, 1./6., 1./6., 1./6., 1./6., 1./6., 1./6., 00174 // pt 1 {0.5, 1/6, 1/6} 00175 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3., 0.0, 00176 // pt 2 {1/6, 0.5, 0.1/6} 00177 0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3., 00178 // pt 3 {1/6, 1/6, 0.5} 00179 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 1./3., 00180 // pt 4 {1/6, 1/6, 1/6} 00181 0.0, 0.0, 0.0, 0.0, 1./3., 0.0, 1./3., 1./3., 0.0, 0.0, 00182 // pt 5 00183 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0, 00184 // pt 6 00185 0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.0, 00186 // pt 7 00187 0.0, 0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0, 0.276393202250021030359082633126872376456, 00188 // pt 8 00189 0.0, 0.0, 0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 00190 }; 00191 00192 // GRAD and D1: each row gives the 3x10 correct values of the gradients of the 10 basis functions 00193 double pointBasisGrads[] = { 00194 // point 0 00195 -1./4., -1./4., -1./4., \ 00196 1./4., 0.0, 0.0, \ 00197 0.0, 1./4., 0.0, \ 00198 0.0, 0.0, 1./4., \ 00199 0.0, -3./4., -3./4., \ 00200 3./4., 3./4., 0.0, \ 00201 -3./4., 0.0, -3./4., \ 00202 -3./4., -3./4., 0.0, \ 00203 3./4., 0.0, 3./4., \ 00204 0.0, 3./4., 3./4., \ 00205 00206 // point 1 00207 -1./24., -1./24., -1./24., \ 00208 7./8., 0.0, 0.0, \ 00209 0.0, 1./24., 0.0, \ 00210 0.0, 0.0, 1./24., \ 00211 -35./36., -19./12., -19./12., \ 00212 11./18., 19./12., 0.0, \ 00213 -17./36., 0.0, -1./3., \ 00214 -17./36., -1./3., 0.0, \ 00215 11./18., 0.0, 19./12., \ 00216 -5./36., 1./3., 1./3., \ 00217 00218 // point 2 00219 -1./24., -1./24., -1./24., \ 00220 1./24., 0.0, 0.0, \ 00221 0.0, 7./8., 0.0, \ 00222 0.0, 0.0, 1./24., \ 00223 0.0, -17./36., -1./3., \ 00224 19./12., 11./18., 0.0, \ 00225 -19./12., -35./36., -19./12., \ 00226 -1./3., -17./36., 0.0, \ 00227 1./3., -5./36., 1./3., \ 00228 0.0, 11./18., 19./12., \ 00229 00230 // point 3 00231 -1./24., -1./24., -1./24., \ 00232 1./24., 0.0, 0.0, \ 00233 0.0, 1./24., 0.0, \ 00234 0.0, 0.0, 7./8., \ 00235 0.0, -1./3., -17./36., \ 00236 1./3., 1./3., -5./36., \ 00237 -1./3., 0.0, -17./36., \ 00238 -19./12., -19./12., -35./36., \ 00239 19./12., 0.0, 11./18., \ 00240 0.0, 19./12., 11./18., \ 00241 00242 // point 4 00243 -7./8., -7./8., -7./8., \ 00244 1./24., 0.0, 0.0, \ 00245 0.0, 1./24., 0.0, \ 00246 0.0, 0.0, 1./24., \ 00247 35./36., -11./18., -11./18., \ 00248 17./36., 17./36., 5./36., \ 00249 -11./18., 35./36., -11./18., \ 00250 -11./18., -11./18., 35./36., \ 00251 17./36., 5./36., 17./36., \ 00252 5./36., 17./36., 17./36., \ 00253 00254 // point 5 00255 -1.088525491562421136153440125774228588290, -1.088525491562421136153440125774228588290, -1.088525491562421136153440125774228588290, \ 00256 -0.029508497187473712051146708591409529430, 0.0, 0.0, \ 00257 0.0, -0.029508497187473712051146708591409529430, 0.0, \ 00258 0.0, 0.0, -0.029508497187473712051146708591409529430, \ 00259 1.30437298687487732290535130675991113734, -0.563661001875017525299235527605726980380, -0.563661001875017525299235527605726980380, \ 00260 0.377322003750035050598471055211453960760, 0.377322003750035050598471055211453960760, 0.186338998124982474700764472394273019620, \ 00261 -0.563661001875017525299235527605726980380, 1.30437298687487732290535130675991113734, -0.563661001875017525299235527605726980380, \ 00262 -0.563661001875017525299235527605726980380, -0.563661001875017525299235527605726980380, 1.30437298687487732290535130675991113734, \ 00263 0.377322003750035050598471055211453960760, 0.186338998124982474700764472394273019620, 0.377322003750035050598471055211453960760, \ 00264 0.186338998124982474700764472394273019620, 0.377322003750035050598471055211453960760, 0.377322003750035050598471055211453960760, \ 00265 00266 // point 6 00267 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \ 00268 1.088525491562421136153440125774228588290, 0.0, 0.0, \ 00269 0.0, -0.029508497187473712051146708591409529430, 0.0, \ 00270 0.0, 0.0, -0.029508497187473712051146708591409529430, \ 00271 -1.30437298687487732290535130675991113734, -1.868033988749894848204586834365638117720, -1.868033988749894848204586834365638117720, \ 00272 0.563661001875017525299235527605726980380, 1.868033988749894848204586834365638117720, 0.0, \ 00273 -0.377322003750035050598471055211453960760, 0.0, -0.190983005625052575897706582817180941140, \ 00274 -0.377322003750035050598471055211453960760, -0.190983005625052575897706582817180941140, 0.0, \ 00275 0.563661001875017525299235527605726980380, 0.0, 1.868033988749894848204586834365638117720, \ 00276 -0.186338998124982474700764472394273019620, 0.190983005625052575897706582817180941140, 0.19098300562505257589770658281718094114, \ 00277 00278 // point 7 00279 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \ 00280 -0.029508497187473712051146708591409529430, 0.0, 0.0, \ 00281 0.0, 1.088525491562421136153440125774228588290, 0.0, \ 00282 0.0, 0.0, -0.029508497187473712051146708591409529430, \ 00283 0.0, -0.377322003750035050598471055211453960760, -0.190983005625052575897706582817180941140, \ 00284 1.868033988749894848204586834365638117720, 0.563661001875017525299235527605726980380, 0.0, \ 00285 -1.868033988749894848204586834365638117720, -1.30437298687487732290535130675991113734, -1.868033988749894848204586834365638117720, \ 00286 -0.190983005625052575897706582817180941140, -0.377322003750035050598471055211453960760, 0.0, \ 00287 0.190983005625052575897706582817180941140, -0.186338998124982474700764472394273019620, 0.190983005625052575897706582817180941140, \ 00288 0.0, 0.563661001875017525299235527605726980380, 1.868033988749894848204586834365638117720, \ 00289 00290 // point 8 00291 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \ 00292 -0.029508497187473712051146708591409529430, 0.0, 0.0, \ 00293 0.0, -0.029508497187473712051146708591409529430, 0.0, \ 00294 0.0, 0.0, 1.088525491562421136153440125774228588290, \ 00295 0.0, -0.190983005625052575897706582817180941140, -0.377322003750035050598471055211453960760, \ 00296 0.190983005625052575897706582817180941140, 0.190983005625052575897706582817180941140, -0.186338998124982474700764472394273019620, \ 00297 -0.190983005625052575897706582817180941140, 0.0, -0.377322003750035050598471055211453960760, \ 00298 -1.868033988749894848204586834365638117720, -1.868033988749894848204586834365638117720, -1.30437298687487732290535130675991113734, 00299 1.868033988749894848204586834365638117720, 0.0, 0.563661001875017525299235527605726980380, \ 00300 0.0, 1.868033988749894848204586834365638117720, 0.563661001875017525299235527605726980380, \ 00301 }; 00302 00303 try{ 00304 00305 // Dimensions for the output arrays: 00306 int numFields = tetBasis.getCardinality(); 00307 int numNodes = tetNodes.dimension(0); 00308 int spaceDim = tetBasis.getBaseCellTopology().getDimension(); 00309 00310 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00311 FieldContainer<double> vals; 00312 00313 // Check VALUE of basis functions at nodes: resize vals to rank-2 container:\n"; 00314 *outStream << " check VALUE of basis functions at nodes\n"; 00315 vals.resize(numFields, numNodes); 00316 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE); 00317 00318 for (int i = 0; i < numFields; i++) { 00319 for (int j = 0; j < numNodes; j++) { 00320 int l = i + j * numFields; 00321 if (std::abs(vals(i,j) - nodalBasisValues[l]) > INTREPID_TOL) { 00322 errorFlag++; 00323 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00324 00325 // Output the multi-index of the value where the error is: 00326 *outStream << " At multi-index { "; 00327 *outStream << i << " ";*outStream << j << " "; 00328 *outStream << "} computed value: " << vals(i,j) 00329 << " but reference value: " << nodalBasisValues[l] << "\n"; 00330 } 00331 } 00332 } 00333 00334 // Check VALUE of basis functions at points: resize vals to rank-2 container:\n"; 00335 *outStream << " check VALUE of basis functions at points\n"; 00336 int numPoints = tetPoints.dimension(0); 00337 vals.resize(numFields, numPoints); 00338 vals.initialize(0.0); 00339 tetBasis.getValues(vals, tetPoints, OPERATOR_VALUE); 00340 00341 for (int i = 0; i < numFields; i++) { 00342 for (int j = 0; j < numPoints; j++) { 00343 int l = i + j * numFields; 00344 if (std::abs(vals(i,j) - pointBasisValues[l]) > INTREPID_TOL) { 00345 errorFlag++; 00346 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00347 00348 // Output the multi-index of the value where the error is: 00349 *outStream << " At multi-index { "; 00350 *outStream << i << " ";*outStream << j << " "; 00351 *outStream << "} computed value: " << vals(i,j) 00352 << " but reference value: " << pointBasisValues[l] << "\n"; 00353 } 00354 } 00355 } 00356 00357 // Check GRAD of basis functions at points: resize vals to rank-3 container:\n"; 00358 numPoints = tetPoints.dimension(0); 00359 vals.resize(numFields, numPoints, spaceDim); 00360 vals.initialize(0.0); 00361 tetBasis.getValues(vals, tetPoints, OPERATOR_GRAD); 00362 00363 for (int i = 0; i < numFields; i++) { 00364 for (int j = 0; j < numPoints; j++) { 00365 for (int k = 0; k < spaceDim; k++) { 00366 int l = k + i * spaceDim + j * spaceDim * numFields; 00367 if (std::abs(vals(i,j,k) - pointBasisGrads[l]) > INTREPID_TOL) { 00368 errorFlag++; 00369 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00370 00371 // Output the multi-index of the value where the error is: 00372 *outStream << " At multi-index { "; 00373 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00374 *outStream << "} computed grad component: " << vals(i,j,k) 00375 << " but reference grad component: " << pointBasisGrads[l] << "\n"; 00376 } 00377 } 00378 } 00379 } 00380 } 00381 // Catch unexpected errors 00382 catch (std::logic_error err) { 00383 *outStream << err.what() << "\n\n"; 00384 errorFlag = -1000; 00385 }; 00386 00387 00388 if (errorFlag != 0) 00389 std::cout << "End Result: TEST FAILED\n"; 00390 else 00391 std::cout << "End Result: TEST PASSED\n"; 00392 00393 // reset format state of std::cout 00394 std::cout.copyfmt(oldFormatState); 00395 00396 return errorFlag; 00397 }
1.7.6.1