50 #include "Intrepid_HGRAD_TRI_Cn_FEM_ORTH.hpp"
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 & ,
84 for (
int cell=0; cell<result.
dimension(0); cell++) {
85 for (
int pt=0; pt<result.
dimension(1); pt++) {
86 result(cell,pt) = - xd*(xd-1)*std::pow(points(cell,pt,x), xd-2) * std::pow(points(cell,pt,y), yd);
93 for (
int cell=0; cell<result.
dimension(0); cell++) {
94 for (
int pt=0; pt<result.
dimension(1); pt++) {
95 result(cell,pt) -= yd*(yd-1)*std::pow(points(cell,pt,y), yd-2) * std::pow(points(cell,pt,x), xd);
101 for (
int cell=0; cell<result.
dimension(0); cell++) {
102 for (
int pt=0; pt<result.
dimension(1); pt++) {
103 result(cell,pt) += std::pow(points(cell,pt,x), xd) * std::pow(points(cell,pt,y), yd);
114 const shards::CellTopology & parentCell,
115 int sideOrdinal,
int xd,
int yd) {
128 for (
int cell=0; cell<numCells; cell++) {
129 for (
int pt=0; pt<numPoints; pt++) {
130 grad_u(cell,pt,x) = xd*std::pow(points(cell,pt,x), xd-1) * std::pow(points(cell,pt,y), yd);
137 for (
int cell=0; cell<numCells; cell++) {
138 for (
int pt=0; pt<numPoints; pt++) {
139 grad_u(cell,pt,y) = yd*std::pow(points(cell,pt,y), yd-1) * std::pow(points(cell,pt,x), xd);
148 FunctionSpaceTools::scalarMultiplyDataData<double>(side_normals, normal_lengths, side_normals,
true);
150 FunctionSpaceTools::dotMultiplyDataData<double>(result, grad_u, side_normals);
157 for (
int cell=0; cell<result.
dimension(0); cell++) {
158 for (
int pt=0; pt<result.
dimension(1); pt++) {
159 result(cell,pt) = std::pow(points(pt,x), xd)*std::pow(points(pt,y), yd);
167 int main(
int argc,
char *argv[]) {
169 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
173 int iprint = argc - 1;
174 Teuchos::RCP<std::ostream> outStream;
175 Teuchos::oblackholestream bhs;
177 outStream = Teuchos::rcp(&std::cout,
false);
179 outStream = Teuchos::rcp(&bhs,
false);
182 Teuchos::oblackholestream oldFormatState;
183 oldFormatState.copyfmt(std::cout);
186 <<
"===============================================================================\n" \
188 <<
"| Unit Test (Basis_HGRAD_TRI_Cn_FEM_ORTH) |\n" \
190 <<
"| 1) Patch test involving mass and stiffness matrices, |\n" \
191 <<
"| for the Neumann problem on a triangular patch |\n" \
192 <<
"| Omega with boundary Gamma. |\n" \
194 <<
"| - div (grad u) + u = f in Omega, (grad u) . n = g on Gamma |\n" \
196 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
197 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
198 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
200 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
201 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
203 <<
"===============================================================================\n"\
204 <<
"| TEST 1: Patch test |\n"\
205 <<
"===============================================================================\n";
210 outStream -> precision(16);
217 shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >());
218 shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >());
219 int cellDim = cell.getDimension();
220 int sideDim = side.getDimension();
223 int numIntervals = 10;
224 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2;
227 for (
int j=0; j<=numIntervals; j++) {
228 for (
int i=0; i<=numIntervals; i++) {
229 if (i <= numIntervals-j) {
230 interp_points_ref(counter,0) = i*(1.0/numIntervals);
231 interp_points_ref(counter,1) = j*(1.0/numIntervals);
240 cell_nodes(0, 0, 0) = 0.1;
241 cell_nodes(0, 0, 1) = -0.1;
242 cell_nodes(0, 1, 0) = 1.1;
243 cell_nodes(0, 1, 1) = -0.1;
244 cell_nodes(0, 2, 0) = 0.1;
245 cell_nodes(0, 2, 1) = 0.9;
256 interp_points.resize(numInterpPoints, cellDim);
258 for (
int x_order=0; x_order <= max_order; x_order++) {
259 for (
int y_order=0; y_order <= max_order-x_order; y_order++) {
263 u_exact(exact_solution, interp_points, x_order, y_order);
265 int total_order = std::max(x_order + y_order, 1);
267 for (
int basis_order=total_order; basis_order <= max_order; basis_order++) {
270 double zero = basis_order*basis_order*100*INTREPID_TOL;
273 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
275 int numFields = basis->getCardinality();
278 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*basis_order);
279 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.
create(side, 2*basis_order);
280 int numCubPointsCell = cellCub->getNumPoints();
281 int numCubPointsSide = sideCub->getNumPoints();
295 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell);
297 FieldContainer<double> transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
298 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim);
305 unsigned numSides = 3;
315 FieldContainer<double> transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
316 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points_side_refcell(1, numFields, numCubPointsSide);
332 cellCub->getCubature(cub_points_cell, cub_weights_cell);
340 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure_cell, jacobian_det_cell, cub_weights_cell);
345 basis->getValues(value_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_VALUE);
348 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_cell,
349 value_of_basis_at_cub_points_cell);
352 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_cell,
353 weighted_measure_cell,
354 transformed_value_of_basis_at_cub_points_cell);
357 FunctionSpaceTools::integrate<double>(fe_matrix,
358 transformed_value_of_basis_at_cub_points_cell,
359 weighted_transformed_value_of_basis_at_cub_points_cell,
366 basis->getValues(grad_of_basis_at_cub_points_cell, cub_points_cell, OPERATOR_GRAD);
369 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points_cell,
371 grad_of_basis_at_cub_points_cell);
374 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points_cell,
375 weighted_measure_cell,
376 transformed_grad_of_basis_at_cub_points_cell);
379 FunctionSpaceTools::integrate<double>(fe_matrix,
380 transformed_grad_of_basis_at_cub_points_cell,
381 weighted_transformed_grad_of_basis_at_cub_points_cell,
392 rhsFunc(rhs_at_cub_points_cell_physical, cub_points_cell_physical, x_order, y_order);
395 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
396 rhs_at_cub_points_cell_physical,
397 weighted_transformed_value_of_basis_at_cub_points_cell,
401 sideCub->getCubature(cub_points_side, cub_weights_side);
402 for (
unsigned i=0; i<numSides; i++) {
409 FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
410 jacobian_side_refcell,
416 basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
418 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
419 value_of_basis_at_cub_points_side_refcell);
422 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
423 weighted_measure_side_refcell,
424 transformed_value_of_basis_at_cub_points_side_refcell);
430 neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
431 cell, (
int)i, x_order, y_order);
433 FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
434 neumann_data_at_cub_points_side_physical,
435 weighted_transformed_value_of_basis_at_cub_points_side_refcell,
446 Teuchos::LAPACK<int, double> solver;
447 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
453 basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
455 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
456 value_of_basis_at_interp_points);
457 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
464 *outStream <<
"\nRelative norm-2 error between exact solution polynomial of order ("
465 << x_order <<
", " << y_order <<
") and finite element interpolant of order " << basis_order <<
": "
471 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
472 << x_order <<
", " << y_order <<
") and basis order " << basis_order <<
"\n\n";
481 catch (
const std::logic_error & err) {
482 *outStream << err.what() <<
"\n\n";
487 std::cout <<
"End Result: TEST FAILED\n";
489 std::cout <<
"End Result: TEST PASSED\n";
492 std::cout.copyfmt(oldFormatState);
Implementation of the default H(grad)-compatible orthogonal basis (Dubiner) of arbitrary degree on tr...
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.