52 #include "Intrepid_HGRAD_TET_C1_FEM.hpp" 
   55 #include "Teuchos_oblackholestream.hpp" 
   56 #include "Teuchos_RCP.hpp" 
   57 #include "Teuchos_GlobalMPISession.hpp" 
   60 using namespace Intrepid;
 
   62 #define INTREPID_TEST_COMMAND( S )                                                                                  \ 
   67   catch (const std::logic_error & err) {                                                                                    \ 
   68       *outStream << "Expected Error ----------------------------------------------------------------\n";            \ 
   69       *outStream << err.what() << '\n';                                                                             \ 
   70       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \ 
   75 int main(
int argc, 
char *argv[]) {
 
   77   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
 
   81   int iprint     = argc - 1;
 
   82   Teuchos::RCP<std::ostream> outStream;
 
   83   Teuchos::oblackholestream bhs; 
 
   85     outStream = Teuchos::rcp(&std::cout, 
false);
 
   87     outStream = Teuchos::rcp(&bhs, 
false);
 
   90   Teuchos::oblackholestream oldFormatState;
 
   91   oldFormatState.copyfmt(std::cout);
 
   94   << 
"===============================================================================\n" \
 
   96   << 
"|                      Unit Test (FunctionSpaceTools)                         |\n" \
 
   98   << 
"|     1) basic operator transformations and integration in HGRAD              |\n" \
 
  100   << 
"|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
 
  101   << 
"|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
 
  103   << 
"|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
 
  104   << 
"|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
 
  106   << 
"===============================================================================\n";
 
  110 #ifdef HAVE_INTREPID_DEBUG 
  111   int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
 
  112   int endThrowNumber = beginThrowNumber + 28;
 
  119   << 
"===============================================================================\n"\
 
  120   << 
"| TEST 1: exceptions                                                          |\n"\
 
  121   << 
"===============================================================================\n";
 
  124 #ifdef HAVE_INTREPID_DEBUG 
  140     *outStream << 
"\n >>>>> TESTING computeCellMeasure:\n";
 
  141     INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2, a_2) );
 
  142     INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2_2, a_2) );
 
  144     *outStream << 
"\n >>>>> TESTING computeFaceMeasure:\n";
 
  145     INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
 
  146     INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2_2_3_3, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
 
  148     *outStream << 
"\n >>>>> TESTING computeEdgeMeasure:\n";
 
  149     INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
 
  150     INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2_2_2_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
 
  152     *outStream << 
"\n >>>>> TESTING integrate:\n";
 
  153     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  154     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
 
  155     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
 
  156     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
 
  157     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  158     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  159     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  160     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  161     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  162     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  164     *outStream << 
"\n >>>>> TESTING operatorIntegral:\n";
 
  165     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  166     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  167     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  168     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  170     *outStream << 
"\n >>>>> TESTING functionalIntegral:\n";
 
  171     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_2_3_3, a_2_2_2, COMP_CPP) );
 
  172     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  173     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  174     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  176     *outStream << 
"\n >>>>> TESTING dataIntegral:\n";
 
  177     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2, a_2_2_2, COMP_CPP) );
 
  178     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
 
  179     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
 
  180     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
 
  182     *outStream << 
"\n >>>>> TESTING applyLeftFieldSigns:\n";
 
  183     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2, a_2) );
 
  184     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2) );
 
  185     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_3_2) );
 
  186     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_3) );
 
  187     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_2) );
 
  189     *outStream << 
"\n >>>>> TESTING applyRightFieldSigns:\n";
 
  190     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2, a_2) );
 
  191     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2) );
 
  192     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_3_2) );
 
  193     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_3) );
 
  194     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_2) );
 
  196     *outStream << 
"\n >>>>> TESTING applyFieldSigns:\n";
 
  197     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2, a_2) );
 
  198     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2) );
 
  199     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_3_2) );
 
  200     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2_3) );
 
  201     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2_2_3_3, a_2_2) );
 
  203     *outStream << 
