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);
67 int iprint = argc - 1;
68 Teuchos::RCP<std::ostream> outStream;
69 Teuchos::oblackholestream bhs;
71 outStream = Teuchos::rcp(&std::cout,
false);
73 outStream = Teuchos::rcp(&bhs,
false);
76 Teuchos::oblackholestream oldFormatState;
77 oldFormatState.copyfmt(std::cout);
80 <<
"===============================================================================\n" \
82 <<
"| Unit Test (FunctionSpaceTools) |\n" \
84 <<
"| 1) basic operator transformations and integration in HCURL |\n" \
86 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
87 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
89 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
90 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
92 <<
"===============================================================================\n";
101 <<
"===============================================================================\n"\
102 <<
"| TEST 1: correctness of math operations |\n"\
103 <<
"===============================================================================\n";
105 outStream->precision(20);
108 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >();
111 DefaultCubatureFactory<double> cubFactory;
113 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);
114 int spaceDim = myCub->getDimension();
115 int numCubPoints = myCub->getNumPoints();
117 Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexBasis;
118 int numFields = hexBasis.getCardinality();
123 int numCellData = numCells*numNodes*spaceDim;
124 int numSignData = numCells*numFields;
126 double hexnodes[] = {
165 double edgesigns[] = {
166 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
167 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1,
168 -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1,
169 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1
173 FieldContainer<double> cub_points(numCubPoints, spaceDim);
174 FieldContainer<double> cub_weights(numCubPoints);
175 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
176 FieldContainer<double> field_signs(numCells, numFields);
177 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
178 FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
179 FieldContainer<double> jacobian_det(numCells, numCubPoints);
180 FieldContainer<double> weighted_measure(numCells, numCubPoints);
182 FieldContainer<double> curl_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
183 FieldContainer<double> transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
184 FieldContainer<double> weighted_transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
185 FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
187 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
188 FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
189 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
190 FieldContainer<double> mass_matrices(numCells, numFields, numFields);
195 myCub->getCubature(cub_points, cub_weights);
198 cell_nodes.setValues(hexnodes, numCellData);
200 field_signs.setValues(edgesigns, numSignData);
207 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
211 hexBasis.getValues(curl_of_basis_at_cub_points, cub_points, OPERATOR_CURL);
214 fst::HCURLtransformCURL<double>(transformed_curl_of_basis_at_cub_points,
217 curl_of_basis_at_cub_points);
219 fst::multiplyMeasure<double>(weighted_transformed_curl_of_basis_at_cub_points,
221 transformed_curl_of_basis_at_cub_points);
224 fst::applyFieldSigns<double>(transformed_curl_of_basis_at_cub_points, field_signs);
225 fst::applyFieldSigns<double>(weighted_transformed_curl_of_basis_at_cub_points, field_signs);
228 fst::integrate<double>(stiffness_matrices,
229 transformed_curl_of_basis_at_cub_points,
230 weighted_transformed_curl_of_basis_at_cub_points,
234 hexBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
237 fst::HCURLtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
239 value_of_basis_at_cub_points);
242 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
244 transformed_value_of_basis_at_cub_points);
247 fst::integrate<double>(mass_matrices,
248 transformed_value_of_basis_at_cub_points,
249 weighted_transformed_value_of_basis_at_cub_points,
253 fst::applyLeftFieldSigns<double>(mass_matrices, field_signs);
254 fst::applyRightFieldSigns<double>(mass_matrices, field_signs);
259 string basedir =
"./testdata";
260 for (
int cell_id = 0; cell_id < numCells-1; cell_id++) {
262 stringstream namestream;
264 namestream << basedir <<
"/mass_HCURL_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
265 namestream >> filename;
267 ifstream massfile(&filename[0]);
268 if (massfile.is_open()) {
269 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
275 std::cout <<
"End Result: TEST FAILED\n";
280 namestream << basedir <<
"/stiff_HCURL_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
281 namestream >> filename;
283 ifstream stifffile(&filename[0]);
284 if (stifffile.is_open())
286 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
292 std::cout <<
"End Result: TEST FAILED\n";
298 for (
int cell_id = 3; cell_id < numCells; cell_id++) {
300 stringstream namestream;
302 namestream << basedir <<
"/mass_fp_HCURL_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
303 namestream >> filename;
305 ifstream massfile(&filename[0]);
306 if (massfile.is_open()) {
307 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
313 std::cout <<
"End Result: TEST FAILED\n";
318 namestream << basedir <<
"/stiff_fp_HCURL_HEX_I1_FEM" <<
"_" <<
"0" << cell_id+1 <<
".dat";
319 namestream >> filename;
321 ifstream stifffile(&filename[0]);
322 if (stifffile.is_open())
324 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
330 std::cout <<
"End Result: TEST FAILED\n";
340 catch (
const std::logic_error & err) {
341 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
342 *outStream << err.what() <<
'\n';
343 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
348 std::cout <<
"End Result: TEST FAILED\n";
350 std::cout <<
"End Result: TEST PASSED\n";
353 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::HCURL_HEX_I1_FEM class.