Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HCURL_TET_In_FEM/test_03.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_HCURL_TET_In_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, int );
00067 void bcFunc( FieldContainer<double> &, const FieldContainer<double> & ,
00068              const shards::CellTopology   & ,
00069              int , int , int , int , int );
00070 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int );
00071 
00072 void u_exact( FieldContainer<double> &result, 
00073               const FieldContainer<double> &points,
00074               int comp, 
00075               int xd,
00076               int yd,
00077               int zd)
00078 {
00079   for (int cell=0;cell<result.dimension(0);cell++){
00080     for (int pt=0;pt<result.dimension(1);pt++) {
00081       result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
00082         *std::pow(points(cell,pt,2),zd);
00083     }
00084   }
00085   return;
00086 }
00087 
00088 void bcFunc( FieldContainer<double>       & result, 
00089              const FieldContainer<double> & points , 
00090              const shards::CellTopology   & parentCell ,
00091              int sideOrdinal , int comp , int a, int b, int c )
00092 {
00093 
00094   int numCells  = result.dimension(0);
00095   int numPoints = result.dimension(1);
00096 
00097   // reference face normal
00098   FieldContainer<double> normal(3);
00099   CellTools<double>::getReferenceFaceNormal(normal,sideOrdinal,parentCell);
00100 
00101   result.initialize();
00102 
00103   if (comp == 0) {
00104     for (int cell=0;cell<numCells;cell++) {
00105       for (int pt=0;pt<numPoints;pt++) {
00106         // first component
00107         if (c > 0) {
00108           result(cell,pt,0) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
00109         }
00110         if (b > 0) {
00111           result(cell,pt,0) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
00112         }
00113         // second component
00114         if (b > 0) {
00115           result(cell,pt,1) = b * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
00116         }
00117         // third component
00118         if (c > 0) {
00119           result(cell,pt,2) = c * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
00120         }
00121       }
00122     }
00123   }
00124   else if (comp == 1) {
00125     for (int cell=0;cell<numCells;cell++) {
00126       for (int pt=0;pt<numPoints;pt++) {
00127         // first component
00128         if (a > 0) {
00129           result(cell,pt,0) = a * normal(1) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00130         }
00131         // second component
00132         if (c > 0) {
00133           result(cell,pt,1) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
00134         }
00135         if (a > 0) {
00136           result(cell,pt,1) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00137         }
00138         // third component
00139         if (c > 0) {
00140           result(cell,pt,2) = c * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
00141         }       
00142       }
00143     }
00144   }
00145   else if (comp == 2) {
00146     for (int cell=0;cell<numCells;cell++) {
00147       for (int pt=0;pt<numPoints;pt++) {
00148         // first component
00149         if (a > 0) {
00150           result(cell,pt,0) = a * normal(2) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00151         }
00152         // second component
00153         if (b > 0) {
00154           result(cell,pt,1) = b * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
00155         }
00156         // third component
00157         if (b > 0) {
00158           result(cell,pt,2) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1)* std::pow(points(cell,pt,2),c);
00159         }
00160         if (a > 0) {
00161           result(cell,pt,2) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00162         }
00163       }
00164     }
00165   }
00166   
00167   return;
00168 }
00169 
00170 void rhsFunc( FieldContainer<double> & result , 
00171               const FieldContainer<double> &points ,
00172               int comp,
00173               int a,
00174               int b,
00175               int c )
00176 {
00177   // rhs is curl(curl(E))+E, so load up the exact solution
00178 
00179   u_exact(result,points,comp,a,b,c);
00180 
00181   if (comp == 0) {
00182     // component 0
00183     if (b>=2) {
00184       for (int cell=0;cell<result.dimension(0);cell++) {
00185         for (int pt=0;pt<result.dimension(1);pt++) {
00186           result(cell,pt,0) -= b*(b-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b-2)*std::pow(points(cell,pt,2),c);
00187         }
00188       }
00189     }
00190     if (c>=2) {
00191       for (int cell=0;cell<result.dimension(0);cell++) {
00192         for (int pt=0;pt<result.dimension(1);pt++) {
00193           result(cell,pt,0)-=c*(c-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
00194         }
00195       }
00196     }
00197     // component one
00198     if ( (a>=1) && (b>=1) ) {
00199       for (int cell=0;cell<result.dimension(0);cell++) {
00200         for (int pt=0;pt<result.dimension(1);pt++) {
00201           result(cell,pt,1)+=a*b*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b-1)*std::pow(points(cell,pt,2),c);
00202         }
00203       }
00204     }
00205     // component two
00206     if ( (a>=1) && (c>=1) ) {
00207       for (int cell=0;cell<result.dimension(0);cell++) {
00208         for (int pt=0;pt<result.dimension(1);pt++) {
00209           result(cell,pt,2)+=a*c*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-1);
00210         }
00211       }
00212     }
00213   }
00214   if (comp == 1) {
00215     for (int cell=0;cell<result.dimension(0);cell++) {
00216       for (int pt=0;pt<result.dimension(1);pt++) {
00217         // first component
00218         if (a > 0 && b > 0) {
00219           result(cell,pt,0) += a * b * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
00220         }
00221         // second component
00222         if (c > 1) {
00223           result(cell,pt,1) -= c*(c-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
00224         }
00225         if (a > 1) {
00226           result(cell,pt,1) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00227         }
00228         // third component
00229         if (b > 0 && c > 0) {
00230           result(cell,pt,2) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
00231         }
00232       }
00233     }
00234   }
00235   else if (comp == 2) {
00236     for (int cell=0;cell<result.dimension(0);cell++) {
00237       for (int pt=0;pt<result.dimension(1);pt++) {
00238         // first component
00239         if (a>0 && c>0) {
00240           result(cell,pt,0) += a * c * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
00241         }
00242         // second component
00243         if (b>0 && c>0) {
00244           result(cell,pt,1) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
00245         }
00246         // third component
00247         if (b>1) {
00248           result(cell,pt,2) -= b*(b-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-2) * std::pow(points(cell,pt,2),c);
00249         }
00250         if (a>1) {
00251           result(cell,pt,2) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
00252         }
00253       }
00254     }
00255   }
00256   return;
00257 }
00258 
00259 int main(int argc, char *argv[]) {
00260   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00261   
00262   // This little trick lets us print to std::cout only if
00263   // a (dummy) command-line argument is provided.
00264   int iprint = argc - 1;
00265   Teuchos::RCP<std::ostream> outStream;
00266   Teuchos::oblackholestream bhs; // outputs nothing
00267   if (iprint > 0)
00268     outStream = Teuchos::rcp(&std::cout, false);
00269   else
00270     outStream = Teuchos::rcp(&bhs, false);
00271   
00272   // Save the format state of the original std::cout.
00273   Teuchos::oblackholestream oldFormatState;
00274   oldFormatState.copyfmt(std::cout);
00275   
00276   *outStream \
00277     << "===============================================================================\n" \
00278     << "|                                                                             |\n" \
00279     << "|                  Unit Test (Basis_HGRAD_TRI_In_FEM)                         |\n" \
00280     << "|                                                                             |\n" \
00281     << "| 1) Patch test involving H(curl) matrices                                    |\n" \
00282     << "|                                                                             |\n" \
00283     << "|   Questions? Contact Pavel Bochev (pbboche@sandia.gov),                     |\n" \
00284     << "|                      Robert Kirby (robert.c.kirby@ttu.edu),                 |\n" \
00285     << "|                      Denis Ridzal (dridzal@sandia.gov),                     |\n" \
00286     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00287     << "|                                                                             |\n" \
00288     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00289     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00290     << "|                                                                             |\n" \
00291     << "===============================================================================\n" \
00292     << "| TEST 1: Patch test for mass + curl-curl matrices                            |\n" \
00293     << "===============================================================================\n";
00294   
00295   
00296   int errorFlag = 0;
00297   
00298   outStream -> precision(16);
00299   
00300   try {
00301     DefaultCubatureFactory<double> cubFactory;                                          // create cubature factory
00302     shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());  // create parent cell topology
00303 
00304     shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());
00305     
00306     int cellDim = cell.getDimension();
00307     int sideDim = side.getDimension();
00308     
00309     int min_order = 1;
00310     int max_order = 4;
00311     
00312     int numIntervals = max_order;
00313     int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
00314     FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
00315     int counter = 0;
00316     for (int j=0; j<=numIntervals; j++) {
00317       for (int i=0; i<=numIntervals-j; i++) {
00318         for (int k=0;k<numIntervals-j-i;k++) {
00319           interp_points_ref(counter,0) = i*(1.0/numIntervals);
00320           interp_points_ref(counter,1) = j*(1.0/numIntervals);
00321           interp_points_ref(counter,2) = k*(1.0/numIntervals);
00322           counter++;
00323         }
00324       }
00325     }
00326     
00327     for (int basis_order=min_order;basis_order<=max_order;basis_order++) {
00328       // create basis
00329       Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
00330         Teuchos::rcp(new Basis_HCURL_TET_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_EQUISPACED) );
00331       
00332       int numFields = basis->getCardinality();
00333       
00334       // create cubatures
00335       Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
00336       Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
00337       
00338       int numCubPointsCell = cellCub->getNumPoints();
00339       int numCubPointsSide = sideCub->getNumPoints();
00340       
00341       // hold cubature information
00342       FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
00343       FieldContainer<double> cub_weights_cell(numCubPointsCell);
00344       FieldContainer<double> cub_points_side(numCubPointsSide,sideDim);
00345       FieldContainer<double> cub_weights_side(numCubPointsSide);
00346       FieldContainer<double> cub_points_side_refcell(numCubPointsSide,cellDim);
00347 
00348       // hold basis function information on refcell
00349       FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
00350       FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
00351       FieldContainer<double> curl_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
00352       FieldContainer<double> w_curl_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
00353       FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields,numCubPointsSide,cellDim );
00354       FieldContainer<double> w_value_of_basis_at_cub_points_side_refcell(1,numFields,numCubPointsSide,cellDim );
00355 
00356 
00357       // holds rhs data
00358       FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim);
00359       FieldContainer<double> bc_func_at_cub_points_side_refcell(1,numCubPointsSide,cellDim);
00360       FieldContainer<double> bc_fields_per_side(1,numFields);
00361 
00362       // FEM mass matrix
00363       FieldContainer<double> fe_matrix_bak(1,numFields,numFields);
00364       FieldContainer<double> fe_matrix(1,numFields,numFields);
00365       FieldContainer<double> rhs_and_soln_vec(1,numFields);
00366       
00367       FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim);
00368       FieldContainer<double> interpolant( 1, numInterpPoints , cellDim );
00369 
00370       int info = 0;
00371       Teuchos::LAPACK<int, double> solver;
00372       
00373       // set test tolerance
00374       double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
00375       
00376       // build matrices outside the loop, and then just do the rhs
00377       // for each iteration
00378       cellCub->getCubature(cub_points_cell, cub_weights_cell);
00379       sideCub->getCubature(cub_points_side, cub_weights_side);
00380 
00381       // need the vector basis
00382       basis->getValues(value_of_basis_at_cub_points_cell,
00383                        cub_points_cell,
00384                        OPERATOR_VALUE);
00385       basis->getValues(curl_of_basis_at_cub_points_cell,
00386                        cub_points_cell,
00387                        OPERATOR_CURL);
00388 
00389       basis->getValues( value_of_basis_at_interp_points ,
00390                         interp_points_ref ,
00391                         OPERATOR_VALUE );
00392 
00393 
00394       // construct mass and curl-curl matrices
00395       cub_weights_cell.resize(1,numCubPointsCell);
00396       FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
00397                                                   cub_weights_cell ,
00398                                                   value_of_basis_at_cub_points_cell ); 
00399       FunctionSpaceTools::multiplyMeasure<double>(w_curl_of_basis_at_cub_points_cell ,
00400                                                   cub_weights_cell ,
00401                                                   curl_of_basis_at_cub_points_cell ); 
00402       cub_weights_cell.resize(numCubPointsCell);
00403       
00404       
00405       value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
00406       curl_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
00407       FunctionSpaceTools::integrate<double>(fe_matrix_bak,
00408                                             w_value_of_basis_at_cub_points_cell ,
00409                                             value_of_basis_at_cub_points_cell ,
00410                                             COMP_BLAS );
00411       FunctionSpaceTools::integrate<double>(fe_matrix_bak,
00412                                             w_curl_of_basis_at_cub_points_cell ,
00413                                             curl_of_basis_at_cub_points_cell ,
00414                                             COMP_BLAS ,
00415                                             true );
00416       value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
00417       curl_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
00418 
00419       for (int comp=0;comp<cellDim;comp++) {
00420         for (int x_order=0;x_order<basis_order;x_order++) {
00421           for (int y_order=0;y_order<basis_order-x_order;y_order++) {
00422             for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
00423               fe_matrix.initialize();
00424               // copy mass matrix 
00425               for (int i=0;i<numFields;i++) {
00426                 for (int j=0;j<numFields;j++) {
00427                   fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
00428                 }
00429               }
00430               
00431               // clear old vector data
00432               rhs_and_soln_vec.initialize();
00433               
00434               // now get rhs vector
00435               cub_points_cell.resize(1,numCubPointsCell,cellDim);
00436              
00437               rhs_at_cub_points_cell.initialize();
00438               rhsFunc(rhs_at_cub_points_cell,
00439                       cub_points_cell,
00440                       comp, 
00441                       x_order,
00442                       y_order,
00443                       z_order);
00444             
00445               cub_points_cell.resize(numCubPointsCell,cellDim);
00446 
00447               cub_weights_cell.resize(numCubPointsCell);
00448 
00449               FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
00450                                                     rhs_at_cub_points_cell,
00451                                                     w_value_of_basis_at_cub_points_cell,
00452                                                     COMP_BLAS);
00453               
00454               // now I need to get the boundary condition terms in place
00455 
00456               for (unsigned i=0;i<4;i++) {
00457                 // map side quadrature to reference cell side
00458                 CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell,cub_points_side,sideDim,
00459                                                          (int) i, cell);
00460 
00461                 //evaluate basis at these points
00462                 basis->getValues( value_of_basis_at_cub_points_side_refcell ,
00463                                   cub_points_side_refcell ,
00464                                   OPERATOR_VALUE );
00465 
00466                 // evaluate imposed current on surface, which is n x curl(u_exact), on the quad points
00467                 cub_points_side_refcell.resize( 1 , numCubPointsSide , cellDim );
00468                 bcFunc(bc_func_at_cub_points_side_refcell,
00469                        cub_points_side_refcell, 
00470                        cell ,
00471                        i,
00472                        comp ,
00473                        x_order,
00474                        y_order,
00475                        z_order);
00476                 cub_points_side_refcell.resize( numCubPointsSide , cellDim );
00477 
00478                 // now I need to integrate the bc function against the test basis
00479                 // need to weight the basis functions with quadrature weights
00480                 // the size of the face is embedded in the normal
00481                 cub_weights_side.resize(1,numCubPointsSide);
00482                 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_side_refcell ,
00483                                                             cub_weights_side ,
00484                                                             value_of_basis_at_cub_points_side_refcell ); 
00485                 cub_weights_side.resize(numCubPointsSide);
00486 
00487                 FunctionSpaceTools::integrate<double>( bc_fields_per_side , 
00488                                                        bc_func_at_cub_points_side_refcell ,
00489                                                        w_value_of_basis_at_cub_points_side_refcell ,
00490                                                        COMP_BLAS );
00491 
00492                 RealSpaceTools<double>::subtract(rhs_and_soln_vec, bc_fields_per_side );
00493                 
00494 
00495               }
00496 
00497               // solve linear system
00498               solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info);
00499               solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
00500               
00501               interp_points_ref.resize(1,numInterpPoints,cellDim);
00502               // get exact solution for comparison
00503               FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
00504               exact_solution.initialize();
00505               u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
00506               interp_points_ref.resize(numInterpPoints,cellDim);
00507 
00508               // compute interpolant
00509               // first evaluate basis at interpolation points
00510               value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
00511               FunctionSpaceTools::evaluate<double>( interpolant , 
00512                                                     rhs_and_soln_vec ,
00513                                                     value_of_basis_at_interp_points );
00514               value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
00515               
00516               RealSpaceTools<double>::subtract(interpolant,exact_solution);
00517               
00518               double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO);
00519               
00520               *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
00521                          << x_order << ", " << y_order << ", " << z_order
00522                          << ") in component " << comp 
00523                          << " and finite element interpolant of order " << basis_order << ": "
00524                          << nrm << "\n";
00525               
00526               if (nrm > zero) {
00527                 *outStream << "\n\nPatch test failed for solution polynomial order ("
00528                            << x_order << ", " << y_order << ", " << z_order << ") and basis order "
00529                            << basis_order << "\n\n";
00530                 errorFlag++;
00531               }
00532             }
00533           }
00534         }
00535       }
00536     }
00537     
00538   }
00539   
00540   catch (std::logic_error err) {
00541     *outStream << err.what() << "\n\n";
00542     errorFlag = -1000;
00543   };
00544   
00545   if (errorFlag != 0)
00546     std::cout << "End Result: TEST FAILED\n";
00547   else
00548     std::cout << "End Result: TEST PASSED\n";
00549   
00550   // reset format state of std::cout
00551   std::cout.copyfmt(oldFormatState);
00552   
00553   return errorFlag;
00554 }