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