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"
67 using namespace Intrepid;
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);
95 u_exact( result , points , comp , xd , yd , zd );
98 int 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);
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 =
166 int numFields = basis->getCardinality();
170 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create( cell , 2* basis_order );
173 int numCubPointsCell = cellCub->getNumPoints();
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);
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);
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);
int dimension(const int whichDim) const
Returns the specified dimension.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::CubatureTensor class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
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.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Header file for the Intrepid::CubaturePolylib class.
Header file for the Intrepid::HCURL_HEX_In_FEM class.