Intrepid
example_04.cpp
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
66// Intrepid includes
71#include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
77#include "Intrepid_Utils.hpp"
78
79// Epetra includes
80#include "Epetra_Time.h"
81#include "Epetra_Map.h"
82#include "Epetra_FECrsMatrix.h"
83#include "Epetra_FEVector.h"
84#include "Epetra_SerialComm.h"
85
86// Teuchos includes
87#include "Teuchos_oblackholestream.hpp"
88#include "Teuchos_RCP.hpp"
89#include "Teuchos_BLAS.hpp"
90
91// Shards includes
92#include "Shards_CellTopology.hpp"
93
94// EpetraExt includes
95#include "EpetraExt_RowMatrixOut.h"
96#include "EpetraExt_MultiVectorOut.h"
97
98using namespace std;
99using namespace Intrepid;
100
101int main(int argc, char *argv[]) {
102
103 //Check number of arguments
104 if (argc < 2) {
105 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
106 std::cout <<"Usage:\n\n";
107 std::cout <<" ./Intrepid_example_Drivers_Example_05.exe min_deg max_deg verbose\n\n";
108 std::cout <<" where \n";
109 std::cout <<" int min_deg - beginning polynomial degree to check \n";
110 std::cout <<" int max_deg - final polynomial degree to check \n";
111 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n";
112 exit(1);
113 }
114
115 // This little trick lets us print to std::cout only if
116 // a (dummy) command-line argument is provided.
117 int iprint = argc - 1;
118 Teuchos::RCP<std::ostream> outStream;
119 Teuchos::oblackholestream bhs; // outputs nothing
120 if (iprint > 1)
121 outStream = Teuchos::rcp(&std::cout, false);
122 else
123 outStream = Teuchos::rcp(&bhs, false);
124
125 // Save the format state of the original std::cout.
126 Teuchos::oblackholestream oldFormatState;
127 oldFormatState.copyfmt(std::cout);
128
129 *outStream \
130 << "===============================================================================\n" \
131 << "| |\n" \
132 << "| Example: Check diagonalization of reference mass matrix |\n" \
133 << "| on line, quad, and hex |\n" \
134 << "| |\n" \
135 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
136 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
137 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
138 << "| |\n" \
139 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
140 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
141 << "| |\n" \
142 << "===============================================================================\n";
143
144
145// ************************************ GET INPUTS **************************************
146 int min_degree = atoi(argv[1]);
147 int max_degree = atoi(argv[2]);
148
149
150// *********************************** CELL TOPOLOGY **********************************
151
152 // Get cell topology for base line
153 typedef shards::CellTopology CellTopology;
154 CellTopology line_2(shards::getCellTopologyData<shards::Line<2> >() );
155 CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
156 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
157
158 std::vector<CellTopology> cts(3);
159
160 // loop over degrees
161 for (int deg=min_degree;deg<=max_degree;deg++) {
162 std::vector<Teuchos::RCP<Basis<double,FieldContainer<double> > > > bases;
163 bases.push_back( Teuchos::rcp(new Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
164 bases.push_back( Teuchos::rcp(new Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
165 bases.push_back( Teuchos::rcp(new Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> >(deg,POINTTYPE_SPECTRAL)) );
166
167 // ************************************ CUBATURE **************************************
168 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
169 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
170
171 std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > >
172 cub_to_tensor;
173
174 // now we loop over spatial dimensions
175 for (int sd=1;sd<=3;sd++) {
176 // ************************************ CUBATURE **************************************
177 cub_to_tensor.push_back( glcub );
178 CubatureTensor<double,FieldContainer<double>,FieldContainer<double> > cubcur( cub_to_tensor );
179
180 int cubDim = cubcur.getDimension();
181 int numCubPoints = cubcur.getNumPoints();
182
183 FieldContainer<double> cubPoints(numCubPoints, cubDim);
184 FieldContainer<double> cubWeights(numCubPoints);
185
186 cubcur.getCubature(cubPoints, cubWeights);
187
188 // ************************************ BASIS AT QP**************************************
189 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis_cur = bases[sd-1];
190 const int numFields = basis_cur->getCardinality();
191 FieldContainer<double> bf_at_cub_pts( numFields, numCubPoints );
192 FieldContainer<double> trans_bf_at_cub_pts( 1 , numFields,numCubPoints );
193 FieldContainer<double> w_trans_bf_at_cub_pts( 1, numFields , numCubPoints );
194 basis_cur->getValues( bf_at_cub_pts , cubPoints , OPERATOR_VALUE );
195
196 // *********************************** FORM MASS MATRIX *****************************
197 FunctionSpaceTools::HGRADtransformVALUE<double>( trans_bf_at_cub_pts ,
198 bf_at_cub_pts );
199 cubWeights.resize(1,numCubPoints);
200 FunctionSpaceTools::multiplyMeasure<double>( w_trans_bf_at_cub_pts ,
201 cubWeights ,
202 trans_bf_at_cub_pts );
203 cubWeights.resize(numCubPoints);
204
205 FieldContainer<double> mass_matrix( 1 , numFields, numFields );
206 FunctionSpaceTools::integrate<double>( mass_matrix ,
207 trans_bf_at_cub_pts ,
208 w_trans_bf_at_cub_pts ,
209 COMP_BLAS );
210
211 // now we loop over the mass matrix and compare the nondiagonal entries to zero
212 double max_offdiag = 0.0;
213 for (int i=0;i<numFields;i++) {
214 for (int j=0;j<numFields;j++) {
215 if (i != j) {
216 if ( abs(mass_matrix(0,i,j)) >= max_offdiag) {
217 max_offdiag = abs(mass_matrix(0,i,j));
218 }
219 }
220 }
221 }
222 double min_diag = mass_matrix(0,0,0);
223 for (int i=0;i<numFields;i++) {
224 if ( mass_matrix(0,i,i) <= min_diag ) {
225 min_diag = mass_matrix(0,i,i);
226 }
227 }
228 *outStream << "Degree = " << deg << " and dimension = " << sd << "; Max offdiagonal"
229 << " element is " << max_offdiag << " and min diagonal element is " << min_diag
230 << ".\n";
231
232 }
233
234 }
235
236 std::cout << "End Result: TEST PASSED\n";
237
238 std::cout.copyfmt(oldFormatState);
239
240return 0;
241}
242
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for the Intrepid::CellTools class.
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::CubatureTensor class.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Header file for the Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for the Intrepid::HGRAD_QUAD_Cn_FEM class.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid utilities.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
Utilizes cubature (integration) rules contained in the library Polylib (Spencer Sherwin,...