51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_GlobalMPISession.hpp"
56using 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"; \
72int 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_HCURL_HEX_In_FEM) |\n" \
95 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
96 <<
"| 2) Basis values for VALUE and CURL 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);
117 Basis_HCURL_HEX_In_FEM<double, FieldContainer<double> > hexBasis(deg,closedPts,openPts);
122 int throwCounter = 0;
125 FieldContainer<double> hexNodes(8, 3);
126 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
127 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
128 hexNodes(2,0) = -1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
129 hexNodes(3,0) = 1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
130 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
131 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
132 hexNodes(6,0) = -1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
133 hexNodes(7,0) = 1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
138 FieldContainer<double> vals;
142 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
143 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
147 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
148 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), 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
166 FieldContainer<double> badPoints1(4, 5, 3);
167 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
170 FieldContainer<double> badPoints2(4, 2);
171 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
174 FieldContainer<double> badVals1(4, 3);
175 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
178 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
181 FieldContainer<double> badVals2(hexBasis.getCardinality() + 1, hexNodes.dimension(0), 3);
182 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
185 FieldContainer<double> badVals3(hexBasis.getCardinality(), hexNodes.dimension(0) + 1, 3);
186 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
189 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
190 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
193 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
197 vals.resize(hexBasis.getCardinality(),
198 hexNodes.dimension(0),
199 Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension()));
200 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
205 catch (
const std::logic_error & err) {
206 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207 *outStream << err.what() <<
'\n';
208 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
214 if (throwCounter != nException) {
216 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
222 <<
"===============================================================================\n"\
223 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
224 <<
"===============================================================================\n";
227 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
230 for (
unsigned i = 0; i < allTags.size(); i++) {
231 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
235 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
236 if( !( (myTag[0] == allTags[i][0]) &&
237 (myTag[1] == allTags[i][1]) &&
238 (myTag[2] == allTags[i][2]) &&
239 (myTag[3] == allTags[i][3]) ) ) {
241 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
242 *outStream <<
" getDofOrdinal( {"
243 << allTags[i][0] <<
", "
244 << allTags[i][1] <<
", "
245 << allTags[i][2] <<
", "
246 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
247 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
251 << myTag[3] <<
"}\n";
256 for(
int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
257 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
258 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
259 if( bfOrd != myBfOrd) {
261 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
262 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
266 << myTag[3] <<
"} but getDofOrdinal({"
270 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
274 catch (
const std::logic_error & err){
275 *outStream << err.what() <<
"\n\n";
281 <<
"===============================================================================\n"\
282 <<
"| TEST 3: correctness of basis function values |\n"\
283 <<
"===============================================================================\n";
285 outStream -> precision(20);
288 double basisValues[] = {
289 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
290 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
291 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0,
292 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0,
293 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
294 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
295 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0,
296 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0,
297 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0,
298 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0,
299 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0,
300 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1
305 double basisCurls[] = {
306 0,-0.5,0.5, 0,-0.5,0.5, 0,0,0.5, 0,0,0.5, 0,-0.5,0, 0,-0.5,0, 0,0,0, 0,0,0,
307 0,0,-0.5, 0,0,-0.5, 0,-0.5,-0.5, 0,-0.5,-0.5, 0,0,0, 0,0,0, 0,-0.5,0, 0,-0.5,0,
308 0,0.5,0, 0,0.5,0, 0,0,0, 0,0,0, 0,0.5,0.5, 0,0.5,0.5, 0,0,0.5, 0,0,0.5,
309 0,0,0, 0,0,0, 0,0.5,0, 0,0.5,0, 0,0,-0.5, 0,0,-0.5, 0,0.5,-0.5, 0,0.5,-0.5,
313 0.5,0,-0.5, 0,0,-0.5, 0.5,0,-0.5, 0,0,-0.5, 0.5,0,0, 0,0,0, 0.5,0,0, 0,0,0,
316 0,0,0.5, 0.5,0,0.5, 0,0,0.5, 0.5,0,0.5, 0,0,0, 0.5,0,0, 0,0,0, 0.5,0,0,
319 -0.5,0,0, 0,0,0, -0.5,0,0, 0,0,0, -0.5,0,-0.5, 0,0,-0.5, -0.5,0,-0.5, 0,0,-0.5,
322 0.0,0,0, -0.5,0,0, 0.0,0,0, -0.5,0,0, 0.0,0,0.5, -0.5,0,0.5, 0.0,0,0.5, -0.5,0,0.5,
325 -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0, -0.5,0.5,0, 0,0.5,0, -0.5,0,0, 0,0,0,
328 0.0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0, 0,-0.5,0, -0.5,-0.5,0, 0,0,0, -0.5,0,0,
331 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0, 0.5,0,0, 0,0,0, 0.5,0.5,0, 0,0.5,0,
334 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0, 0,0,0, 0.5,0,0, 0,-0.5,0, 0.5,-0.5,0
341 int numFields = hexBasis.getCardinality();
342 int numPoints = hexNodes.dimension(0);
343 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
346 FieldContainer<double> vals;
349 vals.resize(numFields, numPoints, spaceDim);
350 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
351 for (
int i = 0; i < numFields; i++) {
352 for (
int j = 0; j < numPoints; j++) {
353 for (
int k = 0; k < spaceDim; k++) {
356 int l = k + i * spaceDim * numPoints + j * spaceDim;
357 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
359 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
362 *outStream <<
" At multi-index { ";
363 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
364 *outStream <<
"} computed value: " << vals(i,j,k)
365 <<
" but reference value: " << basisValues[l] <<
"\n";
372 vals.resize(numFields, numPoints, spaceDim);
373 hexBasis.getValues(vals, hexNodes, OPERATOR_CURL);
374 for (
int i = 0; i < numFields; i++) {
375 for (
int j = 0; j < numPoints; j++) {
376 for (
int k = 0; k < spaceDim; k++) {
378 int l = k + i * spaceDim * numPoints + j * spaceDim;
379 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
381 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
384 *outStream <<
" At multi-index { ";
385 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
386 *outStream <<
"} computed curl component: " << vals(i,j,k)
387 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
396 catch (
const std::logic_error & err) {
397 *outStream << err.what() <<
"\n\n";
402 std::cout <<
"End Result: TEST FAILED\n";
404 std::cout <<
"End Result: TEST PASSED\n";
407 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_HEX_In_FEM class.