53 #include "Intrepid_HGRAD_TET_C1_FEM.hpp" 
   56 #include "Teuchos_oblackholestream.hpp" 
   57 #include "Teuchos_RCP.hpp" 
   58 #include "Teuchos_GlobalMPISession.hpp" 
   61 using namespace Intrepid;
 
   63 #define INTREPID_TEST_COMMAND( S )                                                                                  \ 
   68   catch (std::logic_error err) {                                                                                    \ 
   69       *outStream << "Expected Error ----------------------------------------------------------------\n";            \ 
   70       *outStream << err.what() << '\n';                                                                             \ 
   71       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \ 
   76 int main(
int argc, 
char *argv[]) {
 
   78   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
 
   82   int iprint     = argc - 1;
 
   83   Teuchos::RCP<std::ostream> outStream;
 
   84   Teuchos::oblackholestream bhs; 
 
   86     outStream = Teuchos::rcp(&std::cout, 
false);
 
   88     outStream = Teuchos::rcp(&bhs, 
false);
 
   91   Teuchos::oblackholestream oldFormatState;
 
   92   oldFormatState.copyfmt(std::cout);
 
   95   << 
"===============================================================================\n" \
 
   97   << 
"|                      Unit Test (FunctionSpaceTools)                         |\n" \
 
   99   << 
"|     1) basic operator transformations and integration in HGRAD              |\n" \
 
  101   << 
"|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
 
  102   << 
"|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
 
  104   << 
"|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
 
  105   << 
"|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
 
  107   << 
"===============================================================================\n";
 
  111 #ifdef HAVE_INTREPID_DEBUG 
  112   int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
 
  113   int endThrowNumber = beginThrowNumber + 28;
 
  120   << 
"===============================================================================\n"\
 
  121   << 
"| TEST 1: exceptions                                                          |\n"\
 
  122   << 
"===============================================================================\n";
 
  125 #ifdef HAVE_INTREPID_DEBUG 
  141     *outStream << 
"\n >>>>> TESTING computeCellMeasure:\n";
 
  142     INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2, a_2) );
 
  143     INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2_2, a_2) );
 
  145     *outStream << 
"\n >>>>> TESTING computeFaceMeasure:\n";
 
  146     INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
 
  147     INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2_2_3_3, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
 
  149     *outStream << 
"\n >>>>> TESTING computeEdgeMeasure:\n";
 
  150     INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
 
  151     INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2_2_2_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
 
  153     *outStream << 
"\n >>>>> TESTING integrate:\n";
 
  154     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  155     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
 
  156     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
 
  157     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
 
  158     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  159     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  160     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  161     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  162     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  163     INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  165     *outStream << 
"\n >>>>> TESTING operatorIntegral:\n";
 
  166     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  167     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
 
  168     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  169     INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  171     *outStream << 
"\n >>>>> TESTING functionalIntegral:\n";
 
  172     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_2_3_3, a_2_2_2, COMP_CPP) );
 
  173     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
 
  174     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
 
  175     INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
 
  177     *outStream << 
"\n >>>>> TESTING dataIntegral:\n";
 
  178     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2, a_2_2_2, COMP_CPP) );
 
  179     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
 
  180     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
 
  181     INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
 
  183     *outStream << 
"\n >>>>> TESTING applyLeftFieldSigns:\n";
 
  184     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2, a_2) );
 
  185     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2) );
 
  186     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_3_2) );
 
  187     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_3) );
 
  188     INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_2) );
 
  190     *outStream << 
"\n >>>>> TESTING applyRightFieldSigns:\n";
 
  191     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2, a_2) );
 
  192     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2) );
 
  193     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_3_2) );
 
  194     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_3) );
 
  195     INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_2) );
 
  197     *outStream << 
"\n >>>>> TESTING applyFieldSigns:\n";
 
  198     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2, a_2) );
 
  199     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2) );
 
  200     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_3_2) );
 
  201     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2_3) );
 
  202     INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2_2_3_3, a_2_2) );
 
  204     *outStream << 
