49#include "Intrepid_HGRAD_HEX_C1_FEM.hpp"
50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
55using namespace Intrepid;
57#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
63 catch (const std::logic_error & err) { \
65 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
66 *outStream << err.what() << '\n'; \
67 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
71int main(
int argc,
char *argv[]) {
73 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77 int iprint = argc - 1;
78 Teuchos::RCP<std::ostream> outStream;
79 Teuchos::oblackholestream bhs;
81 outStream = Teuchos::rcp(&std::cout,
false);
83 outStream = Teuchos::rcp(&bhs,
false);
86 Teuchos::oblackholestream oldFormatState;
87 oldFormatState.copyfmt(std::cout);
90 <<
"===============================================================================\n" \
92 <<
"| Unit Test (Basis_HGRAD_HEX_C1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE, GRAD, CURL, and Dk operators |\n" \
97 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
98 <<
"| Denis Ridzal (dridzal@sandia.gov), |\n" \
99 <<
"| Kara Peterson (kjpeter@sandia.gov). |\n" \
101 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
104 <<
"===============================================================================\n"\
105 <<
"| TEST 1: Basis creation, exception testing |\n"\
106 <<
"===============================================================================\n";
109 Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> > hexBasis;
114 int throwCounter = 0;
117 FieldContainer<double> hexNodes(15, 3);
118 hexNodes(0,0) = -1.0; hexNodes(0,1) = -1.0; hexNodes(0,2) = -1.0;
119 hexNodes(1,0) = 1.0; hexNodes(1,1) = -1.0; hexNodes(1,2) = -1.0;
120 hexNodes(2,0) = 1.0; hexNodes(2,1) = 1.0; hexNodes(2,2) = -1.0;
121 hexNodes(3,0) = -1.0; hexNodes(3,1) = 1.0; hexNodes(3,2) = -1.0;
123 hexNodes(4,0) = -1.0; hexNodes(4,1) = -1.0; hexNodes(4,2) = 1.0;
124 hexNodes(5,0) = 1.0; hexNodes(5,1) = -1.0; hexNodes(5,2) = 1.0;
125 hexNodes(6,0) = 1.0; hexNodes(6,1) = 1.0; hexNodes(6,2) = 1.0;
126 hexNodes(7,0) = -1.0; hexNodes(7,1) = 1.0; hexNodes(7,2) = 1.0;
128 hexNodes(8,0) = 0.0; hexNodes(8,1) = 0.0; hexNodes(8,2) = 0.0;
130 hexNodes(9,0) = 1.0; hexNodes(9,1) = 0.0; hexNodes(9,2) = 0.0;
131 hexNodes(10,0)= -1.0; hexNodes(10,1)= 0.0; hexNodes(10,2)= 0.0;
133 hexNodes(11,0)= 0.0; hexNodes(11,1)= 1.0; hexNodes(11,2)= 0.0;
134 hexNodes(12,0)= 0.0; hexNodes(12,1)= -1.0; hexNodes(12,2)= 0.0;
136 hexNodes(13,0)= 0.0; hexNodes(13,1)= 0.0; hexNodes(13,2)= 1.0;
137 hexNodes(14,0)= 0.0; hexNodes(14,1)= 0.0; hexNodes(14,2)= -1.0;
141 FieldContainer<double> vals;
146 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
147 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
151 vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
152 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
157 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
159 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
161 INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
163 INTREPID_TEST_COMMAND( hexBasis.getDofTag(8), throwCounter, nException );
165 INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
167#ifdef HAVE_INTREPID_DEBUG
170 FieldContainer<double> badPoints1(4, 5, 3);
171 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
174 FieldContainer<double> badPoints2(4, 2);
175 INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
178 FieldContainer<double> badVals1(4, 3, 1);
179 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
182 FieldContainer<double> badVals2(4, 3);
183 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException );
186 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException );
189 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException );
192 FieldContainer<double> badVals3(hexBasis.getCardinality() + 1, hexNodes.dimension(0));
193 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
196 FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0) + 1);
197 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException );
200 FieldContainer<double> badVals5(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
201 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException );
204 FieldContainer<double> badVals6(hexBasis.getCardinality(), hexNodes.dimension(0), 40);
205 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException );
208 FieldContainer<double> badVals7(hexBasis.getCardinality(), hexNodes.dimension(0), 50);
209 INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException );
213 catch (
const std::logic_error & err) {
214 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
215 *outStream << err.what() <<
'\n';
216 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
222 if (throwCounter != nException) {
224 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
229 <<
"===============================================================================\n"\
230 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
231 <<
"===============================================================================\n";
234 std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
237 for (
unsigned i = 0; i < allTags.size(); i++) {
238 int bfOrd = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
240 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
241 if( !( (myTag[0] == allTags[i][0]) &&
242 (myTag[1] == allTags[i][1]) &&
243 (myTag[2] == allTags[i][2]) &&
244 (myTag[3] == allTags[i][3]) ) ) {
246 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
247 *outStream <<
" getDofOrdinal( {"
248 << allTags[i][0] <<
", "
249 << allTags[i][1] <<
", "
250 << allTags[i][2] <<
", "
251 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
252 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
256 << myTag[3] <<
"}\n";
261 for(
int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
262 std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
263 int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
264 if( bfOrd != myBfOrd) {
266 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
267 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
271 << myTag[3] <<
"} but getDofOrdinal({"
275 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
279 catch (
const std::logic_error & err){
280 *outStream << err.what() <<
"\n\n";
286 <<
"===============================================================================\n"\
287 <<
"| TEST 3: correctness of basis function values |\n"\
288 <<
"===============================================================================\n";
290 outStream -> precision(20);
293 double basisValues[] = {
295 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
296 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
297 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
298 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
300 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
301 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
302 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
303 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
305 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125,
307 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0,
308 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25,
310 0.0, 0.0, 0.25, 0.25, 0.0, 0.0, 0.25, 0.25,
311 0.25, 0.25, 0.0, 0.0, 0.25, 0.25, 0.0, 0.0,
313 0.0, 0.0, 0.0, 0.0, 0.25, 0.25, 0.25, 0.25,
314 0.25, 0.25, 0.25, 0.25, 0.0, 0.0, 0.0, 0.0,
318 double basisGrads[] = {
320 -0.5,-0.5,-0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
321 -0.5, 0.0, 0.0, 0.5,-0.5,-0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
322 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5,-0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,
323 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5,
325 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5,-0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0,
326 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.5, 0.0, 0.0, 0.5,-0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0,
327 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0, 0.5, 0.5, 0.5, -0.5, 0.0, 0.0,
328 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-0.5, 0.0,-0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.5, 0.5, 0.5,
330 -0.125,-0.125,-0.125, 0.125,-0.125,-0.125, 0.125, 0.125,-0.125, \
331 -0.125, 0.125,-0.125, -0.125,-0.125, 0.125, 0.125,-0.125, 0.125, \
332 0.125, 0.125, 0.125, -0.125, 0.125, 0.125,
334 -0.125, 0.0, 0.0, 0.125,-0.25, -0.25, 0.125, 0.25, -0.25, -0.125, 0.0, 0.0, \
335 -0.125, 0.0, 0.0, 0.125,-0.25, 0.25, 0.125, 0.25, 0.25, -0.125, 0.0, 0.0,
337 -0.125,-0.25, -0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, -0.25,\
338 -0.125,-0.25, 0.25, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, -0.125, 0.25, 0.25,
340 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125,-0.25, -0.25, 0.125,-0.25,\
341 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.25, 0.125, 0.25, -0.25, 0.125, 0.25,
343 -0.25, -0.125,-0.25, 0.25, -0.125,-0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, \
344 -0.25, -0.125, 0.25, 0.25, -0.125, 0.25, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0,
346 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, 0.0, 0.0, -0.125, \
347 -0.25, -0.25, 0.125, 0.25, -0.25, 0.125, 0.25, 0.25, 0.125, -0.25, 0.25, 0.125,
349 -0.25, -0.25, -0.125, 0.25, -0.25, -0.125, 0.25, 0.25, -0.125, -0.25, 0.25, -0.125, \
350 0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125, 0.0, 0.0, 0.125
356 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0, 0, 0.25, 0., 0, \
357 0., 0, 0, -0.25, 0., 0, -0.25, 0, 0, 0., -0.25, 0, -0.25, 0, 0, 0., \
358 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0., \
360 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0, 0, 0.25, 0., 0, \
361 -0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, 0, 0., \
362 0.25, 0, -0.25, 0, 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0., \
364 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0, 0, 0.25, -0.25, 0, \
365 -0.25, 0, 0, -0.25, 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., \
366 0, -0.25, 0, 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0., \
368 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0.25, -0.25, 0, \
369 0., 0, 0, -0.25, 0.25, 0, -0.25, 0, 0, 0., 0., 0, -0.25, 0, 0, 0., \
370 0., 0, 0., 0, 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0.,\
372 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, \
373 0, 0., 0., 0, -0.25, 0, 0, 0.25, -0.25, 0, -0.25, 0, 0, -0.25, 0.25, \
374 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0., \
376 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0, 0, 0., 0., 0, -0.25, \
377 0, 0, 0., 0., 0, 0., 0, 0, 0.25, -0.25, 0, 0., 0, 0, -0.25, 0.25, 0, \
378 -0.25, 0, 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0., \
380 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0, 0, 0., -0.25, 0, -0.25, \
381 0, 0, 0., 0.25, 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, \
382 -0.25, 0, 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0., \
384 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, \
385 0, 0., 0.25, 0, -0.25, 0, 0, 0.25, 0., 0, -0.25, 0, 0, -0.25, 0., 0, \
386 0., 0, 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0., \
388 0, 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0, 0, \
389 0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
390 0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
391 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0., \
393 0, 0.125, 0.125, 0, 0., 0, 0, -0.125, -0.125, 0, 0.25, 0, 0, 0.125, \
394 -0.125, 0, -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, -0.125, 0, \
395 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, 0.125, 0, 0.25, 0, 0, \
396 -0.125, -0.125, 0, 0., 0., \
398 0, 0.125, 0.125, 0, 0.25, 0, 0, -0.125, -0.125, 0, 0., 0, 0, 0.125, \
399 -0.125, 0, 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, -0.125, 0, \
400 -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, 0.125, 0, 0., 0, 0, \
401 -0.125, -0.125, 0, 0.25, 0., \
403 0, 0.125, 0., 0, 0.125, 0, 0, -0.125, 0., 0, 0.125, 0, 0, 0.125, \
404 -0.25, 0, -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, \
405 -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, 0.25, 0, 0.125, 0, \
406 0, -0.125, -0.25, 0, 0.125, 0., \
408 0, 0.125, 0.25, 0, 0.125, 0, 0, -0.125, -0.25, 0, 0.125, 0, 0, 0.125, \
409 0., 0, -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, -0.25, 0, \
410 -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, 0.125, 0, \
411 0, -0.125, 0., 0, 0.125, 0., \
413 0, 0., 0.125, 0, 0.125, 0, 0, 0., -0.125, 0, 0.125, 0, 0, 0., -0.125, \
414 0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0.25, -0.125, 0, -0.125, \
415 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0.25, 0.125, 0, 0.125, 0, 0, \
416 -0.25, -0.125, 0, 0.125, 0., \
418 0, 0.25, 0.125, 0, 0.125, 0, 0, -0.25, -0.125, 0, 0.125, 0, 0, 0.25, \
419 -0.125, 0, -0.125, 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0., -0.125, \
420 0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0., 0.125, 0, 0.125, 0, \
421 0, 0., -0.125, 0, 0.125, 0.
427 int numFields = hexBasis.getCardinality();
428 int numPoints = hexNodes.dimension(0);
429 int spaceDim = hexBasis.getBaseCellTopology().getDimension();
433 FieldContainer<double> vals;
436 vals.resize(numFields, numPoints);
437 hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
438 for (
int i = 0; i < numFields; i++) {
439 for (
int j = 0; j < numPoints; j++) {
440 int l = i + j * numFields;
441 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
443 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
446 *outStream <<
" At multi-index { ";
447 *outStream << i <<
" ";*outStream << j <<
" ";
448 *outStream <<
"} computed value: " << vals(i,j)
449 <<
" but reference value: " << basisValues[l] <<
"\n";
455 vals.resize(numFields, numPoints, spaceDim);
456 hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD);
457 for (
int i = 0; i < numFields; i++) {
458 for (
int j = 0; j < numPoints; j++) {
459 for (
int k = 0; k < spaceDim; k++) {
460 int l = k + i * spaceDim + j * spaceDim * numFields;
461 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
463 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
466 *outStream <<
" At multi-index { ";
467 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
468 *outStream <<
"} computed grad component: " << vals(i,j,k)
469 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
476 hexBasis.getValues(vals, hexNodes, OPERATOR_D1);
477 for (
int i = 0; i < numFields; i++) {
478 for (
int j = 0; j < numPoints; j++) {
479 for (
int k = 0; k < spaceDim; k++) {
480 int l = k + i * spaceDim + j * spaceDim * numFields;
481 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
483 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
486 *outStream <<
" At multi-index { ";
487 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
488 *outStream <<
"} computed D1 component: " << vals(i,j,k)
489 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
497 vals.resize(numFields, numPoints, D2Cardin);
498 hexBasis.getValues(vals, hexNodes, OPERATOR_D2);
499 for (
int i = 0; i < numFields; i++) {
500 for (
int j = 0; j < numPoints; j++) {
501 for (
int k = 0; k < D2Cardin; k++) {
502 int l = k + i * D2Cardin + j * D2Cardin * numFields;
503 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
505 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
508 *outStream <<
" At multi-index { ";
509 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
510 *outStream <<
"} computed D2 component: " << vals(i,j,k)
511 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
518 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
522 vals.resize(numFields, numPoints, DkCardin);
524 hexBasis.getValues(vals, hexNodes, op);
525 for (
int i = 0; i < vals.size(); i++) {
526 if (std::abs(vals[i]) > INTREPID_TOL) {
528 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
531 std::vector<int> myIndex;
532 vals.getMultiIndex(myIndex,i);
534 *outStream <<
" At multi-index { ";
535 for(
int j = 0; j < vals.rank(); j++) {
536 *outStream << myIndex[j] <<
" ";
538 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
539 <<
" but reference D" << ord <<
" component: 0 \n";
546 catch (
const std::logic_error & err) {
547 *outStream << err.what() <<
"\n\n";
553 <<
"===============================================================================\n"\
554 <<
"| TEST 4: correctness of DoF locations |\n"\
555 <<
"===============================================================================\n";
558 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
560 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
561 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
563 FieldContainer<double> cvals;
564 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
567#ifdef HAVE_INTREPID_DEBUG
569 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
571 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
573 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
576 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
578 if (throwCounter != nException) {
580 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
584 basis->getValues(bvals, cvals, OPERATOR_VALUE);
586 for (
int i=0; i<bvals.dimension(0); i++) {
587 for (
int j=0; j<bvals.dimension(1); j++) {
588 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
590 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 0.0);
591 *outStream << buffer;
593 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
595 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 1.0);
596 *outStream << buffer;
602 catch (
const std::logic_error & err){
603 *outStream << err.what() <<
"\n\n";
608 std::cout <<
"End Result: TEST FAILED\n";
610 std::cout <<
"End Result: TEST PASSED\n";
613 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
int getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Hexahedron cell.