|
Intrepid
|
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
1.7.6.1