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