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;
68 const shards::CellTopology & ,
69 int ,
int ,
int ,
int ,
int );
79 for (
int cell=0;cell<result.
dimension(0);cell++){
80 for (
int pt=0;pt<result.
dimension(1);pt++) {
81 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd)
82 *std::pow(points(cell,pt,2),zd);
90 const shards::CellTopology & parentCell ,
91 int sideOrdinal ,
int comp ,
int a,
int b,
int c )
104 for (
int cell=0;cell<numCells;cell++) {
105 for (
int pt=0;pt<numPoints;pt++) {
108 result(cell,pt,0) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
111 result(cell,pt,0) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
115 result(cell,pt,1) = b * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
119 result(cell,pt,2) = c * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
124 else if (comp == 1) {
125 for (
int cell=0;cell<numCells;cell++) {
126 for (
int pt=0;pt<numPoints;pt++) {
129 result(cell,pt,0) = a * normal(1) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
133 result(cell,pt,1) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
136 result(cell,pt,1) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
140 result(cell,pt,2) = c * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
145 else if (comp == 2) {
146 for (
int cell=0;cell<numCells;cell++) {
147 for (
int pt=0;pt<numPoints;pt++) {
150 result(cell,pt,0) = a * normal(2) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
154 result(cell,pt,1) = b * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
158 result(cell,pt,2) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1)* std::pow(points(cell,pt,2),c);
161 result(cell,pt,2) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
179 u_exact(result,points,comp,a,b,c);
184 for (
int cell=0;cell<result.
dimension(0);cell++) {
185 for (
int pt=0;pt<result.
dimension(1);pt++) {
186 result(cell,pt,0) -= b*(b-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b-2)*std::pow(points(cell,pt,2),c);
191 for (
int cell=0;cell<result.
dimension(0);cell++) {
192 for (
int pt=0;pt<result.
dimension(1);pt++) {
193 result(cell,pt,0)-=c*(c-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
198 if ( (a>=1) && (b>=1) ) {
199 for (
int cell=0;cell<result.
dimension(0);cell++) {
200 for (
int pt=0;pt<result.
dimension(1);pt++) {
201 result(cell,pt,1)+=a*b*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b-1)*std::pow(points(cell,pt,2),c);
206 if ( (a>=1) && (c>=1) ) {
207 for (
int cell=0;cell<result.
dimension(0);cell++) {
208 for (
int pt=0;pt<result.
dimension(1);pt++) {
209 result(cell,pt,2)+=a*c*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-1);
215 for (
int cell=0;cell<result.
dimension(0);cell++) {
216 for (
int pt=0;pt<result.
dimension(1);pt++) {
218 if (a > 0 && b > 0) {
219 result(cell,pt,0) += a * b * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c);
223 result(cell,pt,1) -= c*(c-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2);
226 result(cell,pt,1) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
229 if (b > 0 && c > 0) {
230 result(cell,pt,2) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
235 else if (comp == 2) {
236 for (
int cell=0;cell<result.
dimension(0);cell++) {
237 for (
int pt=0;pt<result.
dimension(1);pt++) {
240 result(cell,pt,0) += a * c * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1);
244 result(cell,pt,1) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1);
248 result(cell,pt,2) -= b*(b-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-2) * std::pow(points(cell,pt,2),c);
251 result(cell,pt,2) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c);
259 int main(
int argc,
char *argv[]) {
260 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
264 int iprint = argc - 1;
265 Teuchos::RCP<std::ostream> outStream;
266 Teuchos::oblackholestream bhs;
268 outStream = Teuchos::rcp(&std::cout,
false);
270 outStream = Teuchos::rcp(&bhs,
false);
273 Teuchos::oblackholestream oldFormatState;
274 oldFormatState.copyfmt(std::cout);
277 <<
"===============================================================================\n" \
279 <<
"| Unit Test (Basis_HGRAD_TRI_In_FEM) |\n" \
281 <<
"| 1) Patch test involving H(curl) matrices |\n" \
283 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
284 <<
"| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \
285 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
286 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
288 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
289 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
291 <<
"===============================================================================\n" \
292 <<
"| TEST 1: Patch test for mass + curl-curl matrices |\n" \
293 <<
"===============================================================================\n";
298 outStream -> precision(16);
302 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());
304 shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());
306 int cellDim = cell.getDimension();
307 int sideDim = side.getDimension();
312 int numIntervals = max_order;
313 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6;
316 for (
int j=0; j<=numIntervals; j++) {
317 for (
int i=0; i<=numIntervals-j; i++) {
318 for (
int k=0;k<numIntervals-j-i;k++) {
319 interp_points_ref(counter,0) = i*(1.0/numIntervals);
320 interp_points_ref(counter,1) = j*(1.0/numIntervals);
321 interp_points_ref(counter,2) = k*(1.0/numIntervals);
327 for (
int basis_order=min_order;basis_order<=max_order;basis_order++) {
329 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis =
332 int numFields = basis->getCardinality();
335 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.
create(cell, 2*(basis_order+1));
336 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.
create(side, 2*(basis_order+1));
338 int numCubPointsCell = cellCub->getNumPoints();
339 int numCubPointsSide = sideCub->getNumPoints();
354 FieldContainer<double> w_value_of_basis_at_cub_points_side_refcell(1,numFields,numCubPointsSide,cellDim );
371 Teuchos::LAPACK<int, double> solver;
374 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL;
378 cellCub->getCubature(cub_points_cell, cub_weights_cell);
379 sideCub->getCubature(cub_points_side, cub_weights_side);
382 basis->getValues(value_of_basis_at_cub_points_cell,
385 basis->getValues(curl_of_basis_at_cub_points_cell,
389 basis->getValues( value_of_basis_at_interp_points ,
395 cub_weights_cell.resize(1,numCubPointsCell);
396 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell ,
398 value_of_basis_at_cub_points_cell );
399 FunctionSpaceTools::multiplyMeasure<double>(w_curl_of_basis_at_cub_points_cell ,
401 curl_of_basis_at_cub_points_cell );
402 cub_weights_cell.resize(numCubPointsCell);
405 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
406 curl_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim );
407 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
408 w_value_of_basis_at_cub_points_cell ,
409 value_of_basis_at_cub_points_cell ,
411 FunctionSpaceTools::integrate<double>(fe_matrix_bak,
412 w_curl_of_basis_at_cub_points_cell ,
413 curl_of_basis_at_cub_points_cell ,
416 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
417 curl_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim );
419 for (
int comp=0;comp<cellDim;comp++) {
420 for (
int x_order=0;x_order<basis_order;x_order++) {
421 for (
int y_order=0;y_order<basis_order-x_order;y_order++) {
422 for (
int z_order=0;z_order<basis_order-x_order-y_order;z_order++) {
423 fe_matrix.initialize();
425 for (
int i=0;i<numFields;i++) {
426 for (
int j=0;j<numFields;j++) {
427 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j);
432 rhs_and_soln_vec.initialize();
435 cub_points_cell.resize(1,numCubPointsCell,cellDim);
437 rhs_at_cub_points_cell.initialize();
438 rhsFunc(rhs_at_cub_points_cell,
445 cub_points_cell.resize(numCubPointsCell,cellDim);
447 cub_weights_cell.resize(numCubPointsCell);
449 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec,
450 rhs_at_cub_points_cell,
451 w_value_of_basis_at_cub_points_cell,
456 for (
unsigned i=0;i<4;i++) {
462 basis->getValues( value_of_basis_at_cub_points_side_refcell ,
463 cub_points_side_refcell ,
467 cub_points_side_refcell.resize( 1 , numCubPointsSide , cellDim );
468 bcFunc(bc_func_at_cub_points_side_refcell,
469 cub_points_side_refcell,
476 cub_points_side_refcell.resize( numCubPointsSide , cellDim );
481 cub_weights_side.resize(1,numCubPointsSide);
482 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_side_refcell ,
484 value_of_basis_at_cub_points_side_refcell );
485 cub_weights_side.resize(numCubPointsSide);
487 FunctionSpaceTools::integrate<double>( bc_fields_per_side ,
488 bc_func_at_cub_points_side_refcell ,
489 w_value_of_basis_at_cub_points_side_refcell ,
498 solver.POTRF(
'L',numFields,&fe_matrix[0],numFields,&info);
499 solver.POTRS(
'L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info);
501 interp_points_ref.resize(1,numInterpPoints,cellDim);
504 exact_solution.initialize();
505 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order);
506 interp_points_ref.resize(numInterpPoints,cellDim);
510 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim);
511 FunctionSpaceTools::evaluate<double>( interpolant ,
513 value_of_basis_at_interp_points );
514 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim);
520 *outStream <<
"\nNorm-2 error between scalar components of exact solution of order ("
521 << x_order <<
", " << y_order <<
", " << z_order
522 <<
") in component " << comp
523 <<
" and finite element interpolant of order " << basis_order <<
": "
527 *outStream <<
"\n\nPatch test failed for solution polynomial order ("
528 << x_order <<
", " << y_order <<
", " << z_order <<
") and basis order "
529 << basis_order <<
"\n\n";
540 catch (
const std::logic_error & err) {
541 *outStream << err.what() <<
"\n\n";
546 std::cout <<
"End Result: TEST FAILED\n";
548 std::cout <<
"End Result: TEST PASSED\n";
551 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::HCURL_TET_In_FEM class.
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.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Te...
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
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.