"\n >>>>> TESTING evaluate:\n";
 
  205     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2) );
 
  206     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2_2_3_3) );
 
  207     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2_2, a_2_2_2_3_3) );
 
  208     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_3_2, a_2_2_2_3_3) );
 
  209     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_2_3, a_2_2_2_3_3) );
 
  210     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_3_2_2_2, a_2_2, a_2_2_2_2_2) );
 
  211     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_3_2_2, a_2_2, a_2_2_2_2_2) );
 
  212     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_2, a_2_2, a_2_2_2_2_2) );
 
  213     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_3, a_2_2, a_2_2_2_2_2) );
 
  214     INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_2, a_2_2, a_2_2_2_2_2) );
 
  217   catch (std::logic_error err) {
 
  218     *outStream << 
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
 
  219     *outStream << err.what() << 
'\n';
 
  220     *outStream << 
"-------------------------------------------------------------------------------" << 
"\n\n";
 
  224 #ifdef HAVE_INTREPID_DEBUG 
  225   if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
 
  231   << 
"===============================================================================\n"\
 
  232   << 
"| TEST 2: correctness of math operations                                      |\n"\
 
  233   << 
"===============================================================================\n";
 
  235   outStream->precision(20);
 
  239       shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >();
 
  247       Teuchos::RCP<Cubature<double, FieldContainer<double, 1>, 
FieldContainer<double, 2> > > myCub = cubFactory.create(cellType, cubDegree);
 
  249       int spaceDim = myCub->getDimension();
 
  251       int numCubPoints = myCub->getNumPoints();
 
  262       int numCellData = numCells*numNodes*spaceDim;
 
  263       double tetnodes[] = {
 
  308       myCub->getCubature(cub_points, cub_weights);
 
  311       cell_nodes.setValues(tetnodes, numCellData);
 
  319       fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
 
  323       tetBasis.
getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
 
  326       fst::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
 
  328                                       grad_of_basis_at_cub_points);
 
  331       fst::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
 
  333                                    transformed_grad_of_basis_at_cub_points);
 
  336       fst::integrate<double>(stiffness_matrices,
 
  337                              transformed_grad_of_basis_at_cub_points,
 
  338                              weighted_transformed_grad_of_basis_at_cub_points,
 
  344       tetBasis.
getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
 
  347       fst::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
 
  348                                        value_of_basis_at_cub_points);
 
  351       fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
 
  353                                    transformed_value_of_basis_at_cub_points);
 
  356       fst::integrate<double>(mass_matrices,
 
  357                              transformed_value_of_basis_at_cub_points,
 
  358                              weighted_transformed_value_of_basis_at_cub_points,
 
  365       string basedir = 
"./testdata";
 
  366       for (
int cell_id = 0; cell_id < numCells; cell_id++) {
 
  368         stringstream namestream;
 
  370         namestream <<  basedir << 
"/mass_TET_FEM_P1" << 
"_" << 
"0" << cell_id+1 << 
".dat";
 
  371         namestream >> filename;
 
  373         ifstream massfile(&filename[0]);
 
  374         if (massfile.is_open()) {
 
  375           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, 0) > 0)
 
  381           std::cout << 
"End Result: TEST FAILED\n";
 
  386         namestream << basedir << 
"/stiff_TET_FEM_P1" << 
"_" << 
"0" << cell_id+1 << 
".dat";
 
  387         namestream >> filename;
 
  389         ifstream stifffile(&filename[0]);
 
  390         if (stifffile.is_open())
 
  392           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, 0) > 0)
 
  398           std::cout << 
"End Result: TEST FAILED\n";
 
  408   catch (std::logic_error err) {
 
  409     *outStream << 
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
 
  410     *outStream << err.what() << 
'\n';
 
  411     *outStream << 
"-------------------------------------------------------------------------------" << 
"\n\n";
 
  417     std::cout << 
"End Result: TEST FAILED\n";
 
  419     std::cout << 
"End Result: TEST PASSED\n";
 
  422   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. 
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
A factory class that generates specific instances of cubatures. 
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell...