58 #include "Teuchos_oblackholestream.hpp"
59 #include "Teuchos_RCP.hpp"
60 #include "Teuchos_GlobalMPISession.hpp"
61 #include "Teuchos_SerialDenseMatrix.hpp"
62 #include "Teuchos_SerialDenseVector.hpp"
63 #include "Teuchos_LAPACK.hpp"
66 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) = 0.0;
83 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd)
84 *pow(points(cell,pt,2),zd);
87 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2)
88 *pow(points(cell,pt,2),zd);
91 result(cell,pt) += zd*(zd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd)
92 *pow(points(cell,pt,2),zd-2);
105 for (
int cell=0;cell<result.
dimension(0);cell++){
106 for (
int pt=0;pt<result.
dimension(1);pt++) {
107 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
108 *std::pow(points(cell,pt,2),zd);
114 int main(
int argc,
char *argv[]) {
115 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
119 int iprint = argc - 1;
120 Teuchos::RCP<std::ostream> outStream;
121 Teuchos::oblackholestream bhs;
123 outStream = Teuchos::rcp(&std::cout,
false);
125 outStream = Teuchos::rcp(&bhs,
false);
128 Teuchos::oblackholestream oldFormatState;
129 oldFormatState.copyfmt(std::cout);
132 <<
"===============================================================================\n" \
134 <<
"| Unit Test (Basis_HGRAD_HEX_In_FEM) |\n" \
136 <<
"| 1) Patch test involving H(div) matrices |\n" \
137 <<
"| for the Dirichlet problem on a hexahedron |\n" \
138 <<
"| Omega with boundary Gamma. |\n" \
140 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
141 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
142 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
143 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
145 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
146 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
148 <<
"===============================================================================\n" \
149 <<
"| TEST 1: Patch test |\n" \
150 <<
"===============================================================================\n";
155 outStream -> precision(16);
159 shards::CellTopology cell(shards::getCellTopologyData< shards::Hexahedron<> >());
160 shards::CellTopology side(shards::getCellTopologyData< shards::Quadrilateral<> >());
161 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >() );
163 int cellDim = cell.getDimension();
164 int sideDim = side.getDimension();
169 int numIntervals = 2;
170 int numInterpPoints = (numIntervals + 1)*(numIntervals + 1)*(numIntervals+1);
173 for (
int k=0;k<numIntervals;k++) {
174 for (
int j=0; j<=numIntervals; j++) {
175 for (
int i=0; i<=numIntervals; i++) {
176 interp_points_ref(counter,0) = i*(1.0/numIntervals);
177 interp_points_ref(counter,1) = j*(1.0/numIntervals);
178 interp_points_ref(counter,2) = k*(1.0/numIntervals);
184 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
187 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
190 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
193 int numVectorFields = vectorBasis->getCardinality();
194 int numScalarFields = scalarBasis->getCardinality();
195 int numTotalFields = numVectorFields + numScalarFields;
198 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*(basis_order+1));
199 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.
create(side, 2*(basis_order+1));
201 int numCubPointsCell = cellCub->getNumPoints();
202 int numCubPointsSide = sideCub->getNumPoints();
213 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
245 double zero = (basis_order+1)*(basis_order+1)*1000.0*INTREPID_TOL;
250 cellCub->getCubature(cub_points_cell, cub_weights_cell);
251 sideCub->getCubature(cub_points_side, cub_weights_side);
254 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
257 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
262 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
267 cub_weights_cell.resize(1,numCubPointsCell);
268 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
270 value_of_v_basis_at_cub_points_cell );
271 cub_weights_cell.resize(numCubPointsCell);
274 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
275 FunctionSpaceTools::integrate<double>(fe_matrix_M,
276 w_value_of_v_basis_at_cub_points_cell ,
277 value_of_v_basis_at_cub_points_cell ,
279 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
282 cub_weights_cell.resize(1,numCubPointsCell);
283 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
285 div_of_v_basis_at_cub_points_cell);
286 cub_weights_cell.resize(numCubPointsCell);
288 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
289 FunctionSpaceTools::integrate<double>(fe_matrix_B,
290 w_div_of_v_basis_at_cub_points_cell ,
291 value_of_s_basis_at_cub_points_cell ,
293 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
295 for (
int x_order=0;x_order<=basis_order;x_order++) {
296 for (
int y_order=0;y_order<=basis_order;y_order++) {
297 for (
int z_order=0;z_order<=basis_order;z_order++) {
301 fe_matrix.initialize();
303 for (
int i=0;i<numVectorFields;i++) {
304 for (
int j=0;j<numVectorFields;j++) {
305 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
310 for (
int i=0;i<numVectorFields;i++) {
311 for (
int j=0;j<numScalarFields;j++) {
312 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
313 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
318 rhs_vector_vec.initialize();
319 rhs_vector_scal.initialize();
320 rhs_and_soln_vec.initialize();
325 cub_points_cell.resize(1,numCubPointsCell,cellDim);
326 rhsFunc(rhs_at_cub_points_cell,
332 cub_points_cell.resize(numCubPointsCell,cellDim);
334 cub_weights_cell.resize(1,numCubPointsCell);
335 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
337 value_of_s_basis_at_cub_points_cell);
338 cub_weights_cell.resize(numCubPointsCell);
339 FunctionSpaceTools::integrate<double>(rhs_vector_scal,
340 rhs_at_cub_points_cell,
341 w_value_of_s_basis_at_cub_points_cell,
344 for (
int i=0;i<numScalarFields;i++) {
345 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
350 for (
unsigned side_cur=0;side_cur<6;side_cur++) {
358 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
359 u_exact(diri_data_at_cub_points_side,
360 cub_points_side_refcell,x_order,y_order,z_order);
362 cub_points_side_refcell.resize(numCubPointsSide,cellDim);
366 (
int)side_cur,cell );
369 vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
370 cub_points_side_refcell ,
373 for (
int i=0;i<numVectorFields;i++) {
374 for (
int j=0;j<numCubPointsSide;j++) {
375 n_of_v_basis_at_cub_points_side(i,j) = 0.0;
376 for (
int k=0;k<cellDim;k++) {
377 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) *
378 value_of_v_basis_at_cub_points_side(i,j,k);
383 cub_weights_side.resize(1,numCubPointsSide);
384 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
386 n_of_v_basis_at_cub_points_side);
387 cub_weights_side.resize(numCubPointsSide);
389 FunctionSpaceTools::integrate<double>(rhs_vector_vec,
390 diri_data_at_cub_points_side,
391 w_n_of_v_basis_at_cub_points_side,
395 for (
int i=0;i<numVectorFields;i++) {
396 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
403 Teuchos::LAPACK<int, double> solver;
404 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0],
405 numTotalFields, &info);
408 scalarBasis->getValues(value_of_s_basis_at_interp_points,
411 for (
int pt=0;pt<numInterpPoints;pt++) {
412 interpolant(0,pt)=0.0;
413 for (
int i=0;i<numScalarFields;i++) {
414 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
415 * value_of_s_basis_at_interp_points(i,pt);
419 interp_points_ref.resize(1,numInterpPoints,cellDim);
422 u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
423 interp_points_ref.resize(numInterpPoints,cellDim);
429 *outStream <<
"\nNorm-2 error between scalar components of exact solution of order ("
430 << x_order <<
", " << y_order <<
", " << z_order
431 <<
") and finite element interpolant of order " << basis_order <<
": "
435 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
436 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order (scalar, vector) ("
437 << basis_order <<
", " << basis_order+1 <<
")\n\n";
448 catch (
const std::logic_error & err) {
449 *outStream << err.what() <<
"\n\n";
454 std::cout <<
"End Result: TEST FAILED\n";
456 std::cout <<
"End Result: TEST PASSED\n";
459 std::cout.copyfmt(oldFormatState);
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Header file for the Intrepid::HDIV_HEX_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 Intrepid::HGRAD_HEX_Cn_FEM class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell...
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.