56 #include "Teuchos_oblackholestream.hpp"
57 #include "Teuchos_RCP.hpp"
58 #include "Teuchos_GlobalMPISession.hpp"
59 #include "Teuchos_SerialDenseMatrix.hpp"
60 #include "Teuchos_SerialDenseVector.hpp"
61 #include "Teuchos_LAPACK.hpp"
64 using namespace Intrepid;
76 for (
int cell=0;cell<result.
dimension(0);cell++){
77 for (
int pt=0;pt<result.
dimension(1);pt++) {
78 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
79 *std::pow(points(cell,pt,2),zd);
92 u_exact( result , points , comp , xd , yd , zd );
95 int main(
int argc,
char *argv[]) {
96 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
100 int iprint = argc - 1;
101 Teuchos::RCP<std::ostream> outStream;
102 Teuchos::oblackholestream bhs;
104 outStream = Teuchos::rcp(&std::cout,
false);
106 outStream = Teuchos::rcp(&bhs,
false);
109 Teuchos::oblackholestream oldFormatState;
110 oldFormatState.copyfmt(std::cout);
113 <<
"===============================================================================\n" \
115 <<
"| Unit Test (Basis_HCURL_TET_In_FEM) |\n" \
117 <<
"| 1) Patch test involving H(curl) matrices |\n" \
119 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
120 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
121 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
122 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
124 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
125 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
127 <<
"===============================================================================\n" \
128 <<
"| TEST 2: Patch test for mass matrices |\n" \
129 <<
"===============================================================================\n";
134 outStream -> precision(16);
138 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());
140 int cellDim = cell.getDimension();
145 int numIntervals = max_order;
146 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
149 for (
int j=0; j<=numIntervals; j++) {
150 for (
int i=0; i<=numIntervals-j; i++) {
151 for (
int k=0;k<numIntervals-j-i;k++) {
152 interp_points_ref(counter,0) = i*(1.0/numIntervals);
153 interp_points_ref(counter,1) = j*(1.0/numIntervals);
154 interp_points_ref(counter,2) = k*(1.0/numIntervals);
160 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
162 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
165 int numFields = basis->getCardinality();
168 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*(basis_order+1));
170 int numCubPointsCell = cellCub->getNumPoints();
193 Teuchos::LAPACK<int, double> solver;
196 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
200 cellCub->getCubature(cub_points_cell, cub_weights_cell);
203 basis->getValues(value_of_basis_at_cub_points_cell,
206 basis->getValues( value_of_basis_at_interp_points ,
214 cub_weights_cell.resize(1,numCubPointsCell);
215 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
217 value_of_basis_at_cub_points_cell );
218 cub_weights_cell.resize(numCubPointsCell);
221 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
222 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
223 w_value_of_basis_at_cub_points_cell ,
224 value_of_basis_at_cub_points_cell ,
226 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
231 for (
int x_order=0;x_order<basis_order;x_order++) {
232 for (
int y_order=0;y_order<basis_order-x_order;y_order++) {
233 for (
int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
234 for (
int comp=0;comp<cellDim;comp++) {
235 fe_matrix.initialize();
237 for (
int i=0;i<numFields;i++) {
238 for (
int j=0;j<numFields;j++) {
239 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
244 rhs_and_soln_vec.initialize();
248 cub_points_cell.resize(1,numCubPointsCell,cellDim);
250 rhs_at_cub_points_cell.initialize();
251 rhsFunc(rhs_at_cub_points_cell,
258 cub_points_cell.resize(numCubPointsCell,cellDim);
260 cub_weights_cell.resize(numCubPointsCell);
262 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
263 rhs_at_cub_points_cell,
264 w_value_of_basis_at_cub_points_cell,
271 solver.POTRF(
'L',numFields,&fe_matrix[0],numFields,&info);
272 solver.POTRS(
'L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
274 interp_points_ref.resize(1,numInterpPoints,cellDim);
277 exact_solution.initialize();
278 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
279 interp_points_ref.resize(numInterpPoints,cellDim);
283 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
284 FunctionSpaceTools::evaluate<double>( interpolant ,
286 value_of_basis_at_interp_points );
287 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
293 *outStream <<
"\nNorm-2 error between scalar components of exact solution of order ("
294 << x_order <<
", " << y_order <<
", " << z_order
295 <<
") in component " << comp
296 <<
" and finite element interpolant of order " << basis_order <<
": "
300 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
301 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order (scalar, vector) ("
302 << basis_order <<
", " << basis_order+1 <<
")\n\n";
313 catch (
const std::logic_error & err) {
314 *outStream << err.what() <<
"\n\n";
319 std::cout <<
"End Result: TEST FAILED\n";
321 std::cout <<
"End Result: TEST PASSED\n";
324 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::HCURL_TET_In_FEM class.
int dimension(const int whichDim) const
Returns the specified dimension.
Header file for utility class to provide multidimensional containers.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
A factory class that generates specific instances of cubatures.
Teuchos::RCP< Cubature< Scalar, ArrayPoint, ArrayWeight > > create(const shards::CellTopology &cellTopology, const std::vector< int > °ree)
Factory method.