50 #include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
57 #include "Teuchos_oblackholestream.hpp"
58 #include "Teuchos_RCP.hpp"
59 #include "Teuchos_GlobalMPISession.hpp"
60 #include "Teuchos_SerialDenseMatrix.hpp"
61 #include "Teuchos_SerialDenseVector.hpp"
62 #include "Teuchos_LAPACK.hpp"
65 using namespace Intrepid;
74 for (
int cell=0; cell<result.
dimension(0); cell++) {
75 for (
int pt=0; pt<result.
dimension(1); pt++) {
76 result(cell,pt) = 1.0;
81 else if (degree == 1) {
82 for (
int cell=0; cell<result.
dimension(0); cell++) {
83 for (
int pt=0; pt<result.
dimension(1); pt++) {
84 result(cell,pt) = points(cell,pt,0);
89 for (
int cell=0; cell<result.
dimension(0); cell++) {
90 for (
int pt=0; pt<result.
dimension(1); pt++) {
91 result(cell,pt) = pow(points(cell,pt,0), degree) - degree*(degree-1)*pow(points(cell,pt,0), degree-2);
99 double g_at_one, g_at_minus_one;
104 g_at_minus_one = 0.0;
107 g_at_one = degree*pow(1.0, degree-1);
108 g_at_minus_one = degree*pow(-1.0, degree-1);
113 for (
int i=0; i<num_fields; i++) {
114 g_phi(0,i) = g_at_minus_one*phi1(i,0);
115 g_phi(1,i) = g_at_one*phi2(i,0);
121 for (
int cell=0; cell<result.
dimension(0); cell++) {
122 for (
int pt=0; pt<result.
dimension(1); pt++) {
123 result(cell,pt) = pow(points(pt,0), degree);
131 int main(
int argc,
char *argv[]) {
133 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
137 int iprint = argc - 1;
138 Teuchos::RCP<std::ostream> outStream;
139 Teuchos::oblackholestream bhs;
141 outStream = Teuchos::rcp(&std::cout,
false);
143 outStream = Teuchos::rcp(&bhs,
false);
146 Teuchos::oblackholestream oldFormatState;
147 oldFormatState.copyfmt(std::cout);
150 <<
"===============================================================================\n" \
152 <<
"| Unit Test (Basis_HGRAD_LINE_Cn_FEM) |\n" \
154 <<
"| 1) Patch test involving mass and stiffness matrices, |\n" \
155 <<
"| for the Neumann problem on a REFERENCE line: |\n" \
157 <<
"| - u'' + u = f in (-1,1), u' = g at -1,1 |\n" \
159 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
160 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
161 <<
"| Robert Kirby (robert.c.kirby@ttu.gov), |\n" \
162 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
164 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
165 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
167 <<
"===============================================================================\n"\
168 <<
"| TEST 1: Patch test |\n"\
169 <<
"===============================================================================\n";
173 double zero = 100*INTREPID_TOL;
174 outStream -> precision(20);
182 int numIntervals = 100;
183 int numInterpPoints = numIntervals + 1;
185 for (
int i=0; i<numInterpPoints; i++) {
186 interp_points(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
190 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
193 for (
int basis_order=1; basis_order <= max_order; basis_order++) {
195 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
198 for (
int soln_order=1; soln_order <= basis_order; soln_order++) {
201 u_exact(exact_solution, interp_points, soln_order);
203 int numFields = lineBasis->getCardinality();
206 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.
create(line, 2*basis_order-2);
207 int numCubPoints = lineCub->getNumPoints();
208 int spaceDim = lineCub->getDimension();
225 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
245 lineCub->getCubature(cub_points, cub_weights);
248 cell_nodes(0, 0, 0) = -1.0;
249 cell_nodes(0, 1, 0) = 1.0;
257 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
262 lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
265 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
266 value_of_basis_at_cub_points);
269 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
271 transformed_value_of_basis_at_cub_points);
274 FunctionSpaceTools::integrate<double>(fe_matrix,
275 transformed_value_of_basis_at_cub_points,
276 weighted_transformed_value_of_basis_at_cub_points,
283 lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
286 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
288 grad_of_basis_at_cub_points);
291 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
293 transformed_grad_of_basis_at_cub_points);
296 FunctionSpaceTools::integrate<double>(fe_matrix,
297 transformed_grad_of_basis_at_cub_points,
298 weighted_transformed_grad_of_basis_at_cub_points,
309 rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
312 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
313 rhs_at_cub_points_physical,
314 weighted_transformed_value_of_basis_at_cub_points,
318 one_point(0,0) = 1.0; lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
319 one_point(0,0) = -1.0; lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
320 neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
321 for (
int i=0; i<numFields; i++) {
322 rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
323 rhs_and_soln_vector(0, i) += bc_neumann(1, i);
330 Teuchos::LAPACK<int, double> solver;
332 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
338 lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
340 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
341 value_of_basis_at_interp_points);
342 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
349 *outStream <<
"\nNorm-2 difference between exact solution polynomial of order "
350 << soln_order <<
" and finite element interpolant of order " << basis_order <<
": "
354 *outStream <<
"\n\nPatch test failed for solution polynomial order "
355 << soln_order <<
" and basis order " << basis_order <<
"\n\n";
365 catch (
const std::logic_error & err) {
366 *outStream << err.what() <<
"\n\n";
371 std::cout <<
"End Result: TEST FAILED\n";
373 std::cout <<
"End Result: TEST PASSED\n";
376 std::cout.copyfmt(oldFormatState);
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
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.