Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HGRAD_TRI_C2_FEM/test_02.cpp
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 
00049 #include "Intrepid_FieldContainer.hpp"
00050 #include "Intrepid_HGRAD_TRI_C2_FEM.hpp"
00051 #include "Intrepid_DefaultCubatureFactory.hpp"
00052 #include "Intrepid_RealSpaceTools.hpp"
00053 #include "Intrepid_ArrayTools.hpp"
00054 #include "Intrepid_FunctionSpaceTools.hpp"
00055 #include "Intrepid_CellTools.hpp"
00056 #include "Teuchos_oblackholestream.hpp"
00057 #include "Teuchos_RCP.hpp"
00058 #include "Teuchos_GlobalMPISession.hpp"
00059 #include "Teuchos_SerialDenseMatrix.hpp"
00060 #include "Teuchos_SerialDenseVector.hpp"
00061 #include "Teuchos_LAPACK.hpp"
00062 
00063 using namespace std;
00064 using namespace Intrepid;
00065 
00066 void rhsFunc(FieldContainer<double> &, const FieldContainer<double> &, int, int);
00067 void neumann(FieldContainer<double>       & ,
00068              const FieldContainer<double> & ,
00069              const FieldContainer<double> & ,
00070              const shards::CellTopology   & ,
00071              int, int, int);
00072 void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int, int);
00073 
00075 void rhsFunc(FieldContainer<double> & result,
00076              const FieldContainer<double> & points,
00077              int xd,
00078              int yd) {
00079 
00080   int x = 0, y = 1;
00081 
00082   // second x-derivatives of u
00083   if (xd > 1) {
00084     for (int cell=0; cell<result.dimension(0); cell++) {
00085       for (int pt=0; pt<result.dimension(1); pt++) {
00086         result(cell,pt) = - xd*(xd-1)*std::pow(points(cell,pt,x), xd-2) * std::pow(points(cell,pt,y), yd);
00087       }
00088     }
00089   }
00090 
00091   // second y-derivatives of u
00092   if (yd > 1) {
00093     for (int cell=0; cell<result.dimension(0); cell++) {
00094       for (int pt=0; pt<result.dimension(1); pt++) {
00095         result(cell,pt) -=  yd*(yd-1)*std::pow(points(cell,pt,y), yd-2) * std::pow(points(cell,pt,x), xd);
00096       }
00097     }
00098   }
00099 
00100   // add u
00101   for (int cell=0; cell<result.dimension(0); cell++) {
00102     for (int pt=0; pt<result.dimension(1); pt++) {
00103       result(cell,pt) +=  std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
00104     }
00105   }
00106 
00107 }
00108 
00109 
00111 void neumann(FieldContainer<double>       & result,
00112              const FieldContainer<double> & points,
00113              const FieldContainer<double> & jacs,
00114              const shards::CellTopology   & parentCell,
00115              int sideOrdinal, int xd, int yd) {
00116 
00117   int x = 0, y = 1;
00118 
00119   int numCells  = result.dimension(0);
00120   int numPoints = result.dimension(1);
00121 
00122   FieldContainer<double> grad_u(numCells, numPoints, 2);
00123   FieldContainer<double> side_normals(numCells, numPoints, 2);
00124   FieldContainer<double> normal_lengths(numCells, numPoints);
00125 
00126   // first x-derivatives of u
00127   if (xd > 0) {
00128     for (int cell=0; cell<numCells; cell++) {
00129       for (int pt=0; pt<numPoints; pt++) {
00130         grad_u(cell,pt,x) = xd*std::pow(points(cell,pt,x), xd-1) * std::pow(points(cell,pt,y), yd);
00131       }
00132     }
00133   }
00134 
00135   // first y-derivatives of u
00136   if (yd > 0) {
00137     for (int cell=0; cell<numCells; cell++) {
00138       for (int pt=0; pt<numPoints; pt++) {
00139         grad_u(cell,pt,y) = yd*std::pow(points(cell,pt,y), yd-1) * std::pow(points(cell,pt,x), xd);
00140       }
00141     }
00142   }
00143   
00144   CellTools<double>::getPhysicalSideNormals(side_normals, jacs, sideOrdinal, parentCell);
00145 
00146   // scale normals
00147   RealSpaceTools<double>::vectorNorm(normal_lengths, side_normals, NORM_TWO);
00148   FunctionSpaceTools::scalarMultiplyDataData<double>(side_normals, normal_lengths, side_normals, true); 
00149 
00150   FunctionSpaceTools::dotMultiplyDataData<double>(result, grad_u, side_normals);
00151 
00152 }
00153 
00155 void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int xd, int yd) {
00156   int x = 0, y = 1;
00157   for (int cell=0; cell<result.dimension(0); cell++) {
00158     for (int pt=0; pt<result.dimension(1); pt++) {
00159       result(cell,pt) = std::pow(points(pt,x), xd)*std::pow(points(pt,y), yd);
00160     }
00161   }
00162 }
00163 
00164 
00165 
00166 
00167 int main(int argc, char *argv[]) {
00168 
00169   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00170 
00171   // This little trick lets us print to std::cout only if
00172   // a (dummy) command-line argument is provided.
00173   int iprint     = argc - 1;
00174   Teuchos::RCP<std::ostream> outStream;
00175   Teuchos::oblackholestream bhs; // outputs nothing
00176   if (iprint > 0)
00177     outStream = Teuchos::rcp(&std::cout, false);
00178   else
00179     outStream = Teuchos::rcp(&bhs, false);
00180 
00181   // Save the format state of the original std::cout.
00182   Teuchos::oblackholestream oldFormatState;
00183   oldFormatState.copyfmt(std::cout);
00184 
00185   *outStream \
00186     << "===============================================================================\n" \
00187     << "|                                                                             |\n" \
00188     << "|                    Unit Test (Basis_HGRAD_TRI_C2_FEM)                       |\n" \
00189     << "|                                                                             |\n" \
00190     << "|     1) Patch test involving mass and stiffness matrices,                    |\n" \
00191     << "|        for the Neumann problem on a triangular patch                        |\n" \
00192     << "|        Omega with boundary Gamma.                                           |\n" \
00193     << "|                                                                             |\n" \
00194     << "|        - div (grad u) + u = f  in Omega,  (grad u) . n = g  on Gamma        |\n" \
00195     << "|                                                                             |\n" \
00196     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00197     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00198     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00199     << "|                                                                             |\n" \
00200     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00201     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00202     << "|                                                                             |\n" \
00203     << "===============================================================================\n"\
00204     << "| TEST 1: Patch test                                                          |\n"\
00205     << "===============================================================================\n";
00206 
00207   
00208   int errorFlag = 0;
00209 
00210   outStream -> precision(16);
00211 
00212 
00213   try {
00214 
00215     int max_order = 2;                                                               // max total order of polynomial solution
00216     DefaultCubatureFactory<double>  cubFactory;                                      // create cubature factory
00217     shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >());  // create parent cell topology
00218     shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >());      // create relevant subcell (side) topology
00219     int cellDim = cell.getDimension();
00220     int sideDim = side.getDimension();
00221 
00222     // Define array containing points at which the solution is evaluated, in reference cell.
00223     int numIntervals = 10;
00224     int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2;
00225     FieldContainer<double> interp_points_ref(numInterpPoints, 2);
00226     int counter = 0;
00227     for (int j=0; j<=numIntervals; j++) {
00228       for (int i=0; i<=numIntervals; i++) {
00229         if (i <= numIntervals-j) {
00230           interp_points_ref(counter,0) = i*(1.0/numIntervals);
00231           interp_points_ref(counter,1) = j*(1.0/numIntervals);
00232           counter++;
00233         }
00234       }
00235     }
00236 
00237     /* Parent cell definition. */
00238     FieldContainer<double> cell_nodes(1, 3, cellDim);
00239     // Perturbed reference triangle.
00240     cell_nodes(0, 0, 0) = 0.1;
00241     cell_nodes(0, 0, 1) = -0.1;
00242     cell_nodes(0, 1, 0) = 1.1;
00243     cell_nodes(0, 1, 1) = -0.1;
00244     cell_nodes(0, 2, 0) = 0.1;
00245     cell_nodes(0, 2, 1) = 0.9;
00246     // Reference triangle. 
00247     /*cell_nodes(0, 0, 0) = 0.0;
00248     cell_nodes(0, 0, 1) = 0.0;
00249     cell_nodes(0, 1, 0) = 1.0;
00250     cell_nodes(0, 1, 1) = 0.0;
00251     cell_nodes(0, 2, 0) = 0.0;
00252     cell_nodes(0, 2, 1) = 1.0;*/
00253 
00254     FieldContainer<double> interp_points(1, numInterpPoints, cellDim);
00255     CellTools<double>::mapToPhysicalFrame(interp_points, interp_points_ref, cell_nodes, cell);
00256     interp_points.resize(numInterpPoints, cellDim);
00257 
00258     for (int x_order=0; x_order <= max_order; x_order++) {
00259       for (int y_order=0; y_order <= max_order-x_order; y_order++) {
00260 
00261         // evaluate exact solution
00262         FieldContainer<double> exact_solution(1, numInterpPoints);
00263         u_exact(exact_solution, interp_points, x_order, y_order);
00264 
00265         int basis_order = 2;
00266 
00267         // set test tolerance
00268         double zero = basis_order*basis_order*100*INTREPID_TOL;
00269 
00270         //create basis
00271         Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
00272           Teuchos::rcp(new Basis_HGRAD_TRI_C2_FEM<double,FieldContainer<double> >() );
00273         int numFields = basis->getCardinality();
00274 
00275         // create cubatures
00276         Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*basis_order);
00277         Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*basis_order);
00278         int numCubPointsCell = cellCub->getNumPoints();
00279         int numCubPointsSide = sideCub->getNumPoints();
00280 
00281         /* Computational arrays. */
00282         /* Section 1: Related to parent cell integration. */
00283         FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
00284         FieldContainer<double> cub_points_cell_physical(1, numCubPointsCell, cellDim);
00285         FieldContainer<double> cub_weights_cell(numCubPointsCell);
00286         FieldContainer<double> jacobian_cell(1, numCubPointsCell, cellDim, cellDim);
00287         FieldContainer<double> jacobian_inv_cell(1, numCubPointsCell, cellDim, cellDim);
00288         FieldContainer<double> jacobian_det_cell(1, numCubPointsCell);
00289         FieldContainer<double> weighted_measure_cell(1, numCubPointsCell);
00290 
00291         FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell);
00292         FieldContainer<double> transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
00293         FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
00294         FieldContainer<double> grad_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim);
00295         FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
00296         FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
00297         FieldContainer<double> fe_matrix(1, numFields, numFields);
00298 
00299         FieldContainer<double> rhs_at_cub_points_cell_physical(1, numCubPointsCell);
00300         FieldContainer<double> rhs_and_soln_vector(1, numFields);
00301 
00302         /* Section 2: Related to subcell (side) integration. */
00303         unsigned numSides = 3;
00304         FieldContainer<double> cub_points_side(numCubPointsSide, sideDim);
00305         FieldContainer<double> cub_weights_side(numCubPointsSide);
00306         FieldContainer<double> cub_points_side_refcell(numCubPointsSide, cellDim);
00307         FieldContainer<double> cub_points_side_physical(1, numCubPointsSide, cellDim);
00308         FieldContainer<double> jacobian_side_refcell(1, numCubPointsSide, cellDim, cellDim);
00309         FieldContainer<double> jacobian_det_side_refcell(1, numCubPointsSide);
00310         FieldContainer<double> weighted_measure_side_refcell(1, numCubPointsSide);
00311 
00312         FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields, numCubPointsSide);
00313         FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
00314         FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
00315         FieldContainer<double> neumann_data_at_cub_points_side_physical(1, numCubPointsSide);
00316         FieldContainer<double> neumann_fields_per_side(1, numFields);
00317 
00318         /* Section 3: Related to global interpolant. */
00319         FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
00320         FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
00321         FieldContainer<double> interpolant(1, numInterpPoints);
00322 
00323         FieldContainer<int> ipiv(numFields);
00324 
00325 
00326 
00327         /******************* START COMPUTATION ***********************/
00328 
00329         // get cubature points and weights
00330         cellCub->getCubature(cub_points_cell, cub_weights_cell);
00331 
00332         // compute geometric cell information
00333         CellTools<double>::setJacobian(jacobian_cell, cub_points_cell, cell_nodes, cell);
00334         CellTools<double>::setJacobianInv(jacobian_inv_cell, jacobian_cell);
00335         CellTools<double>::setJacobianDet(jacobian_det_cell, jacobian_cell);
00336 
00337         // compute weighted measure
00338         FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
00339 
00341         // Computing mass matrices:
00342         // tabulate values of basis functions at (reference) cubature points
00343         basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
00344 
00345         // transform values of basis functions
00346         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
00347                                                         value_of_basis_at_cub_points_cell);
00348 
00349         // multiply with weighted measure
00350         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
00351                                                     weighted_measure_cell,
00352                                                     transformed_value_of_basis_at_cub_points_cell);
00353 
00354         // compute mass matrices
00355         FunctionSpaceTools::integrate<double>(fe_matrix,
00356                                               transformed_value_of_basis_at_cub_points_cell,
00357                                               weighted_transformed_value_of_basis_at_cub_points_cell,
00358                                               COMP_BLAS);
00360 
00362         // Computing stiffness matrices:
00363         // tabulate gradients of basis functions at (reference) cubature points
00364         basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
00365 
00366         // transform gradients of basis functions
00367         FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
00368                                                        jacobian_inv_cell,
00369                                                        grad_of_basis_at_cub_points_cell);
00370 
00371         // multiply with weighted measure
00372         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
00373                                                     weighted_measure_cell,
00374                                                     transformed_grad_of_basis_at_cub_points_cell);
00375 
00376         // compute stiffness matrices and sum into fe_matrix
00377         FunctionSpaceTools::integrate<double>(fe_matrix,
00378                                               transformed_grad_of_basis_at_cub_points_cell,
00379                                               weighted_transformed_grad_of_basis_at_cub_points_cell,
00380                                               COMP_BLAS,
00381                                               true);
00383 
00385         // Computing RHS contributions:
00386         // map cell (reference) cubature points to physical space
00387         CellTools<double>::mapToPhysicalFrame(cub_points_cell_physical, cub_points_cell, cell_nodes, cell);
00388 
00389         // evaluate rhs function
00390         rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order);
00391 
00392         // compute rhs
00393         FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
00394                                               rhs_at_cub_points_cell_physical,
00395                                               weighted_transformed_value_of_basis_at_cub_points_cell,
00396                                               COMP_BLAS);
00397 
00398         // compute neumann b.c. contributions and adjust rhs
00399         sideCub->getCubature(cub_points_side, cub_weights_side);
00400         for (unsigned i=0; i<numSides; i++) {
00401           // compute geometric cell information
00402           CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
00403           CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes, cell);
00404           CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);
00405 
00406           // compute weighted edge measure
00407           FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
00408                                                          jacobian_side_refcell,
00409                                                          cub_weights_side,
00410                                                          i,
00411                                                          cell);
00412 
00413           // tabulate values of basis functions at side cubature points, in the reference parent cell domain
00414           basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
00415           // transform
00416           FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
00417                                                           value_of_basis_at_cub_points_side_refcell);
00418 
00419           // multiply with weighted measure
00420           FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
00421                                                       weighted_measure_side_refcell,
00422                                                       transformed_value_of_basis_at_cub_points_side_refcell);
00423 
00424           // compute Neumann data
00425           // map side cubature points in reference parent cell domain to physical space
00426           CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes, cell);
00427           // now compute data
00428           neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
00429                   cell, (int)i, x_order, y_order);
00430 
00431           FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
00432                                                 neumann_data_at_cub_points_side_physical,
00433                                                 weighted_transformed_value_of_basis_at_cub_points_side_refcell,
00434                                                 COMP_BLAS);
00435 
00436           // adjust RHS
00437           RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
00438         }
00440 
00442         // Solution of linear system:
00443         int info = 0;
00444         Teuchos::LAPACK<int, double> solver;
00445         solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
00447 
00449         // Building interpolant:
00450         // evaluate basis at interpolation points
00451         basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
00452         // transform values of basis functions
00453         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
00454                                                         value_of_basis_at_interp_points);
00455         FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
00457 
00458         /******************* END COMPUTATION ***********************/
00459     
00460         RealSpaceTools<double>::subtract(interpolant, exact_solution);
00461 
00462         *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
00463                    << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
00464                    << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
00465                       RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";
00466 
00467         if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
00468             RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
00469           *outStream << "\n\nPatch test failed for solution polynomial order ("
00470                      << x_order << ", " << y_order << ") and basis order " << basis_order << "\n\n";
00471           errorFlag++;
00472         }
00473       } // end for y_order
00474     } // end for x_order
00475 
00476   }
00477   // Catch unexpected errors
00478   catch (std::logic_error err) {
00479     *outStream << err.what() << "\n\n";
00480     errorFlag = -1000;
00481   };
00482 
00483   if (errorFlag != 0)
00484     std::cout << "End Result: TEST FAILED\n";
00485   else
00486     std::cout << "End Result: TEST PASSED\n";
00487 
00488   // reset format state of std::cout
00489   std::cout.copyfmt(oldFormatState);
00490 
00491   return errorFlag;
00492 }