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 (std::logic_error err) { \
66 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
71 int 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";
115 int throwCounter = 0;
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;
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
156 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
160 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
164 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
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 );
178 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
182 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
186 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
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 (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";
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 (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
312 int numPoints = triNodes.dimension(0);
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";
415 int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
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++) {
442 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
443 vals.
resize(numFields, numPoints, DkCardin);
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;
454 int ord = Intrepid::getOperatorOrder(op);
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 (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);
int size() const
Returns size of the FieldContainer defined as the product of its dimensions.
virtual const std::vector< int > & getDofTag(const int dofOrd)
DoF ordinal to DoF tag lookup.
virtual int getCardinality() const
Returns cardinality of the basis.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
void getMultiIndex(int &i0, const int valueEnum) const
Returns the multi-index of a value, based on its enumeration, as a list, for rank-1 containers...
int rank() const
Return rank of the FieldContainer = number of indices used to tag the multi-indexed value...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Header file for the Intrepid::HGRAD_TRI_C2_FEM class.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
virtual int getDofOrdinal(const int subcDim, const int subcOrd, const int subcDofOrd)
DoF tag to ordinal lookup.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
Implementation of the default H(grad)-compatible FEM basis of degree 2 on Triangle cell...