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...