50 #include "Teuchos_oblackholestream.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
55 using namespace Intrepid;
57 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
63 catch (std::logic_error err) { \
65 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
66 *outStream << err.what() << '\n'; \
67 *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_HDIV_WEDGE_I1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE and DIV 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";
114 int throwCounter = 0;
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.25; wedgeNodes(8,2) = 1.0;
128 wedgeNodes(9,0) = 0.25; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
129 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.5; wedgeNodes(10,2)= -0.25;
130 wedgeNodes(11,0)= 0.5; wedgeNodes(11,1)= 0.5; wedgeNodes(11,2)= 0.0;
139 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
142 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_CURL), throwCounter, nException );
147 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(3,0,0), throwCounter, nException );
149 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(1,1,1), throwCounter, nException );
151 INTREPID_TEST_COMMAND( wedgeBasis.
getDofOrdinal(0,4,1), throwCounter, nException );
153 INTREPID_TEST_COMMAND( wedgeBasis.
getDofTag(11), throwCounter, nException );
155 INTREPID_TEST_COMMAND( wedgeBasis.
getDofTag(-1), throwCounter, nException );
157 #ifdef HAVE_INTREPID_DEBUG
161 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
165 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
169 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
173 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals2, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
177 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
181 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals4, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
185 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals5, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
189 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals6, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
193 INTREPID_TEST_COMMAND( wedgeBasis.
getValues(badVals7, wedgeNodes, OPERATOR_VALUE), 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 = wedgeBasis.
getAllDofTags();
220 for (
unsigned i = 0; i < allTags.size(); i++) {
221 int bfOrd = wedgeBasis.
getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
223 std::vector<int> myTag = wedgeBasis.
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 < wedgeBasis.
getCardinality(); bfOrd++) {
245 std::vector<int> myTag = wedgeBasis.
getDofTag(bfOrd);
246 int myBfOrd = wedgeBasis.
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);
276 double basisValues[] = {
277 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -2.00000, 0, 0, 0, \
278 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, -2.00000, 0, \
279 0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, 0, 0, \
280 -2.00000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, \
281 0, 0, 0, 2.00000, 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, \
282 0, 0, 0, 0, 2.00000, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, \
283 0, 0, 0, 0, 0, 2.00000, 0.125000, -0.250000, 0, 0.125000, 0.250000, \
284 0, -0.375000, 0.250000, 0, 0, 0, -2.00000, 0, 0, 0, 0.250000, \
285 -0.375000, 0, 0.250000, 0.125000, 0, -0.250000, 0.125000, 0, 0, 0, \
286 -1.00000, 0, 0, 1.00000, 0.125000, -0.375000, 0, 0.125000, 0.125000, \
287 0, -0.375000, 0.125000, 0, 0, 0, 0, 0, 0, 2.00000, 0.125000, \
288 -0.500000, 0, 0.125000, 0, 0, -0.375000, 0, 0, 0, 0, -0.250000, 0, 0, \
289 1.75000, 0, -0.250000, 0, 0, 0.250000, 0, -0.500000, 0.250000, 0, 0, \
290 0, -1.25000, 0, 0, 0.750000, 0.250000, -0.250000, 0, 0.250000, \
291 0.250000, 0, -0.250000, 0.250000, 0, 0, 0, -1.00000, 0, 0, 1.00000};
294 double basisDivs[] = {
296 1.0, 1.0, 1.0, 1.0, 1.0,
297 1.0, 1.0, 1.0, 1.0, 1.0,
298 1.0, 1.0, 1.0, 1.0, 1.0,
299 1.0, 1.0, 1.0, 1.0, 1.0,
300 1.0, 1.0, 1.0, 1.0, 1.0,
301 1.0, 1.0, 1.0, 1.0, 1.0,
303 1.0, 1.0, 1.0, 1.0, 1.0,
304 1.0, 1.0, 1.0, 1.0, 1.0,
305 1.0, 1.0, 1.0, 1.0, 1.0,
306 1.0, 1.0, 1.0, 1.0, 1.0,
307 1.0, 1.0, 1.0, 1.0, 1.0,
308 1.0, 1.0, 1.0, 1.0, 1.0
315 int numPoints = wedgeNodes.dimension(0);
322 vals.
resize(numFields, numPoints, spaceDim);
323 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_VALUE);
324 for (
int i = 0; i < numFields; i++) {
325 for (
int j = 0; j < numPoints; j++) {
326 for (
int k = 0; k < spaceDim; k++) {
327 int l = k + i * spaceDim + j * spaceDim * numFields;
328 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
330 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
333 *outStream <<
" At multi-index { ";
334 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
335 *outStream <<
"} computed value: " << vals(i,j,k)
336 <<
" but reference value: " << basisValues[l] <<
"\n";
344 vals.
resize(numFields, numPoints);
345 wedgeBasis.
getValues(vals, wedgeNodes, OPERATOR_DIV);
346 for (
int i = 0; i < numFields; i++) {
347 for (
int j = 0; j < numPoints; j++) {
348 int l = i + j * numFields;
349 if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
351 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
354 *outStream <<
" At multi-index { ";
355 *outStream << i <<
" ";*outStream << j <<
" ";
356 *outStream <<
"} computed divergence component: " << vals(i,j)
357 <<
" but reference divergence component: " << basisDivs[l] <<
"\n";
365 catch (std::logic_error err) {
366 *outStream << err.what() <<
"\n\n";
371 std::cout <<
"End Result: TEST FAILED\n";
373 std::cout <<
"End Result: TEST PASSED\n";
376 std::cout.copyfmt(oldFormatState);
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.
Implementation of the default H(div)-compatible FEM basis of degree 1 on Wedge cell.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Wedge cell.
Header file for the Intrepid::HDIV_WEDGE_I1_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...