59#include "Teuchos_oblackholestream.hpp"
60#include "Teuchos_RCP.hpp"
61#include "Teuchos_GlobalMPISession.hpp"
62#include "Teuchos_SerialDenseMatrix.hpp"
63#include "Teuchos_SerialDenseVector.hpp"
64#include "Teuchos_LAPACK.hpp"
67using namespace Intrepid;
69void rhsFunc( FieldContainer<double> &,
const FieldContainer<double> &,
int,
int,
int,
int );
70void u_exact( FieldContainer<double> &,
const FieldContainer<double> &,
int ,
int,
int,
int );
72void u_exact( FieldContainer<double> &result,
73 const FieldContainer<double> &points,
79 for (
int cell=0;cell<result.dimension(0);cell++){
80 for (
int pt=0;pt<result.dimension(1);pt++) {
81 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
82 *std::pow(points(cell,pt,2),zd);
88void rhsFunc( FieldContainer<double> & result ,
89 const FieldContainer<double> &points ,
95 u_exact( result , points , comp , xd , yd , zd );
98int main(
int argc,
char *argv[]) {
99 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
103 int iprint = argc - 1;
104 Teuchos::RCP<std::ostream> outStream;
105 Teuchos::oblackholestream bhs;
107 outStream = Teuchos::rcp(&std::cout,
false);
109 outStream = Teuchos::rcp(&bhs,
false);
112 Teuchos::oblackholestream oldFormatState;
113 oldFormatState.copyfmt(std::cout);
116 <<
"===============================================================================\n" \
118 <<
"| Unit Test (Basis_HCURL_HEX_In_FEM) |\n" \
120 <<
"| 1) Patch test involving H(curl) matrices |\n" \
122 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
123 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
124 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
125 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
127 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
128 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
130 <<
"===============================================================================\n" \
131 <<
"| TEST 2: Patch test for mass matrices |\n" \
132 <<
"===============================================================================\n";
137 outStream -> precision(16);
140 shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<8> >());
142 int cellDim = cell.getDimension();
147 int numIntervals = max_order;
148 int numInterpPoints = (numIntervals + 1)*(numIntervals +1)*(numIntervals+1);
149 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
151 for (
int k=0; k<=numIntervals; k++) {
152 for (
int j=0; j<=numIntervals; j++) {
153 for (
int i=0;i<numIntervals;i++) {
154 interp_points_ref(counter,0) = i*(1.0/numIntervals);
155 interp_points_ref(counter,1) = j*(1.0/numIntervals);
156 interp_points_ref(counter,2) = k*(1.0/numIntervals);
162 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
164 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
165 Teuchos::rcp(
new Basis_HCURL_HEX_In_FEM<
double,FieldContainer<double> >(basis_order,POINTTYPE_SPECTRAL) );
166 int numFields = basis->getCardinality();
169 DefaultCubatureFactory<double> cubFactory;
170 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create( cell , 2* basis_order );
173 int numCubPointsCell = cellCub->getNumPoints();
176 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
177 FieldContainer<double> cub_weights_cell(numCubPointsCell);
180 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim );
181 FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
185 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim);
188 FieldContainer<double> fe_matrix_bak(1,numFields,numFields);
189 FieldContainer<double> fe_matrix(1,numFields,numFields);
190 FieldContainer<double> rhs_and_soln_vec(1,numFields);
192 FieldContainer<int> ipiv(numFields);
193 FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim);
194 FieldContainer<double> interpolant( 1, numInterpPoints , cellDim );
197 Teuchos::LAPACK<int, double> solver;
200 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
204 cellCub->getCubature(cub_points_cell, cub_weights_cell);
207 basis->getValues(value_of_basis_at_cub_points_cell,
210 basis->getValues( value_of_basis_at_interp_points ,
216 cub_weights_cell.resize(1,numCubPointsCell);
217 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
219 value_of_basis_at_cub_points_cell );
220 cub_weights_cell.resize(numCubPointsCell);
223 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
224 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
225 w_value_of_basis_at_cub_points_cell ,
226 value_of_basis_at_cub_points_cell ,
228 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
230 for (
int x_order=0;x_order<basis_order;x_order++) {
231 for (
int y_order=0;y_order<basis_order;y_order++) {
232 for (
int z_order=0;z_order<basis_order;z_order++) {
233 for (
int comp=0;comp<cellDim;comp++) {
234 fe_matrix.initialize();
236 for (
int i=0;i<numFields;i++) {
237 for (
int j=0;j<numFields;j++) {
238 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
243 rhs_and_soln_vec.initialize();
247 cub_points_cell.resize(1,numCubPointsCell,cellDim);
249 rhs_at_cub_points_cell.initialize();
250 rhsFunc(rhs_at_cub_points_cell,
257 cub_points_cell.resize(numCubPointsCell,cellDim);
258 cub_weights_cell.resize(numCubPointsCell);
260 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
261 rhs_at_cub_points_cell,
262 w_value_of_basis_at_cub_points_cell,
267 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vec[0],
272 interp_points_ref.resize(1,numInterpPoints,cellDim);
274 FieldContainer<double> exact_solution(1,numInterpPoints,cellDim);
275 exact_solution.initialize();
276 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
277 interp_points_ref.resize(numInterpPoints,cellDim);
281 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
282 FunctionSpaceTools::evaluate<double>( interpolant ,
284 value_of_basis_at_interp_points );
285 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
290 interpolant.resize(1,numInterpPoints,cellDim);
291 FieldContainer<double> l2norm(1);
292 FunctionSpaceTools::dataIntegral<double>( l2norm ,
297 const double nrm = sqrtl(l2norm(0));
299 *outStream <<
"\nFunction space L^2 Norm-2 error between scalar components of exact solution of order ("
300 << x_order <<
", " << y_order <<
", " << z_order
301 <<
") in component " << comp
302 <<
" and finite element interpolant of order " << basis_order <<
": "
306 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
307 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order (scalar, vector) ("
308 << basis_order <<
", " << basis_order+1 <<
")\n\n";
319 catch (
const std::logic_error & err) {
320 *outStream << err.what() <<
"\n\n";
325 std::cout <<
"End Result: TEST FAILED\n";
327 std::cout <<
"End Result: TEST PASSED\n";
330 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::CubatureTensor class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_HEX_In_FEM class.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell.