Intrepid
http://trilinos.sandia.gov/packages/docs/r11.0/packages/intrepid/test/Discretization/Basis/HGRAD_TET_Cn_FEM/test_01.cpp
Go to the documentation of this file.
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 
00044 
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Teuchos_oblackholestream.hpp"
00052 #include "Teuchos_RCP.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 #include "Intrepid_PointTools.hpp"
00055 #include "Intrepid_HGRAD_TET_Cn_FEM.hpp"
00056 #include "Shards_CellTopology.hpp"
00057 
00058 #include <iostream>
00059 using namespace Intrepid;
00060 
00061 
00067 int main(int argc, char *argv[]) {
00068 
00069   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00070   
00071   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00072   int iprint     = argc - 1;
00073   
00074   Teuchos::RCP<std::ostream> outStream;
00075   Teuchos::oblackholestream bhs; // outputs nothing
00076   
00077   if (iprint > 0)
00078     outStream = Teuchos::rcp(&std::cout, false);
00079   else
00080     outStream = Teuchos::rcp(&bhs, false);
00081   
00082   // Save the format state of the original std::cout.
00083   Teuchos::oblackholestream oldFormatState;
00084   oldFormatState.copyfmt(std::cout);
00085   
00086   *outStream \
00087     << "===============================================================================\n" \
00088     << "|                                                                             |\n" \
00089     << "|                           Unit Test HGRAD_TET_Cn_FEM                        |\n" \
00090     << "|                                                                             |\n" \
00091     << "|     1) Tests triangular Lagrange basis                                      |\n" \
00092     << "|                                                                             |\n" \
00093     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00094     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00095     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00096     << "|                                                                             |\n" \
00097     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00098     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00099     << "|                                                                             |\n" \
00100     << "===============================================================================\n";
00101   
00102   int errorFlag  = 0;
00103 
00104   // Let's instantiate a basis
00105   try {
00106     const int deg = 3;
00107     Basis_HGRAD_TET_Cn_FEM<double,FieldContainer<double> >  myBasis( deg , POINTTYPE_WARPBLEND );
00108 
00109     // Get a lattice
00110     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00111     const int nbf = myBasis.getCardinality();
00112     FieldContainer<double> lattice( np_lattice , 3 );
00113     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00114                                                             myBasis.getBaseCellTopology() , 
00115                                                             deg , 
00116                                                             0 , 
00117                                                             POINTTYPE_WARPBLEND );         
00118     FieldContainer<double> vals( nbf , np_lattice );
00119 
00120     myBasis.getValues( vals , lattice , OPERATOR_VALUE );
00121 
00122     // test for Kronecker property
00123     for (int i=0;i<nbf;i++) {
00124       for (int j=0;j<np_lattice;j++) {
00125         if ( i==j && std::abs( vals(i,j) - 1.0 ) > 100.0 * INTREPID_TOL ) {
00126           errorFlag++;
00127           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00128           *outStream << " Basis function " << i << " does not have unit value at its node\n";
00129         }
00130         if ( i!=j && std::abs( vals(i,j) ) > 100.0 * INTREPID_TOL ) {
00131           errorFlag++;
00132           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00133           *outStream << " Basis function " << i << " does not vanish at node " << j << "\n";
00134           *outStream << " Basis function value is " << vals(i,j) << "\n";
00135         }
00136       }
00137     }
00138   }
00139   catch (std::exception err) {
00140     *outStream << err.what() << "\n\n";
00141     errorFlag = -1000;
00142   }
00143 
00144   try {
00145     const int deg = 4;
00146     Basis_HGRAD_TET_Cn_FEM<double,FieldContainer<double> >  myBasis( deg , POINTTYPE_WARPBLEND );
00147     const std::vector<std::vector<std::vector<int> > >& dofdata = myBasis.getDofOrdinalData();
00148     for (unsigned d=0;d<dofdata.size();d++) {
00149       std::cout << "Dimension " << d << "\n";
00150       for (unsigned f=0;f<dofdata[d].size();f++) {
00151         std::cout << "\tFacet number " << f << "\n";
00152         std::cout << "\t\tDOFS:\n";
00153         for (unsigned n=0;n<dofdata[d][f].size();n++) {
00154           if ( dofdata[d][f][n] >= 0 ) {
00155             std::cout << "\t\t\t" << dofdata[d][f][n] << "\n";
00156           }
00157         }
00158       }
00159     }
00160   }
00161   catch (std::exception err) {
00162     *outStream << err.what() << "\n\n";
00163     errorFlag = -1000;
00164   }
00165 
00166   try {
00167     const int deg = 3;
00168     Basis_HGRAD_TET_Cn_FEM<double,FieldContainer<double> >  myBasis( deg , POINTTYPE_EQUISPACED );
00169 
00170     // Get a lattice
00171     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00172     const int nbf = myBasis.getCardinality();
00173     FieldContainer<double> lattice( np_lattice , 3 );
00174     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00175                                                             myBasis.getBaseCellTopology() , 
00176                                                             deg , 
00177                                                             0 , 
00178                                                             POINTTYPE_EQUISPACED );         
00179     FieldContainer<double> vals( nbf , np_lattice , 3 );
00180 
00181     myBasis.getValues( vals , lattice , OPERATOR_GRAD );
00182 
00183   }
00184   catch (std::exception err) {
00185     *outStream << err.what() << "\n\n";
00186     errorFlag = -1000;
00187   }
00188 
00189   if (errorFlag != 0)
00190     std::cout << "End Result: TEST FAILED\n";
00191   else
00192     std::cout << "End Result: TEST PASSED\n";
00193   
00194   // reset format state of std::cout
00195   std::cout.copyfmt(oldFormatState);
00196   
00197   return errorFlag;
00198 }