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"; \
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_TRI_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 (dridzal@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";
110 Basis_HGRAD_TRI_C1_FEM<double, FieldContainer<double> > triBasis;
115 int throwCounter = 0;
118 FieldContainer<double> triNodes(5, 2);
119 triNodes(0,0) = 0.0; triNodes(0,1) = 0.0;
120 triNodes(1,0) = 1.0; triNodes(1,1) = 0.0;
121 triNodes(2,0) = 0.0; triNodes(2,1) = 1.0;
122 triNodes(3,0) = 0.5; triNodes(3,1) = 0.5;
123 triNodes(4,0) = 0.0; triNodes(4,1) = 0.75;
128 FieldContainer<double> vals;
132 vals.resize(triBasis.getCardinality(), triNodes.dimension(0) );
133 INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
138 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException );
140 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException );
142 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException );
144 INTREPID_TEST_COMMAND( triBasis.getDofTag(5), throwCounter, nException );
146 INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException );
148#ifdef HAVE_INTREPID_DEBUG
151 FieldContainer<double> badPoints1(4, 5, 3);
152 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
155 FieldContainer<double> badPoints2(4, 3);
156 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
159 FieldContainer<double> badVals1(4, 3, 1);
160 INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
163 FieldContainer<double> badVals2(4, 3);
164 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
167 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
170 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
173 FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0));
174 INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
177 FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1);
178 INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
181 FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), 4);
182 INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
185 FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40);
186 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
189 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException );
193 catch (
const std::logic_error & err) {
194 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
195 *outStream << err.what() <<
'\n';
196 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
201 if (throwCounter != nException) {
203 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
208 <<
"===============================================================================\n"\
209 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
210 <<
"===============================================================================\n";
213 std::vector<std::vector<int> > allTags = triBasis.getAllDofTags();
216 for (
unsigned i = 0; i < allTags.size(); i++) {
217 int bfOrd = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
219 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
220 if( !( (myTag[0] == allTags[i][0]) &&
221 (myTag[1] == allTags[i][1]) &&
222 (myTag[2] == allTags[i][2]) &&
223 (myTag[3] == allTags[i][3]) ) ) {
225 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
226 *outStream <<
" getDofOrdinal( {"
227 << allTags[i][0] <<
", "
228 << allTags[i][1] <<
", "
229 << allTags[i][2] <<
", "
230 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
231 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
235 << myTag[3] <<
"}\n";
240 for(
int bfOrd = 0; bfOrd < triBasis.getCardinality(); bfOrd++) {
241 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
242 int myBfOrd = triBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
243 if( bfOrd != myBfOrd) {
245 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
246 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
250 << myTag[3] <<
"} but getDofOrdinal({"
254 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
258 catch (
const std::logic_error & err){
259 *outStream << err.what() <<
"\n\n";
265 <<
"===============================================================================\n"\
266 <<
"| TEST 3: correctness of basis function values |\n"\
267 <<
"===============================================================================\n";
269 outStream -> precision(20);
272 double basisValues[] = {
281 double basisGrads[] = {
282 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
283 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
284 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
285 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
286 -1.0, -1.0, 1.0, 0.0, 0.0, 1.0,
290 double basisCurls[] = {
291 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
292 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
293 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
294 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0,
295 -1.0, 1.0, 0.0, -1.0, 1.0, 0.0
301 int numFields = triBasis.getCardinality();
302 int numPoints = triNodes.dimension(0);
303 int spaceDim = triBasis.getBaseCellTopology().getDimension();
306 FieldContainer<double> vals;
309 vals.resize(numFields, numPoints);
310 triBasis.getValues(vals, triNodes, OPERATOR_VALUE);
311 for (
int i = 0; i < numFields; i++) {
312 for (
int j = 0; j < numPoints; j++) {
313 int l = i + j * numFields;
314 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
316 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
319 *outStream <<
" At multi-index { ";
320 *outStream << i <<
" ";*outStream << j <<
" ";
321 *outStream <<
"} computed value: " << vals(i,j)
322 <<
" but reference value: " << basisValues[l] <<
"\n";
328 vals.resize(numFields, numPoints, spaceDim);
329 triBasis.getValues(vals, triNodes, OPERATOR_GRAD);
330 for (
int i = 0; i < numFields; i++) {
331 for (
int j = 0; j < numPoints; j++) {
332 for (
int k = 0; k < spaceDim; k++) {
333 int l = k + i * spaceDim + j * spaceDim * numFields;
334 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
336 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
339 *outStream <<
" At multi-index { ";
340 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
341 *outStream <<
"} computed grad component: " << vals(i,j,k)
342 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
349 triBasis.getValues(vals, triNodes, OPERATOR_D1);
350 for (
int i = 0; i < numFields; i++) {
351 for (
int j = 0; j < numPoints; j++) {
352 for (
int k = 0; k < spaceDim; k++) {
353 int l = k + i * spaceDim + j * spaceDim * numFields;
354 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
356 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
359 *outStream <<
" At multi-index { ";
360 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
361 *outStream <<
"} computed D1 component: " << vals(i,j,k)
362 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
370 vals.resize(numFields, numPoints, spaceDim);
371 triBasis.getValues(vals, triNodes, OPERATOR_CURL);
372 for (
int i = 0; i < numFields; i++) {
373 for (
int j = 0; j < numPoints; j++) {
374 for (
int k = 0; k < spaceDim; k++) {
375 int l = k + i * spaceDim + j * spaceDim * numFields;
376 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
378 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
381 *outStream <<
" At multi-index { ";
382 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
383 *outStream <<
"} computed curl component: " << vals(i,j,k)
384 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
391 for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
395 vals.resize(numFields, numPoints, DkCardin);
397 triBasis.getValues(vals, triNodes, op);
398 for (
int i = 0; i < vals.size(); i++) {
399 if (std::abs(vals[i]) > INTREPID_TOL) {
401 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
404 std::vector<int> myIndex;
405 vals.getMultiIndex(myIndex,i);
407 *outStream <<
" At multi-index { ";
408 for(
int j = 0; j < vals.rank(); j++) {
409 *outStream << myIndex[j] <<
" ";
411 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
412 <<
" but reference D" << ord <<
" component: 0 \n";
419 catch (
const std::logic_error & err) {
420 *outStream << err.what() <<
"\n\n";
426 <<
"===============================================================================\n"\
427 <<
"| TEST 4: correctness of DoF locations |\n"\
428 <<
"===============================================================================\n";
431 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
433 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
434 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
436 FieldContainer<double> cvals;
437 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
440#ifdef HAVE_INTREPID_DEBUG
442 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
444 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
446 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
449 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
451 if (throwCounter != nException) {
453 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
457 basis->getValues(bvals, cvals, OPERATOR_VALUE);
459 for (
int i=0; i<bvals.dimension(0); i++) {
460 for (
int j=0; j<bvals.dimension(1); j++) {
461 if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
463 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 0.0);
464 *outStream << buffer;
466 else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
468 sprintf(buffer,
"\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 1.0);
469 *outStream << buffer;
475 catch (
const std::logic_error & err){
476 *outStream << err.what() <<
"\n\n";
481 std::cout <<
"End Result: TEST FAILED\n";
483 std::cout <<
"End Result: TEST PASSED\n";
486 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TRI_C1_FEM class.
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 Triangle cell.