Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_LINE_Cn_FEM_JACOBI/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_JACOBI.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_Cn_FEM_JACOBI)                    |\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 = 10;  // 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     for (int soln_order=1; soln_order <= max_order; soln_order++) {
00190 
00191       // evaluate exact solution
00192       FieldContainer<double> exact_solution(1, numInterpPoints);
00193       u_exact(exact_solution, interp_points, soln_order);
00194 
00195       for (int basis_order=soln_order; basis_order <= max_order; basis_order++) {
00196 
00197         //create basis
00198         Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
00199           Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM_JACOBI<double,FieldContainer<double> >(basis_order) );
00200         int numFields = lineBasis->getCardinality();
00201 
00202         // create cubature
00203         Teuchos::RCP<Cubature<double> > lineCub = cubFactory.create(line, 2*basis_order-2);
00204         int numCubPoints = lineCub->getNumPoints();
00205         int spaceDim = lineCub->getDimension();
00206 
00207         /* Computational arrays. */
00208         FieldContainer<double> cub_points(numCubPoints, spaceDim);
00209         FieldContainer<double> cub_points_physical(1, numCubPoints, spaceDim);
00210         FieldContainer<double> cub_weights(numCubPoints);
00211         FieldContainer<double> cell_nodes(1, 2, spaceDim);
00212         FieldContainer<double> jacobian(1, numCubPoints, spaceDim, spaceDim);
00213         FieldContainer<double> jacobian_inv(1, numCubPoints, spaceDim, spaceDim);
00214         FieldContainer<double> jacobian_det(1, numCubPoints);
00215         FieldContainer<double> weighted_measure(1, numCubPoints);
00216 
00217         FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints);
00218         FieldContainer<double> transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
00219         FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(1, numFields, numCubPoints);
00220         FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
00221         FieldContainer<double> transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
00222         FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
00223         FieldContainer<double> fe_matrix(1, numFields, numFields);
00224 
00225         FieldContainer<double> rhs_at_cub_points_physical(1, numCubPoints);
00226         FieldContainer<double> rhs_and_soln_vector(1, numFields);
00227 
00228         FieldContainer<double> one_point(1, 1);
00229         FieldContainer<double> value_of_basis_at_one(numFields, 1);
00230         FieldContainer<double> value_of_basis_at_minusone(numFields, 1);
00231         FieldContainer<double> bc_neumann(2, numFields);
00232 
00233         FieldContainer<double> value_of_basis_at_interp_points(numFields, numInterpPoints);
00234         FieldContainer<double> transformed_value_of_basis_at_interp_points(1, numFields, numInterpPoints);
00235         FieldContainer<double> interpolant(1, numInterpPoints);
00236 
00237         FieldContainer<int> ipiv(numFields);
00238 
00239         /******************* START COMPUTATION ***********************/
00240 
00241         // get cubature points and weights
00242         lineCub->getCubature(cub_points, cub_weights);
00243 
00244         // fill cell vertex array
00245         cell_nodes(0, 0, 0) = -1.0;
00246         cell_nodes(0, 1, 0) = 1.0;
00247 
00248         // compute geometric cell information
00249         CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, line);
00250         CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
00251         CellTools<double>::setJacobianDet(jacobian_det, jacobian);
00252 
00253         // compute weighted measure
00254         FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
00255 
00257         // Computing mass matrices:
00258         // tabulate values of basis functions at (reference) cubature points
00259         lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
00260 
00261         // transform values of basis functions
00262         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
00263                                                         value_of_basis_at_cub_points);
00264 
00265         // multiply with weighted measure
00266         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
00267                                                     weighted_measure,
00268                                                     transformed_value_of_basis_at_cub_points);
00269 
00270         // compute mass matrices
00271         FunctionSpaceTools::integrate<double>(fe_matrix,
00272                                               transformed_value_of_basis_at_cub_points,
00273                                               weighted_transformed_value_of_basis_at_cub_points,
00274                                               COMP_CPP);
00276 
00278         // Computing stiffness matrices:
00279         // tabulate gradients of basis functions at (reference) cubature points
00280         lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
00281 
00282         // transform gradients of basis functions
00283         FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
00284                                                        jacobian_inv,
00285                                                        grad_of_basis_at_cub_points);
00286 
00287         // multiply with weighted measure
00288         FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
00289                                                     weighted_measure,
00290                                                     transformed_grad_of_basis_at_cub_points);
00291 
00292         // compute stiffness matrices and sum into fe_matrix
00293         FunctionSpaceTools::integrate<double>(fe_matrix,
00294                                               transformed_grad_of_basis_at_cub_points,
00295                                               weighted_transformed_grad_of_basis_at_cub_points,
00296                                               COMP_CPP,
00297                                               true);
00299 
00301         // Computing RHS contributions:
00302         // map (reference) cubature points to physical space
00303         CellTools<double>::mapToPhysicalFrame(cub_points_physical, cub_points, cell_nodes, line);
00304 
00305         // evaluate rhs function
00306         rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
00307 
00308         // compute rhs
00309         FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
00310                                               rhs_at_cub_points_physical,
00311                                               weighted_transformed_value_of_basis_at_cub_points,
00312                                               COMP_CPP);
00313 
00314         // compute neumann b.c. contributions and adjust rhs
00315         one_point(0,0) = 1.0;   lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
00316         one_point(0,0) = -1.0;  lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
00317         neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
00318         for (int i=0; i<numFields; i++) {
00319           rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
00320           rhs_and_soln_vector(0, i) += bc_neumann(1, i);
00321         }
00323 
00325         // Solution of linear system:
00326         int info = 0;
00327         Teuchos::LAPACK<int, double> solver;
00328         //solver.GESV(numRows, 1, &fe_mat(0,0), numRows, &ipiv(0), &fe_vec(0), numRows, &info);
00329         solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
00331 
00333         // Building interpolant:
00334         // evaluate basis at interpolation points
00335         lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
00336         // transform values of basis functions
00337         FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
00338                                                         value_of_basis_at_interp_points);
00339         FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
00341 
00342         /******************* END COMPUTATION ***********************/
00343       
00344         RealSpaceTools<double>::subtract(interpolant, exact_solution);
00345 
00346         *outStream << "\nNorm-2 difference between exact solution polynomial of order "
00347                    << soln_order << " and finite element interpolant of order " << basis_order << ": "
00348                    << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) << "\n";
00349 
00350         if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) > zero) {
00351           *outStream << "\n\nPatch test failed for solution polynomial order "
00352                      << soln_order << " and basis order " << basis_order << "\n\n";
00353           errorFlag++;
00354         }
00355 
00356       } // end for basis_order
00357 
00358    } // end for soln_order
00359 
00360   }
00361   // Catch unexpected errors
00362   catch (std::logic_error err) {
00363     *outStream << err.what() << "\n\n";
00364     errorFlag = -1000;
00365   };
00366 
00367   if (errorFlag != 0)
00368     std::cout << "End Result: TEST FAILED\n";
00369   else
00370     std::cout << "End Result: TEST PASSED\n";
00371 
00372   // reset format state of std::cout
00373   std::cout.copyfmt(oldFormatState);
00374 
00375   return errorFlag;
00376 }