50 #include "Intrepid_HGRAD_TRI_Cn_FEM_ORTH.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;
77 for (
int cell=0;cell<result.
dimension(0);cell++) {
78 for (
int pt=0;pt<result.
dimension(1);pt++) {
79 result(cell,pt) = 0.0;
81 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd);
84 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2);
95 for (
int cell=0;cell<result.
dimension(0);cell++){
96 for (
int pt=0;pt<result.
dimension(1);pt++) {
97 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd);
103 int main(
int argc,
char *argv[]) {
104 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
108 int iprint = argc - 1;
109 Teuchos::RCP<std::ostream> outStream;
110 Teuchos::oblackholestream bhs;
112 outStream = Teuchos::rcp(&std::cout,
false);
114 outStream = Teuchos::rcp(&bhs,
false);
117 Teuchos::oblackholestream oldFormatState;
118 oldFormatState.copyfmt(std::cout);
121 <<
"===============================================================================\n" \
123 <<
"| Unit Test (Basis_HGRAD_TRI_Cn_FEM_ORTH) |\n" \
125 <<
"| 1) Patch test involving H(div) matrices |\n" \
126 <<
"| for the Dirichlet problem on a triangular patch |\n" \
127 <<
"| Omega with boundary Gamma. |\n" \
129 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
130 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
131 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
132 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
134 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
135 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
137 <<
"===============================================================================\n"\
138 <<
"| TEST 1: Patch test |\n"\
139 <<
"===============================================================================\n";
144 outStream -> precision(16);
148 shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >());
149 shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >());
151 int cellDim = cell.getDimension();
152 int sideDim = side.getDimension();
157 int numIntervals = max_order;
158 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2;
161 for (
int j=0; j<=numIntervals; j++) {
162 for (
int i=0; i<=numIntervals; i++) {
163 if (i <= numIntervals-j) {
164 interp_points_ref(counter,0) = i*(1.0/numIntervals);
165 interp_points_ref(counter,1) = j*(1.0/numIntervals);
171 interp_points_ref.resize(numInterpPoints,2);
173 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
175 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
177 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
180 int numVectorFields = vectorBasis->getCardinality();
181 int numScalarFields = scalarBasis->getCardinality();
182 int numTotalFields = numVectorFields + numScalarFields;
185 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*(basis_order+1));
186 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.
create(side, 2*(basis_order+1));
188 int numCubPointsCell = cellCub->getNumPoints();
189 int numCubPointsSide = sideCub->getNumPoints();
200 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
232 double zero = (basis_order+1)*(basis_order+1)*100*INTREPID_TOL;
237 cellCub->getCubature(cub_points_cell, cub_weights_cell);
238 sideCub->getCubature(cub_points_side, cub_weights_side);
241 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
244 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
249 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
254 cub_weights_cell.resize(1,numCubPointsCell);
255 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
257 value_of_v_basis_at_cub_points_cell );
258 cub_weights_cell.resize(numCubPointsCell);
261 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
262 FunctionSpaceTools::integrate<double>(fe_matrix_M,
263 w_value_of_v_basis_at_cub_points_cell ,
264 value_of_v_basis_at_cub_points_cell ,
266 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
269 cub_weights_cell.resize(1,numCubPointsCell);
270 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
272 div_of_v_basis_at_cub_points_cell);
273 cub_weights_cell.resize(numCubPointsCell);
275 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
276 FunctionSpaceTools::integrate<double>(fe_matrix_B,
277 w_div_of_v_basis_at_cub_points_cell ,
278 value_of_s_basis_at_cub_points_cell ,
280 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
285 for (
int x_order=0;x_order<=basis_order;x_order++) {
286 for (
int y_order=0;y_order<=basis_order-x_order;y_order++) {
289 fe_matrix.initialize();
291 for (
int i=0;i<numVectorFields;i++) {
292 for (
int j=0;j<numVectorFields;j++) {
293 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
298 for (
int i=0;i<numVectorFields;i++) {
299 for (
int j=0;j<numScalarFields;j++) {
300 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
301 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
306 rhs_vector_vec.initialize();
307 rhs_vector_scal.initialize();
308 rhs_and_soln_vec.initialize();
313 cub_points_cell.resize(1,numCubPointsCell,cellDim);
314 rhsFunc(rhs_at_cub_points_cell,
319 cub_points_cell.resize(numCubPointsCell,cellDim);
321 cub_weights_cell.resize(1,numCubPointsCell);
322 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
324 value_of_s_basis_at_cub_points_cell);
325 cub_weights_cell.resize(numCubPointsCell);
326 FunctionSpaceTools::integrate<double>(rhs_vector_scal,
327 rhs_at_cub_points_cell,
328 w_value_of_s_basis_at_cub_points_cell,
331 for (
int i=0;i<numScalarFields;i++) {
332 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
337 for (
unsigned side_cur=0;side_cur<3;side_cur++) {
346 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
347 u_exact(diri_data_at_cub_points_side,
348 cub_points_side_refcell,x_order,y_order);
349 cub_points_side_refcell.resize(numCubPointsSide,cellDim);
353 (
int)side_cur,cell );
356 vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
357 cub_points_side_refcell ,
361 for (
int i=0;i<numVectorFields;i++) {
362 for (
int j=0;j<numCubPointsSide;j++) {
363 n_of_v_basis_at_cub_points_side(i,j) = 0.0;
364 for (
int k=0;k<cellDim;k++) {
365 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) *
366 value_of_v_basis_at_cub_points_side(i,j,k);
371 cub_weights_side.resize(1,numCubPointsSide);
372 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
374 n_of_v_basis_at_cub_points_side);
375 cub_weights_side.resize(numCubPointsSide);
377 FunctionSpaceTools::integrate<double>(rhs_vector_vec,
378 diri_data_at_cub_points_side,
379 w_n_of_v_basis_at_cub_points_side,
382 for (
int i=0;i<numVectorFields;i++) {
383 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
390 Teuchos::LAPACK<int, double> solver;
391 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0],
392 numTotalFields, &info);
395 scalarBasis->getValues(value_of_s_basis_at_interp_points,
398 for (
int pt=0;pt<numInterpPoints;pt++) {
399 interpolant(0,pt)=0.0;
400 for (
int i=0;i<numScalarFields;i++) {
401 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
402 * value_of_s_basis_at_interp_points(i,pt);
406 interp_points_ref.resize(1,numInterpPoints,cellDim);
409 u_exact( exact_solution , interp_points_ref , x_order, y_order);
410 interp_points_ref.resize(numInterpPoints,cellDim);
416 *outStream <<
"\nNorm-2 error between scalar components of exact solution polynomial of order ("
417 << x_order <<
", " << y_order <<
") and finite element interpolant of order " << basis_order <<
": "
421 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
422 << x_order <<
", " << y_order <<
") and basis order (scalar / vector) ("
423 << basis_order <<
", " << basis_order + 1 <<
")\n\n";
432 catch (
const std::logic_error & err) {
433 *outStream << err.what() <<
"\n\n";
438 std::cout <<
"End Result: TEST FAILED\n";
440 std::cout <<
"End Result: TEST PASSED\n";
443 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 the Intrepid::HDIV_TRI_In_FEM class.
Header file for utility class to provide multidimensional containers.
Implementation of the default H(div)-compatible Raviart-Thomas basis of arbitrary degree on Triangle ...
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.