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