Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_LINE_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_LINE_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);
00068 void neumann(FieldContainer<double> &, const FieldContainer<double> &, const FieldContainer<double> &, int);
00069 void u_exact(FieldContainer<double> &, const FieldContainer<double> &, int);
00070 
00072 void rhsFunc(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
00073   if (degree == 0) {
00074     for (int cell=0; cell<result.dimension(0); cell++) {
00075       for (int pt=0; pt<result.dimension(1); pt++) {
00076         result(cell,pt) = 1.0;
00077 
00078       }
00079     }
00080   }
00081   else if (degree == 1) {
00082     for (int cell=0; cell<result.dimension(0); cell++) {
00083       for (int pt=0; pt<result.dimension(1); pt++) {
00084         result(cell,pt) = points(cell,pt,0);
00085       }
00086     }
00087   }
00088   else {
00089     for (int cell=0; cell<result.dimension(0); cell++) {
00090       for (int pt=0; pt<result.dimension(1); pt++) {
00091         result(cell,pt) = pow(points(cell,pt,0), degree) - degree*(degree-1)*pow(points(cell,pt,0), degree-2);
00092       }
00093     }
00094   }
00095 }
00096 
00098 void neumann(FieldContainer<double> & g_phi, const FieldContainer<double> & phi1, const FieldContainer<double> & phi2, int degree) {
00099   double g_at_one, g_at_minus_one;
00100   int num_fields;
00101 
00102   if (degree == 0) {
00103     g_at_one = 0.0;
00104     g_at_minus_one = 0.0;
00105   }
00106   else {
00107     g_at_one = degree*pow(1.0, degree-1);
00108     g_at_minus_one = degree*pow(-1.0, degree-1);
00109   }
00110 
00111   num_fields = phi1.dimension(0);
00112 
00113   for (int i=0; i<num_fields; i++) {
00114     g_phi(0,i) = g_at_minus_one*phi1(i,0);
00115     g_phi(1,i) = g_at_one*phi2(i,0);
00116   }
00117 }
00118 
00120 void u_exact(FieldContainer<double> & result, const FieldContainer<double> & points, int degree) {
00121   for (int cell=0; cell<result.dimension(0); cell++) {
00122     for (int pt=0; pt<result.dimension(1); pt++) {
00123       result(cell,pt) = pow(points(pt,0), degree);
00124     }
00125   }
00126 }
00127 
00128 
00129 
00130 
00131 int main(int argc, char *argv[]) {
00132 
00133   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00134 
00135   // This little trick lets us print to std::cout only if
00136   // a (dummy) command-line argument is provided.
00137   int iprint     = argc - 1;
00138   Teuchos::RCP<std::ostream> outStream;
00139   Teuchos::oblackholestream bhs; // outputs nothing
00140   if (iprint > 0)
00141     outStream = Teuchos::rcp(&std::cout, false);
00142   else
00143     outStream = Teuchos::rcp(&bhs, false);
00144 
00145   // Save the format state of the original std::cout.
00146   Teuchos::oblackholestream oldFormatState;
00147   oldFormatState.copyfmt(std::cout);
00148 
00149   *outStream \
00150     << "===============================================================================\n" \
00151     << "|                                                                             |\n" \
00152     << "|               Unit Test (Basis_HGRAD_LINE_Cn_FEM)                           |\n" \
00153     << "|                                                                             |\n" \
00154     << "|     1) Patch test involving mass and stiffness matrices,                    |\n" \
00155     << "|        for the Neumann problem on a REFERENCE line:                         |\n" \
00156     << "|                                                                             |\n" \
00157     << "|            - u'' + u = f  in (-1,1),  u' = g at -1,1                        |\n" \
00158     << "|                                                                             |\n" \
00159     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00160     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00161     << "|                      Robert Kirby  (robert.c.kirby@ttu.gov),                |\n" \
00162     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00163     << "|                                                                             |\n" \
00164     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00165     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00166     << "|                                                                             |\n" \
00167     << "===============================================================================\n"\
00168     << "| TEST 1: Patch test                                                          |\n"\
00169     << "===============================================================================\n";
00170 
00171   
00172   int errorFlag = 0;
00173   double zero = 100*INTREPID_TOL;
00174   outStream -> precision(20);
00175 
00176 
00177   try {
00178 
00179     int max_order = 10;  // max total order of polynomial solution
00180 
00181     // Define array containing points at which the solution is evaluated
00182     int numIntervals = 100;
00183     int numInterpPoints = numIntervals + 1;
00184     FieldContainer<double> interp_points(numInterpPoints, 1);
00185     for (int i=0; i<numInterpPoints; i++) {
00186       interp_points(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
00187     }
00188     
00189     DefaultCubatureFactory<double>  cubFactory;                                   // create factory
00190     shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());   // create cell topology
00191 
00192 
00193     for (int basis_order=1; basis_order <= max_order; basis_order++) {
00194       //create basis
00195       Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
00196         Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_SPECTRAL) );
00197 
00198       for (int soln_order=1; soln_order <= basis_order; soln_order++) {
00199         // evaluate exact solution
00200         FieldContainer<double> exact_solution(1, numInterpPoints);
00201         u_exact(exact_solution, interp_points, soln_order);
00202         
00203         int numFields = lineBasis->getCardinality();
00204 
00205         // create cubature
00206         Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, 2*basis_order-2);
00207         int numCubPoints = lineCub->getNumPoints();
00208         int spaceDim = lineCub->getDimension();
00209 
00210         /* Computational arrays. */
00211         FieldContainer<double> cub_points(numCubPoints, spaceDim);
00212         FieldContainer<double> cub_points_physical(1, numCubPoints, spaceDim);
00213         FieldContainer<double> cub_weights(numCubPoints);
00214         FieldContainer<double> cell_nodes(1, 2, spaceDim);
00215         FieldContainer<double> jacobian(1, numCubPoints, spaceDim, spaceDim);
00216         FieldContainer<double> jacobian_inv(1, numCubPoints, spaceDim, spaceDim);
00217         FieldContainer<double> jacobian_det(1, numCubPoints);
00218         FieldContainer<double> weighted_measure(1, numCubPoints);
00219 
00220         FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints);
00221         FieldContainer<double> transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
00222         FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
00223         FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
00224         FieldContainer<double> transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
00225         FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
00226         FieldContainer<double> fe_matrix(1, numFields, numFields);
00227 
00228         FieldContainer<double> rhs_at_cub_points_physical(1, numCubPoints);
00229         FieldContainer<double> rhs_and_soln_vector(1, numFields);
00230 
00231         FieldContainer<double> one_point(1, 1);
00232         FieldContainer<double> value_of_basis_at_one(numFields, 1);
00233         FieldContainer<double> value_of_basis_at_minusone(numFields, 1);
00234         FieldContainer<double> bc_neumann(2, numFields);
00235 
00236         FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
00237         FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
00238         FieldContainer<double> interpolant(1, numInterpPoints);
00239 
00240         FieldContainer<int> ipiv(numFields);
00241 
00242         /******************* START COMPUTATION ***********************/
00243 
00244         // get cubature points and weights
00245         lineCub->getCubature(cub_points, cub_weights);
00246 
00247         // fill cell vertex array
00248         cell_nodes(0, 0, 0) = -1.0;
00249         cell_nodes(0, 1, 0) = 1.0;
00250 
00251         // compute geometric cell information
00252         CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, line);
00253         CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
00254         CellTools<double>::setJacobianDet(jacobian_det, jacobian);
00255 
00256         // compute weighted measure
00257         FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
00258 
00260         // Computing mass matrices:
00261         // tabulate values of basis functions at (reference) cubature points
00262         lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
00263 
00264         // transform values of basis functions
00265         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
00266                                                         value_of_basis_at_cub_points);
00267 
00268         // multiply with weighted measure
00269         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
00270                                                     weighted_measure,
00271                                                     transformed_value_of_basis_at_cub_points);
00272 
00273         // compute mass matrices
00274         FunctionSpaceTools::integrate<double>(fe_matrix,
00275                                               transformed_value_of_basis_at_cub_points,
00276                                               weighted_transformed_value_of_basis_at_cub_points,
00277                                               COMP_CPP);
00279 
00281         // Computing stiffness matrices:
00282         // tabulate gradients of basis functions at (reference) cubature points
00283         lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
00284 
00285         // transform gradients of basis functions
00286         FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
00287                                                        jacobian_inv,
00288                                                        grad_of_basis_at_cub_points);
00289 
00290         // multiply with weighted measure
00291         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
00292                                                     weighted_measure,
00293                                                     transformed_grad_of_basis_at_cub_points);
00294 
00295         // compute stiffness matrices and sum into fe_matrix
00296         FunctionSpaceTools::integrate<double>(fe_matrix,
00297                                               transformed_grad_of_basis_at_cub_points,
00298                                               weighted_transformed_grad_of_basis_at_cub_points,
00299                                               COMP_CPP,
00300                                               true);
00302 
00304         // Computing RHS contributions:
00305         // map (reference) cubature points to physical space
00306         CellTools<double>::mapToPhysicalFrame(cub_points_physical, cub_points, cell_nodes, line);
00307 
00308         // evaluate rhs function
00309         rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
00310 
00311         // compute rhs
00312         FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
00313                                               rhs_at_cub_points_physical,
00314                                               weighted_transformed_value_of_basis_at_cub_points,
00315                                               COMP_CPP);
00316 
00317         // compute neumann b.c. contributions and adjust rhs
00318         one_point(0,0) = 1.0;   lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
00319         one_point(0,0) = -1.0;  lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
00320         neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
00321         for (int i=0; i<numFields; i++) {
00322           rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
00323           rhs_and_soln_vector(0, i) += bc_neumann(1, i);
00324         }
00326 
00328         // Solution of linear system:
00329         int info = 0;
00330         Teuchos::LAPACK<int, double> solver;
00331         //solver.GESV(numRows, 1, &fe_mat(0,0), numRows, &ipiv(0), &fe_vec(0), numRows, &info);
00332         solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
00334 
00336         // Building interpolant:
00337         // evaluate basis at interpolation points
00338         lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
00339         // transform values of basis functions
00340         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
00341                                                         value_of_basis_at_interp_points);
00342         FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
00344 
00345         /******************* END COMPUTATION ***********************/
00346       
00347         RealSpaceTools<double>::subtract(interpolant, exact_solution);
00348 
00349         *outStream << "\nNorm-2 difference between exact solution polynomial of order "
00350                    << soln_order << " and finite element interpolant of order " << basis_order << ": "
00351                    << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) << "\n";
00352 
00353         if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) > zero) {
00354           *outStream << "\n\nPatch test failed for solution polynomial order "
00355                      << soln_order << " and basis order " << basis_order << "\n\n";
00356           errorFlag++;
00357         }
00358 
00359       } // end for basis_order
00360 
00361    } // end for soln_order
00362 
00363   }
00364   // Catch unexpected errors
00365   catch (std::logic_error err) {
00366     *outStream << err.what() << "\n\n";
00367     errorFlag = -1000;
00368   };
00369 
00370   if (errorFlag != 0)
00371     std::cout << "End Result: TEST FAILED\n";
00372   else
00373     std::cout << "End Result: TEST PASSED\n";
00374 
00375   // reset format state of std::cout
00376   std::cout.copyfmt(oldFormatState);
00377 
00378   return errorFlag;
00379 }