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_C2_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_C2_FEM<double, FieldContainer<double> > triBasis;
115 int throwCounter = 0;
121 FieldContainer<double> triNodes(7, 2);
122 triNodes(0,0) = 0.0; triNodes(0,1) = 0.0;
123 triNodes(1,0) = 1.0; triNodes(1,1) = 0.0;
124 triNodes(2,0) = 0.0; triNodes(2,1) = 1.0;
125 triNodes(3,0) = 0.5; triNodes(3,1) = 0.0;
126 triNodes(4,0) = 0.5; triNodes(4,1) = 0.5;
127 triNodes(5,0) = 0.0; triNodes(5,1) = 0.5;
128 triNodes(6,0) = 0.1732; triNodes(6,1) = 0.6714;
131 FieldContainer<double> vals;
136 vals.resize(triBasis.getCardinality(), triNodes.dimension(0) );
137 INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
142 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException );
144 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException );
146 INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException );
148 INTREPID_TEST_COMMAND( triBasis.getDofTag(6), throwCounter, nException );
150 INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException );
152#ifdef HAVE_INTREPID_DEBUG
155 FieldContainer<double> badPoints1(4, 5, 3);
156 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
159 FieldContainer<double> badPoints2(4, triBasis.getBaseCellTopology().getDimension() + 1);
160 INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
163 FieldContainer<double> badVals1(4, 3, 1);
164 INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
167 FieldContainer<double> badVals2(4, 3);
168 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
171 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
174 INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
177 FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0));
178 INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
181 FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1);
182 INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
185 FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), triBasis.getBaseCellTopology().getDimension() + 1);
186 INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
189 FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40);
190 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
193 INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException );
197 catch (
const std::logic_error & err) {
198 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
199 *outStream << err.what() <<
'\n';
200 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
205 if (throwCounter != nException) {
207 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
212 <<
"===============================================================================\n"\
213 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
214 <<
"===============================================================================\n";
217 std::vector<std::vector<int> > allTags = triBasis.getAllDofTags();
220 for (
unsigned i = 0; i < allTags.size(); i++) {
221 int bfOrd = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
223 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
224 if( !( (myTag[0] == allTags[i][0]) &&
225 (myTag[1] == allTags[i][1]) &&
226 (myTag[2] == allTags[i][2]) &&
227 (myTag[3] == allTags[i][3]) ) ) {
229 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
230 *outStream <<
" getDofOrdinal( {"
231 << allTags[i][0] <<
", "
232 << allTags[i][1] <<
", "
233 << allTags[i][2] <<
", "
234 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
235 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
239 << myTag[3] <<
"}\n";
244 for(
int bfOrd = 0; bfOrd < triBasis.getCardinality(); bfOrd++) {
245 std::vector<int> myTag = triBasis.getDofTag(bfOrd);
246 int myBfOrd = triBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
247 if( bfOrd != myBfOrd) {
249 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
250 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
254 << myTag[3] <<
"} but getDofOrdinal({"
258 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
262 catch (
const std::logic_error & err){
263 *outStream << err.what() <<
"\n\n";
269 <<
"===============================================================================\n"\
270 <<
"| TEST 3: correctness of basis function values |\n"\
271 <<
"===============================================================================\n";
273 outStream -> precision(20);
277 double basisValues[] = {
278 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.10710168,
279 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -0.11320352,
280 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.23015592,
281 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.10766112,
282 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.46514592,
283 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.41734224
288 double basisGrads[] = {
290 -3.0, -3.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0,-1.0, 0.3784, 0.3784,
291 -1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, -0.3072, 0.0,
292 0.0, -1.0, 0.0,-1.0, 0.0, 3.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.6856,
293 4.0, 0.0, -4.0,-4.0, 0.0, 0.0, 0.0, -2.0, -2.0,-2.0, 2.0, 0.0, -0.0712, -0.6928,
294 0.0, 0.0, 0.0, 4.0, 4.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 2.6856, 0.6928,
295 0.0, 4.0, 0.0, 0.0, -4.0,-4.0, 0.0, 2.0, -2.0,-2.0, -2.0, 0.0, -2.6856, -2.064
301 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0,
302 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0,
303 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0,
304 -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0, -8.0,-4.0, 0.0,
305 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0, 0.0, 4.0, 0.0,
306 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0, 0.0,-4.0,-8.0
311 int numFields = triBasis.getCardinality();
312 int numPoints = triNodes.dimension(0);
313 int spaceDim = triBasis.getBaseCellTopology().getDimension();
316 FieldContainer<double> vals;
319 vals.resize(numFields, numPoints);
320 triBasis.getValues(vals, triNodes, OPERATOR_VALUE);
321 for (
int i = 0; i < numFields; i++) {
322 for (
int j = 0; j < numPoints; j++) {
325 int l = j + i * numPoints;
326 if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
328 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
331 *outStream <<
" At multi-index { ";
332 *outStream << i <<
" ";*outStream << j <<
" ";
333 *outStream <<
"} computed value: " << vals(i,j)
334 <<
" but reference value: " << basisValues[l] <<
"\n";
340 vals.resize(numFields, numPoints, spaceDim);
341 triBasis.getValues(vals, triNodes, OPERATOR_GRAD);
342 for (
int i = 0; i < numFields; i++) {
343 for (
int j = 0; j < numPoints; j++) {
344 for (
int k = 0; k < spaceDim; k++) {
346 int l = k + j * spaceDim + i * spaceDim * numPoints;
347 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
349 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
352 *outStream <<
" At multi-index { ";
353 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
354 *outStream <<
"} computed grad component: " << vals(i,j,k)
355 <<
" but reference grad component: " << basisGrads[l] <<
"\n";
362 triBasis.getValues(vals, triNodes, OPERATOR_D1);
363 for (
int i = 0; i < numFields; i++) {
364 for (
int j = 0; j < numPoints; j++) {
365 for (
int k = 0; k < spaceDim; k++) {
367 int l = k + j * spaceDim + i * spaceDim * numPoints;
368 if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
370 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
373 *outStream <<
" At multi-index { ";
374 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
375 *outStream <<
"} computed D1 component: " << vals(i,j,k)
376 <<
" but reference D1 component: " << basisGrads[l] <<
"\n";
383 vals.resize(numFields, numPoints, spaceDim);
384 triBasis.getValues(vals, triNodes, OPERATOR_CURL);
385 for (
int i = 0; i < numFields; i++) {
386 for (
int j = 0; j < numPoints; j++) {
388 int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints;
389 int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints;
391 double curl_value_0 = basisGrads[curl_0];
392 double curl_value_1 =-basisGrads[curl_1];
393 if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
395 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
397 *outStream <<
" At multi-index { ";
398 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << 0 <<
" ";
399 *outStream <<
"} computed curl component: " << vals(i,j,0)
400 <<
" but reference curl component: " << curl_value_0 <<
"\n";
402 if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
404 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
406 *outStream <<
" At multi-index { ";
407 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << 1 <<
" ";
408 *outStream <<
"} computed curl component: " << vals(i,j,1)
409 <<
" but reference curl component: " << curl_value_1 <<
"\n";
416 vals.resize(numFields, numPoints, D2cardinality);
417 triBasis.getValues(vals, triNodes, OPERATOR_D2);
418 for (
int i = 0; i < numFields; i++) {
419 for (
int j = 0; j < numPoints; j++) {
420 for (
int k = 0; k < D2cardinality; k++) {
423 int l = k + j * D2cardinality + i * D2cardinality * numPoints;
424 if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
426 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
429 *outStream <<
" At multi-index { ";
430 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
431 *outStream <<
"} computed D2 component: " << vals(i,j,k)
432 <<
" but reference D2 component: " << basisD2[l] <<
"\n";
439 for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
443 vals.resize(numFields, numPoints, DkCardin);
445 triBasis.getValues(vals, triNodes, op);
446 for (
int i = 0; i < vals.size(); i++) {
447 if (std::abs(vals[i]) > INTREPID_TOL) {
449 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
452 std::vector<int> myIndex;
453 vals.getMultiIndex(myIndex,i);
455 *outStream <<
" At multi-index { ";
456 for(
int j = 0; j < vals.rank(); j++) {
457 *outStream << myIndex[j] <<
" ";
459 *outStream <<
"} computed D"<< ord <<
" component: " << vals[i]
460 <<
" but reference D" << ord <<
" component: 0 \n";
467 catch (
const std::logic_error & err) {
468 *outStream << err.what() <<
"\n\n";
473 std::cout <<
"End Result: TEST FAILED\n";
475 std::cout <<
"End Result: TEST PASSED\n";
478 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TRI_C2_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.