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_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);
122 int throwCounter = 0;
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;
143 vals.
resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
144 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD), throwCounter, nException );
148 vals.
resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
149 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
154 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
156 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
158 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
160 INTREPID_TEST_COMMAND( hexBasis.getDofTag(12), throwCounter, nException );
162 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
164 #ifdef HAVE_INTREPID_DEBUG
168 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
172 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
176 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
179 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_CURL), throwCounter, nException );
183 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
187 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
191 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException ) ;
194 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_CURL), throwCounter, nException ) ;
198 vals.
resize(hexBasis.getCardinality(),
199 hexNodes.dimension(0),
200 Intrepid::getDkCardinality(OPERATOR_D2, hexBasis.getBaseCellTopology().getDimension()));
201 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_D2), throwCounter, nException );
206 catch (
const std::logic_error & err) {
207 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
208 *outStream << err.what() <<
'\n';
209 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
215 if (throwCounter != nException) {
217 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
223 <<
"===============================================================================\n"\
224 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
225 <<
"===============================================================================\n";
228 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
231 for (
unsigned i = 0; i < allTags.size(); i++) {
232 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
236 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
237 if( !( (myTag[0] == allTags[i][0]) &&
238 (myTag[1] == allTags[i][1]) &&
239 (myTag[2] == allTags[i][2]) &&
240 (myTag[3] == allTags[i][3]) ) ) {
242 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
243 *outStream <<
" getDofOrdinal( {"
244 << allTags[i][0] <<
", "
245 << allTags[i][1] <<
", "
246 << allTags[i][2] <<
", "
247 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
248 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
252 << myTag[3] <<
"}\n";
257 for(
int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
258 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
259 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
260 if( bfOrd != myBfOrd) {
262 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
263 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
267 << myTag[3] <<
"} but getDofOrdinal({"
271 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
275 catch (
const std::logic_error & err){
276 *outStream << err.what() <<
"\n\n";
282 <<
"===============================================================================\n"\
283 <<
"| TEST 3: correctness of basis function values |\n"\
284 <<
"===============================================================================\n";
286 outStream -> precision(20);
289 double basisValues[] = {
290 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
291 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
292 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0, 0,0,0, 0,0,0,
293 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 1,0,0, 1,0,0,
294 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
295 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0,
296 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0, 0,0,0,
297 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,1,0, 0,0,0, 0,1,0,
298 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0,
299 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0,
300 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0,
301 0,0,0, 0,0,0, 0,0,0, 0,0,1, 0,0,0, 0,0,0, 0,0,0, 0,0,1
306 double basisCurls[] = {
307 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,
308 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,
309 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,
310 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,
314 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,
317 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,
320 -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,
323 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,
326 -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,
329 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,
332 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,
335 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
342 int numFields = hexBasis.getCardinality();
343 int numPoints = hexNodes.dimension(0);
344 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
350 vals.
resize(numFields, numPoints, spaceDim);
351 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
352 for (
int i = 0; i < numFields; i++) {
353 for (
int j = 0; j < numPoints; j++) {
354 for (
int k = 0; k < spaceDim; k++) {
357 int l = k + i * spaceDim * numPoints + j * spaceDim;
358 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
360 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
363 *outStream <<
" At multi-index { ";
364 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
365 *outStream <<
"} computed value: " << vals(i,j,k)
366 <<
" but reference value: " << basisValues[l] <<
"\n";
373 vals.
resize(numFields, numPoints, spaceDim);
374 hexBasis.getValues(vals, hexNodes, OPERATOR_CURL);
375 for (
int i = 0; i < numFields; i++) {
376 for (
int j = 0; j < numPoints; j++) {
377 for (
int k = 0; k < spaceDim; k++) {
379 int l = k + i * spaceDim * numPoints + j * spaceDim;
380 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
382 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
385 *outStream <<
" At multi-index { ";
386 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
387 *outStream <<
"} computed curl component: " << vals(i,j,k)
388 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
397 catch (
const std::logic_error & err) {
398 *outStream << err.what() <<
"\n\n";
403 std::cout <<
"End Result: TEST FAILED\n";
405 std::cout <<
"End Result: TEST PASSED\n";
408 std::cout.copyfmt(oldFormatState);
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...
Implementation of the default H(div)-compatible FEM basis of degree 1 on Hexahedral cell...
Header file for the Intrepid::HCURL_HEX_In_FEM class.