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 (const 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_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";
115 int throwCounter = 0;
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;
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
152 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
156 INTREPID_TEST_COMMAND( triBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
160 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
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 );
174 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
178 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
182 INTREPID_TEST_COMMAND( triBasis.
getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
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";
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
302 int numPoints = triNodes.dimension(0);
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++) {
394 int DkCardin = Intrepid::getDkCardinality(op, spaceDim);
395 vals.
resize(numFields, numPoints, DkCardin);
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;
406 int ord = Intrepid::getOperatorOrder(op);
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 =
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);
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.
This is an interface class for bases whose degrees of freedom can be associated with spatial location...
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_C1_FEM class.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Header file for utility class to provide multidimensional containers.
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Triangle cell...
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 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...