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;
73 for (
int cell=0; cell<result.
dimension(0); cell++) {
74 for (
int pt=0; pt<result.
dimension(1); pt++) {
75 result(cell,pt) = 1.0;
79 else if (degree == 1) {
80 for (
int cell=0; cell<result.
dimension(0); cell++) {
81 for (
int pt=0; pt<result.
dimension(1); pt++) {
82 result(cell,pt) = points(cell,pt,0);
87 for (
int cell=0; cell<result.
dimension(0); cell++) {
88 for (
int pt=0; pt<result.
dimension(1); pt++) {
89 result(cell,pt) = pow(points(cell,pt,0), degree) - degree*(degree-1)*pow(points(cell,pt,0), degree-2);
97 double g_at_one, g_at_minus_one;
102 g_at_minus_one = 0.0;
105 g_at_one = degree*pow(1.0, degree-1);
106 g_at_minus_one = degree*pow(-1.0, degree-1);
111 for (
int i=0; i<num_fields; i++) {
112 g_phi(0,i) = g_at_minus_one*phi1(i,0);
113 g_phi(1,i) = g_at_one*phi2(i,0);
119 for (
int cell=0; cell<result.
dimension(0); cell++) {
120 for (
int pt=0; pt<result.
dimension(1); pt++) {
121 result(cell,pt) = pow(points(pt,0), degree);
129 int main(
int argc,
char *argv[]) {
131 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
135 int iprint = argc - 1;
136 Teuchos::RCP<std::ostream> outStream;
137 Teuchos::oblackholestream bhs;
139 outStream = Teuchos::rcp(&std::cout,
false);
141 outStream = Teuchos::rcp(&bhs,
false);
144 Teuchos::oblackholestream oldFormatState;
145 oldFormatState.copyfmt(std::cout);
148 <<
"===============================================================================\n" \
150 <<
"| Unit Test (Basis_HGRAD_LINE_Cn_FEM_JACOBI) |\n" \
152 <<
"| 1) Patch test involving mass and stiffness matrices, |\n" \
153 <<
"| for the Neumann problem on a REFERENCE line: |\n" \
155 <<
"| - u'' + u = f in (-1,1), u' = g at -1,1 |\n" \
157 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
158 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
159 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
161 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
162 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
164 <<
"===============================================================================\n"\
165 <<
"| TEST 1: Patch test |\n"\
166 <<
"===============================================================================\n";
170 double zero = 100*INTREPID_TOL;
171 outStream -> precision(20);
179 int numIntervals = 100;
180 int numInterpPoints = numIntervals + 1;
182 for (
int i=0; i<numInterpPoints; i++) {
183 interp_points(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
187 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
189 for (
int soln_order=1; soln_order <= max_order; soln_order++) {
193 u_exact(exact_solution, interp_points, soln_order);
195 for (
int basis_order=soln_order; basis_order <= max_order; basis_order++) {
198 Teuchos::RCP<Basis<double,FieldContainer<double> > > lineBasis =
200 int numFields = lineBasis->getCardinality();
203 Teuchos::RCP<Cubature<double> > lineCub = cubFactory.
create(line, 2*basis_order-2);
204 int numCubPoints = lineCub->getNumPoints();
205 int spaceDim = lineCub->getDimension();
222 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(1, numFields, numCubPoints, spaceDim);
242 lineCub->getCubature(cub_points, cub_weights);
245 cell_nodes(0, 0, 0) = -1.0;
246 cell_nodes(0, 1, 0) = 1.0;
254 FunctionSpaceTools::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
259 lineBasis->getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
262 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
263 value_of_basis_at_cub_points);
266 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
268 transformed_value_of_basis_at_cub_points);
271 FunctionSpaceTools::integrate<double>(fe_matrix,
272 transformed_value_of_basis_at_cub_points,
273 weighted_transformed_value_of_basis_at_cub_points,
280 lineBasis->getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
283 FunctionSpaceTools::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
285 grad_of_basis_at_cub_points);
288 FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
290 transformed_grad_of_basis_at_cub_points);
293 FunctionSpaceTools::integrate<double>(fe_matrix,
294 transformed_grad_of_basis_at_cub_points,
295 weighted_transformed_grad_of_basis_at_cub_points,
306 rhsFunc(rhs_at_cub_points_physical, cub_points_physical, soln_order);
309 FunctionSpaceTools::integrate<double>(rhs_and_soln_vector,
310 rhs_at_cub_points_physical,
311 weighted_transformed_value_of_basis_at_cub_points,
315 one_point(0,0) = 1.0; lineBasis->getValues(value_of_basis_at_one, one_point, OPERATOR_VALUE);
316 one_point(0,0) = -1.0; lineBasis->getValues(value_of_basis_at_minusone, one_point, OPERATOR_VALUE);
317 neumann(bc_neumann, value_of_basis_at_minusone, value_of_basis_at_one, soln_order);
318 for (
int i=0; i<numFields; i++) {
319 rhs_and_soln_vector(0, i) -= bc_neumann(0, i);
320 rhs_and_soln_vector(0, i) += bc_neumann(1, i);
327 Teuchos::LAPACK<int, double> solver;
329 solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
335 lineBasis->getValues(value_of_basis_at_interp_points, interp_points, OPERATOR_VALUE);
337 FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
338 value_of_basis_at_interp_points);
339 FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
346 *outStream <<
"\nNorm-2 difference between exact solution polynomial of order "
347 << soln_order <<
" and finite element interpolant of order " << basis_order <<
": "
351 *outStream <<
"\n\nPatch test failed for solution polynomial order "
352 << soln_order <<
" and basis order " << basis_order <<
"\n\n";
362 catch (
const std::logic_error & err) {
363 *outStream << err.what() <<
"\n\n";
368 std::cout <<
"End Result: TEST FAILED\n";
370 std::cout <<
"End Result: TEST PASSED\n";
373 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::HGRAD_LINE_Cn_FEM class.
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.