51 #include "Teuchos_oblackholestream.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
56 using namespace Intrepid;
58 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
64 catch (const std::logic_error & err) { \
66 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
72 int main(
int argc,
char *argv[]) {
74 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
78 int iprint = argc - 1;
79 Teuchos::RCP<std::ostream> outStream;
80 Teuchos::oblackholestream bhs;
82 outStream = Teuchos::rcp(&std::cout,
false);
84 outStream = Teuchos::rcp(&bhs,
false);
87 Teuchos::oblackholestream oldFormatState;
88 oldFormatState.copyfmt(std::cout);
91 <<
"===============================================================================\n" \
93 <<
"| Unit Test (Basis_HDIV_HEX_In_FEM) |\n" \
95 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 <<
"| 2) Basis values for VALUE and DIV operators |\n" \
98 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
100 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
102 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
105 <<
"===============================================================================\n"\
106 <<
"| TEST 1: Basis creation, exception testing |\n"\
107 <<
"===============================================================================\n";
111 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());
114 PointTools::getLattice<double,FieldContainer<double> >(closedPts,line,deg);
115 PointTools::getLattice<double,FieldContainer<double> >(openPts,line,deg+1,1);
123 int throwCounter = 0;
127 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
128 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
129 hexNodes(2,0) = -1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
130 hexNodes(3,0) = 1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
131 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
132 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
133 hexNodes(6,0) = -1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
134 hexNodes(7,0) = 1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
144 vals.
resize(hexBasis.getCardinality(), hexNodes.dimension(0), 3 );
145 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
148 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
153 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
155 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
157 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
159 INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException );
161 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
163 #ifdef HAVE_INTREPID_DEBUG
167 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
171 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
175 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
179 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_DIV), throwCounter, nException );
183 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
187 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_DIV), throwCounter, nException );
191 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_VALUE), throwCounter, nException );
195 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_DIV), throwCounter, nException );
199 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_VALUE), throwCounter, nException );
203 catch (
const std::logic_error & err) {
204 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
205 *outStream << err.what() <<
'\n';
206 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
211 if (throwCounter != nException) {
213 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
218 <<
"===============================================================================\n"\
219 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
220 <<
"===============================================================================\n";
223 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
226 for (
unsigned i = 0; i < allTags.size(); i++) {
227 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
229 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
230 if( !( (myTag[0] == allTags[i][0]) &&
231 (myTag[1] == allTags[i][1]) &&
232 (myTag[2] == allTags[i][2]) &&
233 (myTag[3] == allTags[i][3]) ) ) {
235 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
236 *outStream <<
" getDofOrdinal( {"
237 << allTags[i][0] <<
", "
238 << allTags[i][1] <<
", "
239 << allTags[i][2] <<
", "
240 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
241 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
245 << myTag[3] <<
"}\n";
250 for(
int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
251 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
252 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
253 if( bfOrd != myBfOrd) {
255 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
256 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
260 << myTag[3] <<
"} but getDofOrdinal({"
264 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
268 catch (
const std::logic_error & err){
269 *outStream << err.what() <<
"\n\n";
275 <<
"===============================================================================\n"\
276 <<
"| TEST 3: correctness of basis function values |\n"\
277 <<
"===============================================================================\n";
279 outStream -> precision(20);
282 double basisValues[] = {
284 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
285 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
287 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,
288 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,
290 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
291 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
293 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
294 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
296 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1,
297 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
299 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
300 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1
304 double basisDivs[] = {
305 -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
306 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
307 -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
308 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
309 -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
310 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
316 int numPoints = hexNodes.dimension(0);
317 int numFields = hexBasis.getCardinality();
318 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
324 vals.
resize(numFields, numPoints, spaceDim);
325 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
326 for (
int i = 0; i < numFields; i++) {
327 for (
int j = 0; j < numPoints; j++) {
328 for (
int k = 0; k < spaceDim; k++) {
329 int l = k + i * numPoints * spaceDim + j * spaceDim;
330 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
332 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
335 *outStream <<
" At multi-index { ";
336 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
337 *outStream <<
"} computed value: " << vals(i,j,k)
338 <<
" but reference value: " << basisValues[l] <<
"\n";
345 vals.
resize(numFields, numPoints);
346 hexBasis.getValues(vals, hexNodes, OPERATOR_DIV);
347 for (
int i = 0; i < numFields; i++) {
348 for (
int j = 0; j < numPoints; j++) {
349 int l = i * numPoints + j;
350 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
352 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
355 *outStream <<
" At multi-index { ";
356 *outStream << i <<
" ";*outStream << j <<
" ";
357 *outStream <<
"} computed divergence component: " << vals(i,j)
358 <<
" but reference divergence component: " << basisDivs[l] <<
"\n";
366 catch (
const std::logic_error & err) {
367 *outStream << err.what() <<
"\n\n";
372 std::cout <<
"End Result: TEST FAILED\n";
374 std::cout <<
"End Result: TEST PASSED\n";
377 std::cout.copyfmt(oldFormatState);
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Header file for the Intrepid::HDIV_HEX_In_FEM class.
Header file for utility class to provide multidimensional containers.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...