49#include "Intrepid_HGRAD_WEDGE_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_WEDGE_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_WEDGE_C1_FEM<double, FieldContainer<double> > wedgeBasis;
114 int throwCounter = 0;
117 FieldContainer<double> wedgeNodes(12, 3);
118 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
119 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
120 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
121 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
122 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
123 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
125 wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
126 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
127 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.30; wedgeNodes(8,2) = 1.0;
128 wedgeNodes(9,0) = 0.3; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.44; wedgeNodes(10,2)= -0.23;
130 wedgeNodes(11,0)= 0.4; wedgeNodes(11,1)= 0.6; wedgeNodes(11,2)= 0.0;
134 FieldContainer<double> vals;
139 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
140 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
144 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
145 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
150 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
152 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
154 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,6,0), throwCounter, nException );
156 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(7), throwCounter, nException );
158 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
160#ifdef HAVE_INTREPID_DEBUG
163 FieldContainer<double> badPoints1(4, 5, 3);
164 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
167 FieldContainer<double> badPoints2(4, 2);
168 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
171 FieldContainer<double> badVals1(4, 3, 1);
172 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
175 FieldContainer<double> badVals2(4, 3);
176 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
179 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
182 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
185 FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
186 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
189 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
190 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
193 FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
194 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
197 FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
198 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
201 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
205 catch (
const std::logic_error & err) {
206 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207 *outStream << err.what() <<
'\n';
208 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
213 if (throwCounter != nException) {
215 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
220 <<
"===============================================================================\n"\
221 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
222 <<
"===============================================================================\n";
225 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
228 for (
unsigned i = 0; i < allTags.size(); i++) {
229 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
231 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
232 if( !( (myTag[0] == allTags[i][0]) &&
233 (myTag[1] == allTags[i][1]) &&
234 (myTag[2] == allTags[i][2]) &&
235 (myTag[3] == allTags[i][3]) ) ) {
237 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
238 *outStream <<
" getDofOrdinal( {"
239 << allTags[i][0] <<
", "
240 << allTags[i][1] <<
", "
241 << allTags[i][2] <<
", "
242 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
243 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
247 << myTag[3] <<
"}\n";
252 for(
int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
253 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
254 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
255 if( bfOrd != myBfOrd) {
257 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
258 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
262 << myTag[3] <<
"} but getDofOrdinal({"
266 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
270 catch (
const std::logic_error & err){
271 *outStream << err.what() <<
"\n\n";
277 <<
"===============================================================================\n"\
278 <<
"| TEST 3: correctness of basis function values |\n"\
279 <<
"===============================================================================\n";
281 outStream -> precision(20);
284 double basisValues[] = {
285 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
286 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
287 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
288 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
289 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
290 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
292 0.25, 0.25, 0.5, 0., 0., 0.,
293 0.125, 0.25, 0.125, 0.125, 0.25, 0.125,
294 0., 0., 0., 0.45, 0.25, 0.3,
295 0.0875, 0.0375, 0, 0.6125, 0.2625, 0,
296 0.3444, 0, 0.2706, 0.2156, 0, 0.1694,
297 0., 0.2, 0.3, 0., 0.2, 0.3
301 double basisGrads[] = {
302 -1., -1., -0.5, 1., 0, 0, 0, 1., 0, 0., 0., 0.5, 0., 0, 0, 0, 0., 0, \
303 -1., -1., 0., 1., 0, -0.5, 0, 1., 0, 0., 0., 0., 0., 0, 0.5, 0, 0., \
304 0, -1., -1., 0., 1., 0, 0, 0, 1., -0.5, 0., 0., 0., 0., 0, 0, 0, 0., \
305 0.5, 0., 0., -0.5, 0., 0, 0, 0, 0., 0, -1., -1., 0.5, 1., 0, 0, 0, \
306 1., 0, 0., 0., 0., 0., 0, -0.5, 0, 0., 0, -1., -1., 0., 1., 0, 0.5, \
307 0, 1., 0, 0., 0., 0., 0., 0, 0, 0, 0., -0.5, -1., -1., 0., 1., 0, 0, \
308 0, 1., 0.5, -1., -1., -0.125, 1., 0, -0.125, 0, 1., -0.25, 0., 0., \
309 0.125, 0., 0, 0.125, 0, 0., 0.25, -0.5, -0.5, -0.125, 0.5, 0, -0.25, \
310 0, 0.5, -0.125, -0.5, -0.5, 0.125, 0.5, 0, 0.25, 0, 0.5, 0.125, 0., \
311 0., -0.225, 0., 0, -0.125, 0, 0., -0.15, -1., -1., 0.225, 1., 0, \
312 0.125, 0, 1., 0.15, -0.125, -0.125, -0.35, 0.125, 0, -0.15, 0, 0.125, \
313 0, -0.875, -0.875, 0.35, 0.875, 0, 0.15, 0, 0.875, 0, -0.615, -0.615, \
314 -0.28, 0.615, 0, 0, 0, 0.615, -0.22, -0.385, -0.385, 0.28, 0.385, 0, \
315 0, 0, 0.385, 0.22, -0.5, -0.5, 0., 0.5, 0, -0.2, 0, 0.5, -0.3, -0.5, \
316 -0.5, 0., 0.5, 0, 0.2, 0, 0.5, 0.3
320 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
321 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
322 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
323 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
324 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
325 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
326 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
327 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, \
328 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, \
329 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, \
330 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, \
331 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, \
332 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, \
333 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, \
334 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, \
335 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
336 -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
337 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
338 -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
339 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
340 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
341 -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
342 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0
347 int numFields = wedgeBasis.getCardinality();
348 int numPoints = wedgeNodes.dimension(0);
349 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
353 FieldContainer<double> vals;
356 vals.resize(numFields, numPoints);
357 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
358 for (
int i = 0; i < numFields; i++) {
359 for (
int j = 0; j < numPoints; j++) {
360 int l = i + j * numFields;
361 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
363 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
366 *outStream <<
" At multi-index { ";
367 *outStream << i <<
" ";*outStream << j <<
" ";
368 *outStream <<
"} computed value: " << vals(i,j)
369 <<
" but reference value: " << basisValues[l] <<
"\n";
375 vals.resize(numFields, numPoints, spaceDim);
376 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
377 for (
int i = 0; i < numFields; i++) {
378 for (
int j = 0; j < numPoints; j++) {
379 for (
int k = 0; k < spaceDim; k++) {
380 int l = k + i * spaceDim + j * spaceDim * numFields;
381 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
383 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
386 *outStream <<
" At multi-index { ";
387 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
388 *outStream <<
"} computed grad component: " << vals(i,j,k)
389 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
396 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
397 for (
int i = 0; i < numFields; i++) {
398 for (
int j = 0; j < numPoints; j++) {
399 for (
int k = 0; k < spaceDim; k++) {
400 int l = k + i * spaceDim + j * spaceDim * numFields;
401 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
403 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
406 *outStream <<
" At multi-index { ";
407 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
408 *outStream <<
"} computed D1 component: " << vals(i,j,k)
409 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
416 vals.resize(numFields, numPoints, D2Cardin);
417 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
418 for (
int i = 0; i < numFields; i++) {
419 for (
int j = 0; j < numPoints; j++) {
420 for (
int k = 0; k < D2Cardin; k++) {
421 int l = k + i * D2Cardin + j * D2Cardin * numFields;
422 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
424 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
427 *outStream <<
" At multi-index { ";
428 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
429 *outStream <<
"} computed D2 component: " << vals(i,j,k)
430 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
437 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
441 vals.resize(numFields, numPoints, DkCardin);
443 wedgeBasis.getValues(vals, wedgeNodes, op);
444 for (
int i = 0; i < vals.size(); i++) {
445 if (std::abs(vals[i]) > INTREPID_TOL) {
447 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
450 std::vector<int> myIndex;
451 vals.getMultiIndex(myIndex,i);
453 *outStream <<
" At multi-index { ";
454 for(
int j = 0; j < vals.rank(); j++) {
455 *outStream << myIndex[j] <<
" ";
457 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
458 <<
" but reference D" << ord <<
" component: 0 \n";
465 catch (
const std::logic_error & err) {
466 *outStream << err.what() <<
"\n\n";
471 std::cout <<
"End Result: TEST FAILED\n";
473 std::cout <<
"End Result: TEST PASSED\n";
476 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.