55 #include "Teuchos_oblackholestream.hpp"
56 #include "Teuchos_RCP.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
60 using namespace Intrepid;
62 int main(
int argc,
char *argv[]) {
64 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
68 int iprint = argc - 1;
69 Teuchos::RCP<std::ostream> outStream;
70 Teuchos::oblackholestream bhs;
72 outStream = Teuchos::rcp(&std::cout,
false);
74 outStream = Teuchos::rcp(&bhs,
false);
77 Teuchos::oblackholestream oldFormatState;
78 oldFormatState.copyfmt(std::cout);
81 <<
"===============================================================================\n" \
83 <<
"| Unit Test (FunctionSpaceTools) |\n" \
85 <<
"| 1) Basic operator transformations and integration in HDIV: |\n" \
87 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
88 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
90 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
91 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
93 <<
"===============================================================================\n";
102 <<
"===============================================================================\n"\
103 <<
"| TEST 1: correctness of math operations |\n"\
104 <<
"===============================================================================\n";
106 outStream->precision(20);
110 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >();
115 Teuchos::RCP<Cubature<double> > myCub = cubFactory.
create(cellType, cubDegree);
116 int spaceDim = myCub->getDimension();
117 int numCubPoints = myCub->getNumPoints();
126 int numCellData = numCells*numNodes*spaceDim;
127 int numSignData = numCells*numFields;
129 double hexnodes[] = {
168 short facesigns[] = {
187 FieldContainer<double> weighted_transformed_div_of_basis_at_cub_points(numCells, numFields, numCubPoints);
191 FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
192 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
198 myCub->getCubature(cub_points, cub_weights);
201 cell_nodes.setValues(hexnodes, numCellData);
204 field_signs.setValues(facesigns, numSignData);
212 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
216 hexBasis.
getValues(div_of_basis_at_cub_points, cub_points, OPERATOR_DIV);
219 fst::HDIVtransformDIV<double>(transformed_div_of_basis_at_cub_points,
221 div_of_basis_at_cub_points);
224 fst::multiplyMeasure<double>(weighted_transformed_div_of_basis_at_cub_points,
226 transformed_div_of_basis_at_cub_points);
229 fst::applyFieldSigns<double>(transformed_div_of_basis_at_cub_points, field_signs);
230 fst::applyFieldSigns<double>(weighted_transformed_div_of_basis_at_cub_points, field_signs);
233 fst::integrate<double>(stiffness_matrices,
234 transformed_div_of_basis_at_cub_points,
235 weighted_transformed_div_of_basis_at_cub_points,
240 hexBasis.
getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
243 fst::HDIVtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
246 value_of_basis_at_cub_points);
249 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
251 transformed_value_of_basis_at_cub_points);
254 fst::integrate<double>(mass_matrices,
255 transformed_value_of_basis_at_cub_points,
256 weighted_transformed_value_of_basis_at_cub_points,
260 fst::applyLeftFieldSigns<double>(mass_matrices, field_signs);
261 fst::applyRightFieldSigns<double>(mass_matrices, field_signs);
267 string basedir =
"./testdata";
268 for (
int cell_id = 0; cell_id < numCells-1; cell_id++) {
270 stringstream namestream;
272 namestream << basedir <<
"/mass_HDIV_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
273 namestream >> filename;
275 ifstream massfile(&filename[0]);
276 if (massfile.is_open()) {
277 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
283 std::cout <<
"End Result: TEST FAILED\n";
288 namestream << basedir <<
"/stiff_HDIV_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
289 namestream >> filename;
291 ifstream stifffile(&filename[0]);
292 if (stifffile.is_open())
294 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
300 std::cout <<
"End Result: TEST FAILED\n";
305 for (
int cell_id = 3; cell_id < numCells; cell_id++) {
307 stringstream namestream;
309 namestream << basedir <<
"/mass_fp_HDIV_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
310 namestream >> filename;
312 ifstream massfile(&filename[0]);
313 if (massfile.is_open()) {
314 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
320 std::cout <<
"End Result: TEST FAILED\n";
325 namestream << basedir <<
"/stiff_fp_HDIV_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
326 namestream >> filename;
328 ifstream stifffile(&filename[0]);
329 if (stifffile.is_open())
331 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
337 std::cout <<
"End Result: TEST FAILED\n";
346 catch (std::logic_error err) {
347 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
348 *outStream << err.what() <<
'\n';
349 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
355 std::cout <<
"End Result: TEST FAILED\n";
357 std::cout <<
"End Result: TEST PASSED\n";
360 std::cout.copyfmt(oldFormatState);
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual int getCardinality() const
Returns cardinality of the basis.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_HEX_I1_FEM class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedron cell...
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.