00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00035 #include "Intrepid_FieldContainer.hpp"
00036 #include "Intrepid_HCURL_WEDGE_I1_FEM.hpp"
00037 #include "Teuchos_oblackholestream.hpp"
00038 #include "Teuchos_RCP.hpp"
00039 #include "Teuchos_GlobalMPISession.hpp"
00040
00041 using namespace std;
00042 using namespace Intrepid;
00043
00044 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
00045 { \
00046 ++nException; \
00047 try { \
00048 S ; \
00049 } \
00050 catch (std::logic_error err) { \
00051 ++throwCounter; \
00052 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00053 *outStream << err.what() << '\n'; \
00054 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
00055 }; \
00056 }
00057
00058 int main(int argc, char *argv[]) {
00059
00060 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00061
00062
00063
00064 int iprint = argc - 1;
00065 Teuchos::RCP<std::ostream> outStream;
00066 Teuchos::oblackholestream bhs;
00067 if (iprint > 0)
00068 outStream = Teuchos::rcp(&std::cout, false);
00069 else
00070 outStream = Teuchos::rcp(&bhs, false);
00071
00072
00073 Teuchos::oblackholestream oldFormatState;
00074 oldFormatState.copyfmt(std::cout);
00075
00076 *outStream \
00077 << "===============================================================================\n" \
00078 << "| |\n" \
00079 << "| Unit Test (Basis_HCURL_WEDGE_I1_FEM) |\n" \
00080 << "| |\n" \
00081 << "| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
00082 << "| 2) Basis values for VALUE and CURL operators |\n" \
00083 << "| |\n" \
00084 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
00085 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
00086 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \
00087 << "| |\n" \
00088 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
00089 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
00090 << "| |\n" \
00091 << "===============================================================================\n"\
00092 << "| TEST 1: Basis creation, exception testing |\n"\
00093 << "===============================================================================\n";
00094
00095
00096 Basis_HCURL_WEDGE_I1_FEM<double, FieldContainer<double> > wedgeBasis;
00097 int errorFlag = 0;
00098
00099
00100 int nException = 0;
00101 int throwCounter = 0;
00102
00103
00104 FieldContainer<double> wedgeNodes(12, 3);
00105 wedgeNodes(0,0) = 0.0; wedgeNodes(0,1) = 0.0; wedgeNodes(0,2) = -1.0;
00106 wedgeNodes(1,0) = 1.0; wedgeNodes(1,1) = 0.0; wedgeNodes(1,2) = -1.0;
00107 wedgeNodes(2,0) = 0.0; wedgeNodes(2,1) = 1.0; wedgeNodes(2,2) = -1.0;
00108 wedgeNodes(3,0) = 0.0; wedgeNodes(3,1) = 0.0; wedgeNodes(3,2) = 1.0;
00109 wedgeNodes(4,0) = 1.0; wedgeNodes(4,1) = 0.0; wedgeNodes(4,2) = 1.0;
00110 wedgeNodes(5,0) = 0.0; wedgeNodes(5,1) = 1.0; wedgeNodes(5,2) = 1.0;
00111
00112 wedgeNodes(6,0) = 0.25; wedgeNodes(6,1) = 0.5; wedgeNodes(6,2) = -1.0;
00113 wedgeNodes(7,0) = 0.5; wedgeNodes(7,1) = 0.25; wedgeNodes(7,2) = 0.0;
00114 wedgeNodes(8,0) = 0.25; wedgeNodes(8,1) = 0.25; wedgeNodes(8,2) = 1.0;
00115 wedgeNodes(9,0) = 0.25; wedgeNodes(9,1) = 0.0; wedgeNodes(9,2) = 0.75;
00116 wedgeNodes(10,0)= 0.0; wedgeNodes(10,1)= 0.5; wedgeNodes(10,2)= -0.25;
00117 wedgeNodes(11,0)= 0.5; wedgeNodes(11,1)= 0.5; wedgeNodes(11,2)= 0.0;
00118
00119
00120
00121
00122 FieldContainer<double> vals;
00123
00124 try{
00125
00126
00127 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
00128 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00129
00130
00131
00132 vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0));
00133 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00134
00135
00136
00137
00138 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00139
00140 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00141
00142 INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,4,1), throwCounter, nException );
00143
00144 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(10), throwCounter, nException );
00145
00146 INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
00147
00148 #ifdef HAVE_INTREPID_DEBUG
00149
00150
00151 FieldContainer<double> badPoints1(4, 5, 3);
00152 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00153
00154
00155 FieldContainer<double> badPoints2(4, 2);
00156 INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00157
00158
00159 FieldContainer<double> badVals1(4, 3);
00160 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00161
00162
00163 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_CURL), throwCounter, nException );
00164
00165
00166 FieldContainer<double> badVals2(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0), 3);
00167 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00168
00169
00170 FieldContainer<double> badVals3(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1, 3);
00171 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00172
00173
00174 FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
00175 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00176
00177
00178 INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_CURL), throwCounter, nException );
00179 #endif
00180
00181 }
00182 catch (std::logic_error err) {
00183 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00184 *outStream << err.what() << '\n';
00185 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00186 errorFlag = -1000;
00187 };
00188
00189
00190
00191 if (throwCounter != nException) {
00192 errorFlag++;
00193 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00194 }
00195
00196 *outStream \
00197 << "\n"
00198 << "===============================================================================\n"\
00199 << "| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
00200 << "===============================================================================\n";
00201
00202 try{
00203 std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
00204
00205
00206 for (unsigned i = 0; i < allTags.size(); i++) {
00207 int bfOrd = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00208
00209 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
00210 if( !( (myTag[0] == allTags[i][0]) &&
00211 (myTag[1] == allTags[i][1]) &&
00212 (myTag[2] == allTags[i][2]) &&
00213 (myTag[3] == allTags[i][3]) ) ) {
00214 errorFlag++;
00215 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00216 *outStream << " getDofOrdinal( {"
00217 << allTags[i][0] << ", "
00218 << allTags[i][1] << ", "
00219 << allTags[i][2] << ", "
00220 << allTags[i][3] << "}) = " << bfOrd <<" but \n";
00221 *outStream << " getDofTag(" << bfOrd << ") = { "
00222 << myTag[0] << ", "
00223 << myTag[1] << ", "
00224 << myTag[2] << ", "
00225 << myTag[3] << "}\n";
00226 }
00227 }
00228
00229
00230 for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
00231 std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
00232 int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00233 if( bfOrd != myBfOrd) {
00234 errorFlag++;
00235 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00236 *outStream << " getDofTag(" << bfOrd << ") = { "
00237 << myTag[0] << ", "
00238 << myTag[1] << ", "
00239 << myTag[2] << ", "
00240 << myTag[3] << "} but getDofOrdinal({"
00241 << myTag[0] << ", "
00242 << myTag[1] << ", "
00243 << myTag[2] << ", "
00244 << myTag[3] << "} ) = " << myBfOrd << "\n";
00245 }
00246 }
00247 }
00248 catch (std::logic_error err){
00249 *outStream << err.what() << "\n\n";
00250 errorFlag = -1000;
00251 };
00252
00253 *outStream \
00254 << "\n"
00255 << "===============================================================================\n"\
00256 << "| TEST 3: correctness of basis function values |\n"\
00257 << "===============================================================================\n";
00258
00259 outStream -> precision(20);
00260
00261
00262 double basisValues[] = {
00263 1.00000, 0, 0, 0, 0, 0, 0, -1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00264 0, 0.500000, 0, 0, 0, 0, 0, 0, 1.00000, 1.00000, 0, 0, 1.00000, 0, 0, \
00265 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.500000, 0, 0, 0, 0, \
00266 0, 0, -1.00000, 0, 0, -1.00000, -1.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00267 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00268 1.00000, 0, 0, 0, 0, 0, 0, -1.00000, 0, 0, 0, 0.500000, 0, 0, 0, 0, \
00269 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00000, 1.00000, 0, 0, 1.00000, 0, \
00270 0, 0, 0, 0, 0, 0, 0, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00271 0, 0, 0, -1.00000, 0, 0, -1.00000, -1.00000, 0, 0, 0, 0, 0, 0, 0, 0, \
00272 0, 0.500000, 0.500000, 0.250000, 0, -0.500000, 0.250000, 0, \
00273 -0.500000, -0.750000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.125000, \
00274 0, 0, 0.125000, 0, 0, 0.250000, 0.375000, 0.250000, 0, -0.125000, \
00275 0.250000, 0, -0.125000, -0.250000, 0, 0.375000, 0.250000, 0, \
00276 -0.125000, 0.250000, 0, -0.125000, -0.250000, 0, 0, 0, 0.125000, 0, \
00277 0, 0.250000, 0, 0, 0.125000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.750000, \
00278 0.250000, 0, -0.250000, 0.250000, 0, -0.250000, -0.750000, 0, 0, 0, \
00279 0.250000, 0, 0, 0.125000, 0, 0, 0.125000, 0.125000, 0.0312500, 0, 0, \
00280 0.0312500, 0, 0, -0.0937500, 0, 0.875000, 0.218750, 0, 0, 0.218750, \
00281 0, 0, -0.656250, 0, 0, 0, 0.375000, 0, 0, 0.125000, 0, 0, 0, \
00282 0.312500, 0, 0, -0.312500, 0, 0, -0.312500, -0.625000, 0, 0.187500, \
00283 0, 0, -0.187500, 0, 0, -0.187500, -0.375000, 0, 0, 0, 0.250000, 0, 0, \
00284 0, 0, 0, 0.250000, 0.250000, 0.250000, 0, -0.250000, 0.250000, 0, \
00285 -0.250000, -0.250000, 0, 0.250000, 0.250000, 0, -0.250000, 0.250000, \
00286 0, -0.250000, -0.250000, 0, 0, 0, 0, 0, 0, 0.250000, 0, 0, 0.250000
00287 };
00288
00289
00290 double basisCurls[] = {
00291 0, -0.500000, 2.00000, 0, 0, 2.00000, -0.500000, 0, 2.00000, 0, \
00292 0.500000, 0, 0, 0, 0, 0.500000, 0, 0, -0.500000, 0.500000, 0, 0, \
00293 -0.500000, 0, 0.500000, 0, 0, 0.500000, -0.500000, 2.00000, 0.500000, \
00294 0, 2.00000, 0, 0, 2.00000, -0.500000, 0.500000, 0, -0.500000, 0, 0, \
00295 0, 0, 0, -0.500000, 0.500000, 0, 0, -0.500000, 0, 0.500000, 0, 0, 0, \
00296 0, 2.00000, 0, 0.500000, 2.00000, -0.500000, 0.500000, 2.00000, 0, 0, \
00297 0, 0, -0.500000, 0, 0.500000, -0.500000, 0, -0.500000, 0.500000, 0, \
00298 0, -0.500000, 0, 0.500000, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, \
00299 0, 0, 0, 0.500000, 2.00000, 0, 0, 2.00000, 0.500000, 0, 2.00000, \
00300 -0.500000, 0.500000, 0, 0, -0.500000, 0, 0.500000, 0, 0, 0.500000, \
00301 -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, -0.500000, 0.500000, 2.00000, \
00302 -0.500000, 0, 2.00000, 0, 0, 2.00000, -0.500000, 0.500000, 0, 0, \
00303 -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, \
00304 0.500000, 0, 0, 0, 2.00000, 0, -0.500000, 2.00000, 0.500000, \
00305 -0.500000, 2.00000, -0.500000, 0.500000, 0, 0, -0.500000, 0, \
00306 0.500000, 0, 0, 0.125000, -0.250000, 2.00000, 0.125000, 0.250000, \
00307 2.00000, -0.375000, 0.250000, 2.00000, -0.125000, 0.250000, 0, \
00308 -0.125000, -0.250000, 0, 0.375000, -0.250000, 0, -0.500000, 0.500000, \
00309 0, 0, -0.500000, 0, 0.500000, 0, 0, 0.250000, -0.375000, 1.00000, \
00310 0.250000, 0.125000, 1.00000, -0.250000, 0.125000, 1.00000, -0.250000, \
00311 0.375000, 1.00000, -0.250000, -0.125000, 1.00000, 0.250000, \
00312 -0.125000, 1.00000, -0.500000, 0.500000, 0, 0, -0.500000, 0, \
00313 0.500000, 0, 0, 0.125000, -0.375000, 0, 0.125000, 0.125000, 0, \
00314 -0.375000, 0.125000, 0, -0.125000, 0.375000, 2.00000, -0.125000, \
00315 -0.125000, 2.00000, 0.375000, -0.125000, 2.00000, -0.500000, \
00316 0.500000, 0, 0, -0.500000, 0, 0.500000, 0, 0, 0.125000, -0.500000, \
00317 0.250000, 0.125000, 0, 0.250000, -0.375000, 0, 0.250000, -0.125000, \
00318 0.500000, 1.75000, -0.125000, 0, 1.75000, 0.375000, 0, 1.75000, \
00319 -0.500000, 0.500000, 0, 0, -0.500000, 0, 0.500000, 0, 0, 0, \
00320 -0.250000, 1.25000, 0, 0.250000, 1.25000, -0.500000, 0.250000, \
00321 1.25000, 0, 0.250000, 0.750000, 0, -0.250000, 0.750000, 0.500000, \
00322 -0.250000, 0.750000, -0.500000, 0.500000, 0, 0, -0.500000, 0, \
00323 0.500000, 0, 0, 0.250000, -0.250000, 1.00000, 0.250000, 0.250000, \
00324 1.00000, -0.250000, 0.250000, 1.00000, -0.250000, 0.250000, 1.00000, \
00325 -0.250000, -0.250000, 1.00000, 0.250000, -0.250000, 1.00000, \
00326 -0.500000, 0.500000, 0, 0, -0.500000, 0, 0.500000, 0, 0
00327 };
00328
00329 try{
00330
00331
00332 int numFields = wedgeBasis.getCardinality();
00333 int numPoints = wedgeNodes.dimension(0);
00334 int spaceDim = wedgeBasis.getBaseCellTopology().getDimension();
00335
00336
00337 FieldContainer<double> vals;
00338
00339
00340 vals.resize(numFields, numPoints, spaceDim);
00341 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
00342 for (int i = 0; i < numFields; i++) {
00343 for (int j = 0; j < numPoints; j++) {
00344 for (int k = 0; k < spaceDim; k++) {
00345
00346
00347 int l = k + i * spaceDim + j * spaceDim * numFields;
00348 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
00349 errorFlag++;
00350 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00351
00352
00353 *outStream << " At multi-index { ";
00354 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00355 *outStream << "} computed value: " << vals(i,j,k)
00356 << " but reference value: " << basisValues[l] << "\n";
00357 }
00358 }
00359 }
00360 }
00361
00362
00363 vals.resize(numFields, numPoints, spaceDim);
00364 wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_CURL);
00365 for (int i = 0; i < numFields; i++) {
00366 for (int j = 0; j < numPoints; j++) {
00367 for (int k = 0; k < spaceDim; k++) {
00368
00369
00370 int l = k + i * spaceDim + j * spaceDim * numFields;
00371 if (std::abs(vals(i,j,k) - basisCurls[l]) > INTREPID_TOL) {
00372 errorFlag++;
00373 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00374
00375
00376 *outStream << " At multi-index { ";
00377 *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00378 *outStream << "} computed curl component: " << vals(i,j,k)
00379 << " but reference curl component: " << basisCurls[l] << "\n";
00380 }
00381 }
00382 }
00383 }
00384
00385 }
00386
00387
00388 catch (std::logic_error err) {
00389 *outStream << err.what() << "\n\n";
00390 errorFlag = -1000;
00391 };
00392
00393 if (errorFlag != 0)
00394 std::cout << "End Result: TEST FAILED\n";
00395 else
00396 std::cout << "End Result: TEST PASSED\n";
00397
00398
00399 std::cout.copyfmt(oldFormatState);
00400
00401 return errorFlag;
00402 }