Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/example/Drivers/example_04.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 
00066 // Intrepid includes
00067 #include "Intrepid_FunctionSpaceTools.hpp"
00068 #include "Intrepid_FieldContainer.hpp"
00069 #include "Intrepid_CellTools.hpp"
00070 #include "Intrepid_ArrayTools.hpp"
00071 #include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
00072 #include "Intrepid_HGRAD_QUAD_Cn_FEM.hpp"
00073 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp"
00074 #include "Intrepid_RealSpaceTools.hpp"
00075 #include "Intrepid_CubaturePolylib.hpp"
00076 #include "Intrepid_CubatureTensor.hpp"
00077 #include "Intrepid_Utils.hpp"
00078 
00079 // Epetra includes
00080 #include "Epetra_Time.h"
00081 #include "Epetra_Map.h"
00082 #include "Epetra_FECrsMatrix.h"
00083 #include "Epetra_FEVector.h"
00084 #include "Epetra_SerialComm.h"
00085 
00086 // Teuchos includes
00087 #include "Teuchos_oblackholestream.hpp"
00088 #include "Teuchos_RCP.hpp"
00089 #include "Teuchos_BLAS.hpp"
00090 
00091 // Shards includes
00092 #include "Shards_CellTopology.hpp"
00093 
00094 // EpetraExt includes
00095 #include "EpetraExt_RowMatrixOut.h"
00096 #include "EpetraExt_MultiVectorOut.h"
00097 
00098 using namespace std;
00099 using namespace Intrepid;
00100 
00101 int main(int argc, char *argv[]) {
00102 
00103   //Check number of arguments
00104     if (argc < 2) {
00105        std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
00106        std::cout <<"Usage:\n\n";
00107        std::cout <<"  ./Intrepid_example_Drivers_Example_05.exe min_deg max_deg verbose\n\n";
00108        std::cout <<" where \n";
00109        std::cout <<"   int min_deg         - beginning polynomial degree to check \n";
00110        std::cout <<"   int max_deg         - final polynomial degree to check \n";
00111        std::cout <<"   verbose (optional)  - any character, indicates verbose output \n\n";
00112        exit(1);
00113     }
00114     
00115   // This little trick lets us print to std::cout only if
00116   // a (dummy) command-line argument is provided.
00117   int iprint     = argc - 1;
00118   Teuchos::RCP<std::ostream> outStream;
00119   Teuchos::oblackholestream bhs; // outputs nothing
00120   if (iprint > 1)
00121     outStream = Teuchos::rcp(&std::cout, false);
00122   else
00123     outStream = Teuchos::rcp(&bhs, false);
00124   
00125   // Save the format state of the original std::cout.
00126   Teuchos::oblackholestream oldFormatState;
00127   oldFormatState.copyfmt(std::cout);
00128   
00129   *outStream \
00130     << "===============================================================================\n" \
00131     << "|                                                                             |\n" \
00132     << "|  Example: Check diagonalization of reference mass matrix                    |\n" \
00133     << "|           on line, quad, and hex                                            |\n" \
00134     << "|                                                                             |\n" \
00135     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00136     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00137     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00138     << "|                                                                             |\n" \
00139     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00140     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00141     << "|                                                                             |\n" \
00142     << "===============================================================================\n";
00143 
00144 
00145 // ************************************ GET INPUTS **************************************
00146   int min_degree = atoi(argv[1]);
00147   int max_degree = atoi(argv[2]);
00148   
00149 
00150 // *********************************** CELL TOPOLOGY **********************************
00151   
00152   // Get cell topology for base line
00153   typedef shards::CellTopology    CellTopology;
00154   CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
00155   CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00156   CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
00157   
00158   std::vector<CellTopology> cts(3);
00159 
00160   // loop over degrees
00161   for (int deg=min_degree;deg<=max_degree;deg++) {
00162     std::vector<Teuchos::RCP<Basis<double,FieldContainer<double> > > > bases;
00163     bases.push_back( Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
00164     bases.push_back( Teuchos::rcp(new Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
00165     bases.push_back( Teuchos::rcp(new Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
00166 
00167     // ************************************ CUBATURE **************************************
00168     Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
00169       = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
00170     
00171     std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
00172       cub_to_tensor;
00173 
00174     // now we loop over spatial dimensions
00175     for (int sd=1;sd<=3;sd++) {
00176       // ************************************ CUBATURE **************************************
00177       cub_to_tensor.push_back( glcub );
00178       CubatureTensor<double,FieldContainer<double>,FieldContainer<double> > cubcur( cub_to_tensor );
00179 
00180       int cubDim       = cubcur.getDimension();
00181       int numCubPoints = cubcur.getNumPoints();
00182     
00183       FieldContainer<double> cubPoints(numCubPoints, cubDim);
00184       FieldContainer<double> cubWeights(numCubPoints);
00185     
00186       cubcur.getCubature(cubPoints, cubWeights);
00187 
00188       // ************************************ BASIS AT QP**************************************
00189       Teuchos::RCP<Basis<double,FieldContainer<double> > > basis_cur = bases[sd-1];
00190       const int numFields = basis_cur->getCardinality();
00191       FieldContainer<double> bf_at_cub_pts( numFields, numCubPoints );
00192       FieldContainer<double> trans_bf_at_cub_pts( 1 , numFields,numCubPoints );
00193       FieldContainer<double> w_trans_bf_at_cub_pts( 1, numFields , numCubPoints );
00194       basis_cur->getValues( bf_at_cub_pts , cubPoints , OPERATOR_VALUE );
00195       
00196       // *********************************** FORM MASS MATRIX *****************************
00197        FunctionSpaceTools::HGRADtransformVALUE<double>( trans_bf_at_cub_pts ,
00198                                                         bf_at_cub_pts );
00199        cubWeights.resize(1,numCubPoints);
00200        FunctionSpaceTools::multiplyMeasure<double>( w_trans_bf_at_cub_pts ,
00201                                                    cubWeights ,
00202                                                    trans_bf_at_cub_pts );
00203        cubWeights.resize(numCubPoints);
00204 
00205       FieldContainer<double> mass_matrix( 1 , numFields, numFields );
00206       FunctionSpaceTools::integrate<double>( mass_matrix ,
00207                                              trans_bf_at_cub_pts ,
00208                                              w_trans_bf_at_cub_pts ,
00209                                              COMP_BLAS );
00210 
00211       // now we loop over the mass matrix and compare the nondiagonal entries to zero
00212       double max_offdiag = 0.0;
00213       for (int i=0;i<numFields;i++) {
00214         for (int j=0;j<numFields;j++) {
00215           if (i != j) {
00216             if ( abs(mass_matrix(0,i,j)) >= max_offdiag) {
00217               max_offdiag = abs(mass_matrix(0,i,j));
00218             }
00219           }
00220         }
00221       }
00222       double min_diag = mass_matrix(0,0,0);
00223       for (int i=0;i<numFields;i++) {
00224         if ( mass_matrix(0,i,i) <= min_diag ) {
00225           min_diag = mass_matrix(0,i,i);
00226         }
00227       }
00228       *outStream << "Degree = " << deg << " and dimension = " << sd << "; Max offdiagonal"
00229                  << " element is " << max_offdiag << " and min diagonal element is " << min_diag
00230                  << ".\n";
00231           
00232     }
00233 
00234   }    
00235     
00236    std::cout << "End Result: TEST PASSED\n";
00237 
00238   std::cout.copyfmt(oldFormatState);
00239 
00240 return 0;
00241 }
00242