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