50#include "Intrepid_HGRAD_TET_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"
65using namespace Intrepid;
67void rhsFunc( FieldContainer<double> &,
const FieldContainer<double> &,
int,
int,
int );
68void u_exact( FieldContainer<double> &,
const FieldContainer<double> &,
int,
int,
int );
72void rhsFunc( FieldContainer<double> &result,
73 const FieldContainer<double> &points,
78 for (
int cell=0;cell<result.dimension(0);cell++) {
79 for (
int pt=0;pt<result.dimension(1);pt++) {
80 result(cell,pt) = 0.0;
82 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd)
83 *pow(points(cell,pt,2),zd);
86 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2)
87 *pow(points(cell,pt,2),zd);
90 result(cell,pt) += zd*(zd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd)
91 *pow(points(cell,pt,2),zd-2);
98void u_exact( FieldContainer<double> &result,
99 const FieldContainer<double> &points,
104 for (
int cell=0;cell<result.dimension(0);cell++){
105 for (
int pt=0;pt<result.dimension(1);pt++) {
106 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
107 *std::pow(points(cell,pt,2),zd);
113int main(
int argc,
char *argv[]) {
114 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
118 int iprint = argc - 1;
119 Teuchos::RCP<std::ostream> outStream;
120 Teuchos::oblackholestream bhs;
122 outStream = Teuchos::rcp(&std::cout,
false);
124 outStream = Teuchos::rcp(&bhs,
false);
127 Teuchos::oblackholestream oldFormatState;
128 oldFormatState.copyfmt(std::cout);
131 <<
"===============================================================================\n" \
133 <<
"| Unit Test (Basis_HGRAD_TET_In_FEM) |\n" \
135 <<
"| 1) Patch test involving H(div) matrices |\n" \
136 <<
"| for the Dirichlet problem on a tetrahedral patch |\n" \
137 <<
"| Omega with boundary Gamma. |\n" \
139 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
140 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
141 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
142 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
144 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
145 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
147 <<
"===============================================================================\n" \
148 <<
"| TEST 1: Patch test |\n" \
149 <<
"===============================================================================\n";
154 outStream -> precision(16);
157 DefaultCubatureFactory<double> cubFactory;
158 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());
159 shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());
161 int cellDim = cell.getDimension();
162 int sideDim = side.getDimension();
167 int numIntervals = max_order;
168 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
169 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim);
171 for (
int j=0; j<=numIntervals; j++) {
172 for (
int i=0; i<=numIntervals-j; i++) {
173 for (
int k=0;k<numIntervals-j-i;k++) {
174 interp_points_ref(counter,0) = i*(1.0/numIntervals);
175 interp_points_ref(counter,1) = j*(1.0/numIntervals);
176 interp_points_ref(counter,2) = k*(1.0/numIntervals);
184 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
186 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis =
187 Teuchos::rcp(
new Basis_HDIV_TET_In_FEM<
double,FieldContainer<double> >(basis_order+1,POINTTYPE_EQUISPACED) );
188 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis =
191 int numVectorFields = vectorBasis->getCardinality();
192 int numScalarFields = scalarBasis->getCardinality();
193 int numTotalFields = numVectorFields + numScalarFields;
196 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1));
197 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1));
199 int numCubPointsCell = cellCub->getNumPoints();
200 int numCubPointsSide = sideCub->getNumPoints();
203 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim);
204 FieldContainer<double> cub_weights_cell(numCubPointsCell);
205 FieldContainer<double> cub_points_side( numCubPointsSide, sideDim );
206 FieldContainer<double> cub_weights_side( numCubPointsSide );
207 FieldContainer<double> cub_points_side_refcell( numCubPointsSide , cellDim );
210 FieldContainer<double> value_of_v_basis_at_cub_points_cell(numVectorFields, numCubPointsCell, cellDim );
211 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim);
212 FieldContainer<double> div_of_v_basis_at_cub_points_cell( numVectorFields, numCubPointsCell );
213 FieldContainer<double> w_div_of_v_basis_at_cub_points_cell( 1, numVectorFields , numCubPointsCell );
214 FieldContainer<double> value_of_s_basis_at_cub_points_cell(numScalarFields,numCubPointsCell);
215 FieldContainer<double> w_value_of_s_basis_at_cub_points_cell(1,numScalarFields,numCubPointsCell);
220 FieldContainer<double> value_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide,cellDim);
221 FieldContainer<double> n_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide);
222 FieldContainer<double> w_n_of_v_basis_at_cub_points_side(1,numVectorFields,numCubPointsSide);
223 FieldContainer<double> diri_data_at_cub_points_side(1,numCubPointsSide);
224 FieldContainer<double> side_normal(cellDim);
227 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell);
230 FieldContainer<double> fe_matrix_M(1,numVectorFields,numVectorFields);
231 FieldContainer<double> fe_matrix_B(1,numVectorFields,numScalarFields);
232 FieldContainer<double> fe_matrix(1,numTotalFields,numTotalFields);
234 FieldContainer<double> rhs_vector_vec(1,numVectorFields);
235 FieldContainer<double> rhs_vector_scal(1,numScalarFields);
236 FieldContainer<double> rhs_and_soln_vec(1,numTotalFields);
238 FieldContainer<int> ipiv(numTotalFields);
239 FieldContainer<double> value_of_s_basis_at_interp_points( numScalarFields , numInterpPoints);
240 FieldContainer<double> interpolant( 1 , numInterpPoints );
243 double zero = (basis_order+1)*(basis_order+1)*100*INTREPID_TOL;
248 cellCub->getCubature(cub_points_cell, cub_weights_cell);
249 sideCub->getCubature(cub_points_side, cub_weights_side);
252 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell,
255 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell,
260 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell,
265 cub_weights_cell.resize(1,numCubPointsCell);
266 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell ,
268 value_of_v_basis_at_cub_points_cell );
269 cub_weights_cell.resize(numCubPointsCell);
272 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim );
273 FunctionSpaceTools::integrate<double>(fe_matrix_M,
274 w_value_of_v_basis_at_cub_points_cell ,
275 value_of_v_basis_at_cub_points_cell ,
277 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim );
280 cub_weights_cell.resize(1,numCubPointsCell);
281 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell,
283 div_of_v_basis_at_cub_points_cell);
284 cub_weights_cell.resize(numCubPointsCell);
286 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell);
287 FunctionSpaceTools::integrate<double>(fe_matrix_B,
288 w_div_of_v_basis_at_cub_points_cell ,
289 value_of_s_basis_at_cub_points_cell ,
291 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell);
296 for (
int x_order=0;x_order<=basis_order;x_order++) {
297 for (
int y_order=0;y_order<=basis_order-x_order;y_order++) {
298 for (
int z_order=0;z_order<=basis_order-x_order-y_order;z_order++) {
302 fe_matrix.initialize();
304 for (
int i=0;i<numVectorFields;i++) {
305 for (
int j=0;j<numVectorFields;j++) {
306 fe_matrix(0,i,j) = fe_matrix_M(0,i,j);
311 for (
int i=0;i<numVectorFields;i++) {
312 for (
int j=0;j<numScalarFields;j++) {
313 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j);
314 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j);
319 rhs_vector_vec.initialize();
320 rhs_vector_scal.initialize();
321 rhs_and_soln_vec.initialize();
326 cub_points_cell.resize(1,numCubPointsCell,cellDim);
327 rhsFunc(rhs_at_cub_points_cell,
333 cub_points_cell.resize(numCubPointsCell,cellDim);
335 cub_weights_cell.resize(1,numCubPointsCell);
336 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell,
338 value_of_s_basis_at_cub_points_cell);
339 cub_weights_cell.resize(numCubPointsCell);
340 FunctionSpaceTools::integrate<double>(rhs_vector_scal,
341 rhs_at_cub_points_cell,
342 w_value_of_s_basis_at_cub_points_cell,
345 for (
int i=0;i<numScalarFields;i++) {
346 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i);
351 for (
unsigned side_cur=0;side_cur<4;side_cur++) {
360 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim);
361 u_exact(diri_data_at_cub_points_side,
362 cub_points_side_refcell,x_order,y_order,z_order);
363 cub_points_side_refcell.resize(numCubPointsSide,cellDim);
367 (
int)side_cur,cell );
370 vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
371 cub_points_side_refcell ,
375 for (
int i=0;i<numVectorFields;i++) {
376 for (
int j=0;j<numCubPointsSide;j++) {
377 n_of_v_basis_at_cub_points_side(i,j) = 0.0;
378 for (
int k=0;k<cellDim;k++) {
379 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) *
380 value_of_v_basis_at_cub_points_side(i,j,k);
385 cub_weights_side.resize(1,numCubPointsSide);
386 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
388 n_of_v_basis_at_cub_points_side);
389 cub_weights_side.resize(numCubPointsSide);
391 FunctionSpaceTools::integrate<double>(rhs_vector_vec,
392 diri_data_at_cub_points_side,
393 w_n_of_v_basis_at_cub_points_side,
396 for (
int i=0;i<numVectorFields;i++) {
397 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
404 Teuchos::LAPACK<int, double> solver;
405 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0],
406 numTotalFields, &info);
409 scalarBasis->getValues(value_of_s_basis_at_interp_points,
412 for (
int pt=0;pt<numInterpPoints;pt++) {
413 interpolant(0,pt)=0.0;
414 for (
int i=0;i<numScalarFields;i++) {
415 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
416 * value_of_s_basis_at_interp_points(i,pt);
420 interp_points_ref.resize(1,numInterpPoints,cellDim);
422 FieldContainer<double> exact_solution(1,numInterpPoints);
423 u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
424 interp_points_ref.resize(numInterpPoints,cellDim);
430 *outStream <<
"\nNorm-2 error between scalar components of exact solution of order ("
431 << x_order <<
", " << y_order <<
", " << z_order
432 <<
") and finite element interpolant of order " << basis_order <<
": "
436 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
437 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order (scalar, vector) ("
438 << basis_order <<
", " << basis_order+1 <<
")\n\n";
449 catch (
const std::logic_error & err) {
450 *outStream << err.what() <<
"\n\n";
455 std::cout <<
"End Result: TEST FAILED\n";
457 std::cout <<
"End Result: TEST PASSED\n";
460 std::cout.copyfmt(oldFormatState);
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_TET_In_FEM class.
Implementation of the default H(div)-compatible Raviart-Thomas basis of arbitrary degree on Tetrahedr...
Implementation of the default H(grad)-compatible orthogonal basis of arbitrary degree on tetrahedron.