"\n >>>>> TESTING evaluate:\n";
 
  204     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2) );
 
  205     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2_2_3_3) );
 
  206     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2_2, a_2_2_2_3_3) );
 
  207     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_3_2, a_2_2_2_3_3) );
 
  208     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_2_3, a_2_2_2_3_3) );
 
  209     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_3_2_2_2, a_2_2, a_2_2_2_2_2) );
 
  210     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_3_2_2, a_2_2, a_2_2_2_2_2) );
 
  211     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_2, a_2_2, a_2_2_2_2_2) );
 
  212     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_3, a_2_2, a_2_2_2_2_2) );
 
  213     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_2, a_2_2, a_2_2_2_2_2) );
 
  216   catch (
const std::logic_error & err) {
 
  217     *outStream << 
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
 
  218     *outStream << err.what() << 
'\n';
 
  219     *outStream << 
"-------------------------------------------------------------------------------" << 
"\n\n";
 
  223 #ifdef HAVE_INTREPID_DEBUG 
  224   if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
 
  230   << 
"===============================================================================\n"\
 
  231   << 
"| TEST 2: correctness of math operations                                      |\n"\
 
  232   << 
"===============================================================================\n";
 
  234   outStream->precision(20);
 
  237       shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >();   
 
  242       Teuchos::RCP<Cubature<double> > myCub = cubFactory.
create(cellType, cubDegree);           
 
  243       int spaceDim = myCub->getDimension();                                                     
 
  244       int numCubPoints = myCub->getNumPoints();                                                 
 
  253       int numCellData = numCells*numNodes*spaceDim;
 
  254       double tetnodes[] = {
 
  287       FieldContainer<double> transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
 
  288       FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
 
  293       FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints);
 
  299       myCub->getCubature(cub_points, cub_weights);
 
  302       cell_nodes.setValues(tetnodes, numCellData);
 
  310       fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
 
  314       tetBasis.
getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
 
  317       fst::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
 
  319                                       grad_of_basis_at_cub_points);
 
  322       fst::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
 
  324                                    transformed_grad_of_basis_at_cub_points);
 
  327       fst::integrate<double>(stiffness_matrices,
 
  328                              transformed_grad_of_basis_at_cub_points,
 
  329                              weighted_transformed_grad_of_basis_at_cub_points,
 
  335       tetBasis.
getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
 
  338       fst::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
 
  339                                        value_of_basis_at_cub_points);
 
  342       fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
 
  344                                    transformed_value_of_basis_at_cub_points);
 
  347       fst::integrate<double>(mass_matrices,
 
  348                              transformed_value_of_basis_at_cub_points,
 
  349                              weighted_transformed_value_of_basis_at_cub_points,
 
  356       string basedir = 
"./testdata";
 
  357       for (
int cell_id = 0; cell_id < numCells; cell_id++) {
 
  359         stringstream namestream;
 
  361         namestream <<  basedir << 
"/mass_TET_FEM_P1" << 
"_" << 
"0" << cell_id+1 << 
".dat";
 
  362         namestream >> filename;
 
  364         ifstream massfile(&filename[0]);
 
  365         if (massfile.is_open()) {
 
  366           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, 0) > 0)
 
  372           std::cout << 
"End Result: TEST FAILED\n";
 
  377         namestream << basedir << 
"/stiff_TET_FEM_P1" << 
"_" << 
"0" << cell_id+1 << 
".dat";
 
  378         namestream >> filename;
 
  380         ifstream stifffile(&filename[0]);
 
  381         if (stifffile.is_open())
 
  383           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, 0) > 0)
 
  389           std::cout << 
"End Result: TEST FAILED\n";
 
  399   catch (
const std::logic_error & err) {
 
  400     *outStream << 
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
 
  401     *outStream << err.what() << 
'\n';
 
  402     *outStream << 
"-------------------------------------------------------------------------------" << 
"\n\n";
 
  408     std::cout << 
"End Result: TEST FAILED\n";
 
  410     std::cout << 
"End Result: TEST PASSED\n";
 
  413   std::cout.copyfmt(oldFormatState);
 
virtual int getCardinality() const 
Returns cardinality of the basis. 
Header file for utility class to provide multidimensional containers. 
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const 
FEM basis evaluation on a reference Tetrahedron cell. 
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. 
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell...