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;
70 const shards::CellTopology & ,
81 int x = 0, y = 1, z = 2;
85 for (
int cell=0; cell<result.
dimension(0); cell++) {
86 for (
int pt=0; pt<result.
dimension(1); pt++) {
87 result(cell,pt) = - xd*(xd-1)*std::pow(points(cell,pt,x), xd-2) *
88 std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
95 for (
int cell=0; cell<result.
dimension(0); cell++) {
96 for (
int pt=0; pt<result.
dimension(1); pt++) {
97 result(cell,pt) -= yd*(yd-1)*std::pow(points(cell,pt,y), yd-2) *
98 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,z), zd);
105 for (
int cell=0; cell<result.
dimension(0); cell++) {
106 for (
int pt=0; pt<result.
dimension(1); pt++) {
107 result(cell,pt) -= zd*(zd-1)*std::pow(points(cell,pt,z), zd-2) *
108 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
114 for (
int cell=0; cell<result.
dimension(0); cell++) {
115 for (
int pt=0; pt<result.
dimension(1); pt++) {
116 result(cell,pt) += std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
127 const shards::CellTopology & parentCell,
128 int sideOrdinal,
int xd,
int yd,
int zd) {
130 int x = 0, y = 1, z = 2;
141 for (
int cell=0; cell<numCells; cell++) {
142 for (
int pt=0; pt<numPoints; pt++) {
143 grad_u(cell,pt,x) = xd*std::pow(points(cell,pt,x), xd-1) *
144 std::pow(points(cell,pt,y), yd) * std::pow(points(cell,pt,z), zd);
151 for (
int cell=0; cell<numCells; cell++) {
152 for (
int pt=0; pt<numPoints; pt++) {
153 grad_u(cell,pt,y) = yd*std::pow(points(cell,pt,y), yd-1) *
154 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,z), zd);
161 for (
int cell=0; cell<numCells; cell++) {
162 for (
int pt=0; pt<numPoints; pt++) {
163 grad_u(cell,pt,z) = zd*std::pow(points(cell,pt,z), zd-1) *
164 std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
173 FunctionSpaceTools::scalarMultiplyDataData<double>(side_normals, normal_lengths, side_normals,
true);
175 FunctionSpaceTools::dotMultiplyDataData<double>(result, grad_u, side_normals);
181 int x = 0, y = 1, z = 2;
182 for (
int cell=0; cell<result.
dimension(0); cell++) {
183 for (
int pt=0; pt<result.
dimension(1); pt++) {
184 result(cell,pt) = std::pow(points(pt,x), xd)*std::pow(points(pt,y), yd)*std::pow(points(pt,z), zd);
192 int main(
int argc,
char *argv[]) {
194 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
198 int iprint = argc - 1;
199 Teuchos::RCP<std::ostream> outStream;
200 Teuchos::oblackholestream bhs;
202 outStream = Teuchos::rcp(&std::cout,
false);
204 outStream = Teuchos::rcp(&bhs,
false);
207 Teuchos::oblackholestream oldFormatState;
208 oldFormatState.copyfmt(std::cout);
211 <<
"===============================================================================\n" \
213 <<
"| Unit Test (Basis_HGRAD_TET_Cn_FEM) |\n" \
215 <<
"| 1) Patch test involving mass and stiffness matrices, |\n" \
216 <<
"| for the Neumann problem on a tetrahedral patch |\n" \
217 <<
"| Omega with boundary Gamma. |\n" \
219 <<
"| - div (grad u) + u = f in Omega, (grad u) . n = g on Gamma |\n" \
221 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
222 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
223 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
225 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
226 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
228 <<
"===============================================================================\n"\
229 <<
"| TEST 1: Patch test |\n"\
230 <<
"===============================================================================\n";
235 outStream -> precision(16);
242 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());
243 shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());
244 int cellDim = cell.getDimension();
245 int sideDim = side.getDimension();
248 int numIntervals = 10;
249 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals + 3))/6;
252 for (
int k=0; k<=numIntervals; k++) {
253 for (
int j=0; j<=numIntervals; j++) {
254 for (
int i=0; i<=numIntervals; i++) {
255 if (i+j+k <= numIntervals) {
256 interp_points_ref(counter,0) = i*(1.0/numIntervals);
257 interp_points_ref(counter,1) = j*(1.0/numIntervals);
258 interp_points_ref(counter,2) = k*(1.0/numIntervals);
268 cell_nodes(0, 0, 0) = -1.0;
269 cell_nodes(0, 0, 1) = -2.0;
270 cell_nodes(0, 0, 2) = 0.0;
271 cell_nodes(0, 1, 0) = 6.0;
272 cell_nodes(0, 1, 1) = 2.0;
273 cell_nodes(0, 1, 2) = 0.0;
274 cell_nodes(0, 2, 0) = -5.0;
275 cell_nodes(0, 2, 1) = 1.0;
276 cell_nodes(0, 2, 2) = 0.0;
277 cell_nodes(0, 3, 0) = -4.0;
278 cell_nodes(0, 3, 1) = -1.0;
279 cell_nodes(0, 3, 2) = 3.0;
309 interp_points.resize(numInterpPoints, cellDim);
312 EPointType pointtype[] = {POINTTYPE_EQUISPACED, POINTTYPE_WARPBLEND};
313 for (
int ptype=0; ptype < 2; ptype++) {
315 *outStream <<
"\nTesting bases with " << EPointTypeToString(pointtype[ptype]) <<
":\n";
317 for (
int x_order=0; x_order <= max_order; x_order++) {
318 for (
int y_order=0; y_order <= max_order-x_order; y_order++) {
319 for (
int z_order=0; z_order <= max_order-x_order-y_order; z_order++) {
323 u_exact(exact_solution, interp_points, x_order, y_order, z_order);
325 int total_order = std::max(x_order + y_order + z_order, 1);
327 for (
int basis_order=total_order; basis_order <= max_order; basis_order++) {
330 double zero = basis_order*basis_order*basis_order*100*INTREPID_TOL;
333 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
335 int numFields = basis->getCardinality();
338 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*basis_order);
339 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.
create(side, 2*basis_order);
340 int numCubPointsCell = cellCub->getNumPoints();
341 int numCubPointsSide = sideCub->getNumPoints();
355 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
357 FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
358 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
365 unsigned numSides = 4;
375 FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
376 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
392 cellCub->getCubature(cub_points_cell, cub_weights_cell);
400 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
405 basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
408 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
409 value_of_basis_at_cub_points_cell);
412 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
413 weighted_measure_cell,
414 transformed_value_of_basis_at_cub_points_cell);
417 FunctionSpaceTools::integrate<double>(fe_matrix,
418 transformed_value_of_basis_at_cub_points_cell,
419 weighted_transformed_value_of_basis_at_cub_points_cell,
426 basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
429 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
431 grad_of_basis_at_cub_points_cell);
434 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
435 weighted_measure_cell,
436 transformed_grad_of_basis_at_cub_points_cell);
439 FunctionSpaceTools::integrate<double>(fe_matrix,
440 transformed_grad_of_basis_at_cub_points_cell,
441 weighted_transformed_grad_of_basis_at_cub_points_cell,
452 rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order, z_order);
455 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
456 rhs_at_cub_points_cell_physical,
457 weighted_transformed_value_of_basis_at_cub_points_cell,
461 sideCub->getCubature(cub_points_side, cub_weights_side);
462 for (
unsigned i=0; i<numSides; i++) {
469 FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_side_refcell,
470 jacobian_side_refcell,
476 basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
478 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
479 value_of_basis_at_cub_points_side_refcell);
482 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
483 weighted_measure_side_refcell,
484 transformed_value_of_basis_at_cub_points_side_refcell);
490 neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
491 cell, (
int)i, x_order, y_order, z_order);
493 FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
494 neumann_data_at_cub_points_side_physical,
495 weighted_transformed_value_of_basis_at_cub_points_side_refcell,
506 Teuchos::LAPACK<int, double> solver;
507 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
513 basis->getValues(value_of_basis_at_interp_points_ref, interp_points_ref, OPERATOR_VALUE);
515 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points_ref,
516 value_of_basis_at_interp_points_ref);
517 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points_ref);
524 *outStream <<
"\nRelative norm-2 error between exact solution polynomial of order ("
525 << x_order <<
", " << y_order <<
", " << z_order
526 <<
") and finite element interpolant of order " << basis_order <<
": "
532 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
533 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order " << basis_order <<
"\n\n";
544 catch (
const std::logic_error & err) {
545 *outStream << err.what() <<
"\n\n";
550 std::cout <<
"End Result: TEST FAILED\n";
552 std::cout <<
"End Result: TEST PASSED\n";
555 std::cout.copyfmt(oldFormatState);
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Tetrahedron ce...
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.
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.
Header file for the Intrepid::HGRAD_TET_Cn_FEM class.