55#include "Teuchos_oblackholestream.hpp"
56#include "Teuchos_RCP.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
60using namespace Intrepid;
62int 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<> >();
113 DefaultCubatureFactory<double> cubFactory;
115 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);
116 int spaceDim = myCub->getDimension();
117 int numCubPoints = myCub->getNumPoints();
120 Basis_HDIV_HEX_I1_FEM<double, FieldContainer<double> > hexBasis;
121 int numFields = hexBasis.getCardinality();
126 int numCellData = numCells*numNodes*spaceDim;
127 int numSignData = numCells*numFields;
129 double hexnodes[] = {
168 short facesigns[] = {
176 FieldContainer<double> cub_points(numCubPoints, spaceDim);
177 FieldContainer<double> cub_weights(numCubPoints);
178 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
179 FieldContainer<short> field_signs(numCells, numFields);
180 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
182 FieldContainer<double> jacobian_det(numCells, numCubPoints);
183 FieldContainer<double> weighted_measure(numCells, numCubPoints);
185 FieldContainer<double> div_of_basis_at_cub_points(numFields, numCubPoints);
186 FieldContainer<double> transformed_div_of_basis_at_cub_points(numCells, numFields, numCubPoints);
187 FieldContainer<double> weighted_transformed_div_of_basis_at_cub_points(numCells, numFields, numCubPoints);
188 FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
190 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
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);
193 FieldContainer<double> mass_matrices(numCells, numFields, numFields);
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 (
const 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);
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HDIV_HEX_I1_FEM class.