50#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
55using namespace Intrepid;
57#define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
63 catch (const std::logic_error & err) { \
65 *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
66 *outStream << err.what() << '\n'; \
67 *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_HCURL_QUAD_I1_FEM) |\n" \
94 <<
"| 1) Conversion of Dof tags into Dof ordinals and back |\n" \
95 <<
"| 2) Basis values for VALUE and CURL 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";
109 Basis_HCURL_QUAD_I1_FEM<double, FieldContainer<double> > quadBasis;
114 int throwCounter = 0;
117 FieldContainer<double> quadNodes(9, 2);
118 quadNodes(0,0) = -1.0; quadNodes(0,1) = -1.0;
119 quadNodes(1,0) = 1.0; quadNodes(1,1) = -1.0;
120 quadNodes(2,0) = 1.0; quadNodes(2,1) = 1.0;
121 quadNodes(3,0) = -1.0; quadNodes(3,1) = 1.0;
123 quadNodes(4,0) = 0.0; quadNodes(4,1) = 0.0;
124 quadNodes(5,0) = 0.0; quadNodes(5,1) = -0.5;
125 quadNodes(6,0) = 0.0; quadNodes(6,1) = 0.5;
126 quadNodes(7,0) = -0.5; quadNodes(7,1) = 0.0;
127 quadNodes(8,0) = 0.5; quadNodes(8,1) = 0.0;
132 FieldContainer<double> vals;
136 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
137 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
141 vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
142 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
147 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
149 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
151 INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
153 INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
155 INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
157#ifdef HAVE_INTREPID_DEBUG
160 FieldContainer<double> badPoints1(4, 5, 3);
161 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
164 FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
165 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
168 FieldContainer<double> badVals1(4, 3);
169 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
171 FieldContainer<double> badCurls1(4,3,2);
173 INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
176 FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
177 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
180 FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
181 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
184 FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
185 INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
189 vals.resize(quadBasis.getCardinality(),
190 quadNodes.dimension(0),
191 Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
192 INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), throwCounter, nException );
196 catch (
const std::logic_error & err) {
197 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
198 *outStream << err.what() <<
'\n';
199 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
205 if (throwCounter != nException) {
207 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
213 <<
"===============================================================================\n"\
214 <<
"| TEST 2: correctness of tag to enum and enum to tag lookups |\n"\
215 <<
"===============================================================================\n";
218 std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
221 for (
unsigned i = 0; i < allTags.size(); i++) {
222 int bfOrd = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
224 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
225 if( !( (myTag[0] == allTags[i][0]) &&
226 (myTag[1] == allTags[i][1]) &&
227 (myTag[2] == allTags[i][2]) &&
228 (myTag[3] == allTags[i][3]) ) ) {
230 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
231 *outStream <<
" getDofOrdinal( {"
232 << allTags[i][0] <<
", "
233 << allTags[i][1] <<
", "
234 << allTags[i][2] <<
", "
235 << allTags[i][3] <<
"}) = " << bfOrd <<
" but \n";
236 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
240 << myTag[3] <<
"}\n";
245 for(
int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
246 std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
247 int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
248 if( bfOrd != myBfOrd) {
250 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
251 *outStream <<
" getDofTag(" << bfOrd <<
") = { "
255 << myTag[3] <<
"} but getDofOrdinal({"
259 << myTag[3] <<
"} ) = " << myBfOrd <<
"\n";
263 catch (
const std::logic_error & err){
264 *outStream << err.what() <<
"\n\n";
270 <<
"===============================================================================\n"\
271 <<
"| TEST 3: correctness of basis function values |\n"\
272 <<
"===============================================================================\n";
274 outStream -> precision(20);
277 double basisValues[] = {
278 0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \
279 0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \
280 -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \
281 0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \
282 0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \
283 0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \
284 -0.250000, 0, 0, -0.125000
288 double basisCurls[] = {
289 0.25, 0.25, 0.25, 0.25,
290 0.25, 0.25, 0.25, 0.25,
291 0.25, 0.25, 0.25, 0.25,
292 0.25, 0.25, 0.25, 0.25,
293 0.25, 0.25, 0.25, 0.25,
294 0.25, 0.25, 0.25, 0.25,
295 0.25, 0.25, 0.25, 0.25,
296 0.25, 0.25, 0.25, 0.25,
297 0.25, 0.25, 0.25, 0.25
304 int numFields = quadBasis.getCardinality();
305 int numPoints = quadNodes.dimension(0);
306 int spaceDim = quadBasis.getBaseCellTopology().getDimension();
309 FieldContainer<double> vals;
312 vals.resize(numFields, numPoints, spaceDim);
313 quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
314 for (
int i = 0; i < numFields; i++) {
315 for (
int j = 0; j < numPoints; j++) {
316 for (
int k = 0; k < spaceDim; k++) {
319 int l = k + i * spaceDim + j * spaceDim * numFields;
320 if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
322 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
325 *outStream <<
" At multi-index { ";
326 *outStream << i <<
" ";*outStream << j <<
" ";*outStream << k <<
" ";
327 *outStream <<
"} computed value: " << vals(i,j,k)
328 <<
" but reference value: " << basisValues[l] <<
"\n";
335 vals.resize(numFields, numPoints);
336 quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
337 for (
int i = 0; i < numFields; i++) {
338 for (
int j = 0; j < numPoints; j++) {
339 int l = i + j * numFields;
340 if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
342 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
345 *outStream <<
" At multi-index { ";
346 *outStream << i <<
" ";*outStream << j <<
" ";
347 *outStream <<
"} computed curl component: " << vals(i,j)
348 <<
" but reference curl component: " << basisCurls[l] <<
"\n";
355 catch (
const std::logic_error & err) {
356 *outStream << err.what() <<
"\n\n";
362 <<
"===============================================================================\n"\
363 <<
"| TEST 4: correctness of DoF locations |\n"\
364 <<
"===============================================================================\n";
367 Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
369 Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
370 Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
373 FieldContainer<double> cvals;
374 FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality(),2);
377#ifdef HAVE_INTREPID_DEBUG
379 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
381 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
383 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
385 cvals.resize(4,spaceDim);
386 INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
388 if (throwCounter != nException) {
390 *outStream << std::setw(70) <<
"^^^^----FAILURE!" <<
"\n";
394 FieldContainer<double> tangents(basis->getCardinality(),spaceDim);
395 tangents(0,0) = 2.0; tangents(0,1) = 0.0;
396 tangents(1,0) = 0.0; tangents(1,1) = 2.0;
397 tangents(2,0) = -2.0; tangents(2,1) = 0.0;
398 tangents(3,0) = 0.0; tangents(3,1) = -2.0;
400 basis->getValues(bvals, cvals, OPERATOR_VALUE);
402 for (
int i=0; i<bvals.dimension(0); i++) {
403 for (
int j=0; j<bvals.dimension(1); j++) {
405 double tangent = 0.0;
406 for(
int d=0;d<spaceDim;d++)
407 tangent += bvals(i,j,d)*tangents(j,d);
409 if ((i != j) && (std::abs(tangent - 0.0) > INTREPID_TOL)) {
411 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), tangent, 0.0);
412 *outStream << buffer;
414 else if ((i == j) && (std::abs(tangent - 1.0) > INTREPID_TOL)) {
416 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), tangent, 1.0);
417 *outStream << buffer;
423 catch (
const std::logic_error & err){
424 *outStream << err.what() <<
"\n\n";
430 std::cout <<
"End Result: TEST FAILED\n";
432 std::cout <<
"End Result: TEST PASSED\n";
435 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HCURL_QUAD_I1_FEM class.
Implementation of the default H(curl)-compatible FEM basis of degree 1 on Quadrilateral cell.