|
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(5, 3); 00130 tetPoints(0,0) = 0.25; tetPoints(0,1) = 0.25; tetPoints(0,2) = 0.25; 00131 tetPoints(1,0) = 0.5; tetPoints(1,1) = (1./6.); tetPoints(1,2) = (1./6.); 00132 tetPoints(2,0) = (1./6.); tetPoints(2,1) = 0.5; tetPoints(2,2) = (1./6.); 00133 tetPoints(3,0) = (1./6.); tetPoints(3,1) = (1./6.); tetPoints(3,2) = 0.5; 00134 tetPoints(4,0) = (1./6.); tetPoints(4,1) = (1./6.); tetPoints(4,2) = (1./6.); 00135 00136 // output precision 00137 outStream -> precision(20); 00138 00139 // VALUE: Each row gives the 10 correct basis set values at an evaluation point 00140 double nodalBasisValues[] = { 00141 // first 4 vertices 00142 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00143 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00144 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00145 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00146 // second 6 vertices 00147 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00148 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 00149 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 00150 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00151 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00152 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 00153 }; 00154 double pointBasisValues[] = { 00155 // pt 0 {0.25, 0.25, 0.25} 00156 0.0, 0.0, 0.0, 0.0, 1./6., 1./6., 1./6., 1./6., 1./6., 1./6., 00157 // pt 1 {0.5, 1/6, 1/6} 00158 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3., 0.0, 00159 // pt 2 {1/6, 0.5, 0.1/6} 00160 0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3., 00161 // pt 3 {1/6, 1/6, 0.5} 00162 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 1./3., 00163 // pt 4 {1/6, 1/6, 1/6} 00164 0.0, 0.0, 0.0, 0.0, 1./3., 0.0, 1./3., 1./3., 0.0, 0.0 00165 }; 00166 00167 // GRAD and D1: each row gives the 3x10 correct values of the gradients of the 10 basis functions 00168 double pointBasisGrads[] = { 00169 // point 0 00170 -0.25, -0.25, -0.25, 0.25, 0.00, 0.00, 0.00, 0.25, 0.0, 0.00, 0.0, 0.25, 0.0, -0.75, -0.75, \ 00171 0.75, 0.75, 0.00, -0.75, 0.00, -0.75, -0.75, -0.75, 0.0, 0.75, 0.0, 0.75, 0.0, 0.75, 0.75, 00172 00173 // point 1 00174 0.15, 0.15, 0.15, 1.45, 0.00, 0.00, 0.00, -0.15, 0.0, 0.00, 0.0, -0.15, -104./60., -2.15, -2.15, \ 00175 5./12., 2.15, 0.00, -17./60., 0.00, -0.15, -17./60., -0.15, 0.0, 5./12., 0.0, 2.15, -2./15., 0.15, 0.15, 00176 00177 // point 2 00178 0.15, 0.15, 0.15, -0.15, 0.0, 0.0, 0.0, 1.45, 0.0, 0.0, 0.0, -0.15, 0.0, -17./60., -0.15, \ 00179 2.15, 5./12., 0.0, -2.15, -104./60., -2.15, -0.15, -17./60., 0.0, 0.15, -2./15., 0.15, 0.0, 5./12., 2.15, 00180 00181 // point 3 00182 0.15, 0.15, 0.15, -0.15, 0.0, 0.0, 0.0, -0.15, 0.0, 0.0, 0.0, 1.45, 0.0, -0.15, -17./60., \ 00183 0.15, 0.15, -2./15., -0.15, 0.0, -17./60., -2.15, -2.15, -104./60., 2.15, 0.0, 5./12., 0.0, 2.15, 5./12., 00184 00185 // point 4 00186 -1.45, -1.45, -1.45, -0.15, 0.0, 0.0, 0.0, -0.15, 0.0, 0.0, 0.0, -0.15, 104./60., -5./12., -5./12., \ 00187 17./60., 17./60., 2./15., -5./12., 104./60., -5./12., -5./12., -5./12., 104./60., 17./60., 2./15., 17./60., 2./15., 17./60., 17./60. 00188 }; 00189 00190 try{ 00191 00192 // Dimensions for the output arrays: 00193 int numFields = tetBasis.getCardinality(); 00194 int numNodes = tetNodes.dimension(0); 00195 int spaceDim = tetBasis.getBaseCellTopology().getDimension(); 00196 00197 // Generic array for values, grads, curls, etc. that will be properly sized before each call 00198 FieldContainer<double> vals; 00199 00200 // Check VALUE of basis functions at nodes: resize vals to rank-2 container:\n"; 00201 vals.resize(numFields, numNodes); 00202 tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE); 00203 00204 for (int i = 0; i < numFields; i++) { 00205 for (int j = 0; j < numNodes; j++) { 00206 int l = i + j * numFields; 00207 if (std::abs(vals(i,j) - nodalBasisValues[l]) > INTREPID_TOL) { 00208 errorFlag++; 00209 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00210 00211 // Output the multi-index of the value where the error is: 00212 *outStream << " At multi-index { "; 00213 *outStream << i << " ";*outStream << j << " "; 00214 *outStream << "} computed value: " << vals(i,j) 00215 << " but reference value: " << nodalBasisValues[l] << "\n"; 00216 } 00217 } 00218 } 00219 00220 // Check VALUE of basis functions at points: resize vals to rank-2 container:\n"; 00221 int numPoints = tetPoints.dimension(0); 00222 vals.resize(numFields, numPoints); 00223 vals.initialize(0.0); 00224 tetBasis.getValues(vals, tetPoints, OPERATOR_VALUE); 00225 00226 for (int i = 0; i < numFields; i++) { 00227 for (int j = 0; j < numPoints; j++) { 00228 int l = i + j * numFields; 00229 if (std::abs(vals(i,j) - pointBasisValues[l]) > INTREPID_TOL) { 00230 errorFlag++; 00231 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00232 00233 // Output the multi-index of the value where the error is: 00234 *outStream << " At multi-index { "; 00235 *outStream << i << " ";*outStream << j << " "; 00236 *outStream << "} computed value: " << vals(i,j) 00237 << " but reference value: " << pointBasisValues[l] << "\n"; 00238 } 00239 } 00240 } 00241 00242 // Check GRAD of basis functions at points: resize vals to rank-3 container:\n"; 00243 numPoints = tetPoints.dimension(0); 00244 vals.resize(numFields, numPoints, spaceDim); 00245 vals.initialize(0.0); 00246 tetBasis.getValues(vals, tetPoints, OPERATOR_GRAD); 00247 00248 for (int i = 0; i < numFields; i++) { 00249 for (int j = 0; j < numPoints; j++) { 00250 for (int k = 0; k < spaceDim; k++) { 00251 int l = k + i * spaceDim + j * spaceDim * numFields; 00252 if (std::abs(vals(i,j,k) - pointBasisGrads[l]) > INTREPID_TOL) { 00253 errorFlag++; 00254 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00255 00256 // Output the multi-index of the value where the error is: 00257 *outStream << " At multi-index { "; 00258 *outStream << i << " ";*outStream << j << " ";*outStream << k << " "; 00259 *outStream << "} computed grad component: " << vals(i,j,k) 00260 << " but reference grad component: " << pointBasisGrads[l] << "\n"; 00261 } 00262 } 00263 } 00264 } 00265 } 00266 // Catch unexpected errors 00267 catch (std::logic_error err) { 00268 *outStream << err.what() << "\n\n"; 00269 errorFlag = -1000; 00270 }; 00271 00272 00273 if (errorFlag != 0) 00274 std::cout << "End Result: TEST FAILED\n"; 00275 else 00276 std::cout << "End Result: TEST PASSED\n"; 00277 00278 // reset format state of std::cout 00279 std::cout.copyfmt(oldFormatState); 00280 00281 return errorFlag; 00282 }
1.7.6.1