Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HCURL_TRI_In_FEM/test

00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or
00025 //                    Denis Ridzal (dridzal@sandia.gov) or
00026 //                    Robert Kirby (robert.c.kirby@ttu.edu)
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00031 
00037 #include "Intrepid_FieldContainer.hpp"
00038 #include "Teuchos_oblackholestream.hpp"
00039 #include "Teuchos_RCP.hpp"
00040 #include "Teuchos_GlobalMPISession.hpp"
00041 #include "Intrepid_PointTools.hpp"
00042 #include "Intrepid_HDIV_TRI_In_FEM.hpp"
00043 #include "Intrepid_HCURL_TRI_In_FEM.hpp"
00044 #include "Shards_CellTopology.hpp"
00045 
00046 #include <iostream>
00047 using namespace Intrepid;
00048 
00053 int main(int argc, char *argv[]) {
00054 
00055   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00056   
00057   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00058   int iprint     = argc - 1;
00059   
00060   Teuchos::RCP<std::ostream> outStream;
00061   Teuchos::oblackholestream bhs; // outputs nothing
00062   
00063   if (iprint > 0)
00064     outStream = Teuchos::rcp(&std::cout, false);
00065   else
00066     outStream = Teuchos::rcp(&bhs, false);
00067   
00068   // Save the format state of the original std::cout.
00069   Teuchos::oblackholestream oldFormatState;
00070   oldFormatState.copyfmt(std::cout);
00071   
00072   *outStream \
00073     << "===============================================================================\n" \
00074     << "|                                                                             |\n" \
00075     << "|                           Unit Test HCURL_TRI_In_FEM                        |\n" \
00076     << "|                                                                             |\n" \
00077     << "|     1) Tests triangular Nedelec-Thomas basis                                |\n" \
00078     << "|                                                                             |\n" \
00079     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00080     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00081     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00082     << "|                                                                             |\n" \
00083     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00084     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00085     << "|                                                                             |\n" \
00086     << "===============================================================================\n";
00087   
00088   int errorFlag  = 0;
00089 
00090   // test for basis values, compare against H(div) basis, for they should be related by rotation
00091   // in the lowest order case
00092   try {
00093     const int deg = 1;
00094     Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >  myBasisDIV( deg , POINTTYPE_EQUISPACED );
00095     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasisCURL( deg , POINTTYPE_EQUISPACED );
00096 
00097     // Get a lattice
00098     const int np_lattice = PointTools::getLatticeSize( myBasisDIV.getBaseCellTopology() , deg , 0 );
00099     FieldContainer<double> lattice( np_lattice , 2 );
00100 
00101     FieldContainer<double> myBasisValuesDIV( myBasisDIV.getCardinality() , np_lattice , 2 );
00102     FieldContainer<double> myBasisValuesCURL( myBasisDIV.getCardinality() , np_lattice , 2 );
00103     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00104                                                             myBasisDIV.getBaseCellTopology() , 
00105                                                             deg , 
00106                                                             0 , 
00107                                                             POINTTYPE_EQUISPACED );    
00108 
00109     myBasisDIV.getValues( myBasisValuesDIV , lattice , OPERATOR_VALUE );
00110     myBasisCURL.getValues( myBasisValuesCURL, lattice , OPERATOR_VALUE );
00111 
00112     for (int i=0;i<myBasisValuesDIV.dimension(0);i++) {
00113       for (int j=0;j<myBasisValuesDIV.dimension(1);j++) {
00114         if (std::abs( myBasisValuesDIV(i,j,1) + myBasisValuesCURL(i,j,0) ) > INTREPID_TOL ) {
00115           errorFlag++;
00116           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00117           // Output the multi-index of the value where the error is:
00118           *outStream << " At multi-index { ";
00119           *outStream << i << " " << j << " and component 0";
00120           *outStream << "}  computed value: " << myBasisValuesCURL(i,j,0)
00121                     << " but correct value: " << -myBasisValuesDIV(i,j,1) << "\n";
00122           *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,0) + myBasisValuesDIV(i,j,1) ) << "\n";
00123         }
00124         if (std::abs( myBasisValuesDIV(i,j,0) - myBasisValuesCURL(i,j,1) ) > INTREPID_TOL ) {
00125           errorFlag++;
00126           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00127           // Output the multi-index of the value where the error is:
00128           *outStream << " At multi-index { ";
00129           *outStream << i << " " << j << " and component 1";
00130           *outStream << "}  computed value: " << myBasisValuesCURL(i,j,1)
00131                     << " but correct value: " << myBasisValuesDIV(i,j,0) << "\n";
00132           *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,1) - myBasisValuesDIV(i,j,1) ) << "\n";
00133         }
00134       }
00135     }
00136   }
00137   catch (std::exception err) {
00138     *outStream << err.what() << "\n\n";
00139     errorFlag = -1000;
00140   }
00141 
00142   // next test: code-code comparison with FIAT 
00143   try {
00144     const int deg = 2;
00145     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00146     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00147     FieldContainer<double> lattice( np_lattice , 2 );
00148     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00149                                                             myBasis.getBaseCellTopology() , 
00150                                                             deg , 
00151                                                             0 , 
00152                                                             POINTTYPE_EQUISPACED );        
00153     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 );
00154 
00155     
00156     myBasis.getValues( myBasisValues, lattice , OPERATOR_VALUE );
00157 
00158     const double fiat_vals[] = {
00159 2.000000000000000e+00, 0.000000000000000e+00,
00160 5.000000000000000e-01, 2.500000000000000e-01,
00161 -1.000000000000000e+00, -1.000000000000000e+00,
00162 2.500000000000000e-01, 0.000000000000000e+00,
00163 -5.000000000000000e-01, -5.000000000000000e-01,
00164 0.000000000000000e+00, 0.000000000000000e+00,
00165 -1.000000000000000e+00, 0.000000000000000e+00,
00166 5.000000000000000e-01, 2.500000000000000e-01,
00167 2.000000000000000e+00, 2.000000000000000e+00,
00168 -5.000000000000000e-01, 0.000000000000000e+00,
00169 2.500000000000000e-01, 2.500000000000000e-01,
00170 0.000000000000000e+00, 0.000000000000000e+00,
00171 0.000000000000000e+00, 0.000000000000000e+00,
00172 0.000000000000000e+00, 2.500000000000000e-01,
00173 0.000000000000000e+00, 2.000000000000000e+00,
00174 5.000000000000000e-01, 0.000000000000000e+00,
00175 -2.500000000000000e-01, 2.500000000000000e-01,
00176 1.000000000000000e+00, 0.000000000000000e+00,
00177 0.000000000000000e+00, 0.000000000000000e+00,
00178 0.000000000000000e+00, -5.000000000000000e-01,
00179 0.000000000000000e+00, -1.000000000000000e+00,
00180 -2.500000000000000e-01, 0.000000000000000e+00,
00181 -2.500000000000000e-01, 2.500000000000000e-01,
00182 -2.000000000000000e+00, 0.000000000000000e+00,
00183 0.000000000000000e+00, 1.000000000000000e+00,
00184 0.000000000000000e+00, 5.000000000000000e-01,
00185 0.000000000000000e+00, 0.000000000000000e+00,
00186 -2.500000000000000e-01, -5.000000000000000e-01,
00187 -2.500000000000000e-01, -2.500000000000000e-01,
00188 -2.000000000000000e+00, -2.000000000000000e+00,
00189 0.000000000000000e+00, -2.000000000000000e+00,
00190 0.000000000000000e+00, -2.500000000000000e-01,
00191 0.000000000000000e+00, 0.000000000000000e+00,
00192 -2.500000000000000e-01, -5.000000000000000e-01,
00193 5.000000000000000e-01, 5.000000000000000e-01,
00194 1.000000000000000e+00, 1.000000000000000e+00,
00195 0.000000000000000e+00, 0.000000000000000e+00,
00196 0.000000000000000e+00, -7.500000000000000e-01,
00197 0.000000000000000e+00, 0.000000000000000e+00,
00198 1.500000000000000e+00, 0.000000000000000e+00,
00199 7.500000000000000e-01, 7.500000000000000e-01,
00200 0.000000000000000e+00, 0.000000000000000e+00,
00201 0.000000000000000e+00, 0.000000000000000e+00,
00202 0.000000000000000e+00, 1.500000000000000e+00,
00203 0.000000000000000e+00, 0.000000000000000e+00,
00204 -7.500000000000000e-01, 0.000000000000000e+00,
00205 7.500000000000000e-01, 7.500000000000000e-01,
00206 0.000000000000000e+00, 0.000000000000000e+00
00207     };
00208 
00209     int cur=0;
00210     for (int i=0;i<myBasisValues.dimension(0);i++) {
00211       for (int j=0;j<myBasisValues.dimension(1);j++) {
00212         for (int k=0;k<myBasisValues.dimension(2);k++) {
00213           if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) {
00214             errorFlag++;
00215             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00216             
00217             // Output the multi-index of the value where the error is:
00218             *outStream << " At multi-index { ";
00219             *outStream << i << " " << j << " " << k;
00220             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00221                       << " but correct value: " << fiat_vals[cur] << "\n";
00222           *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
00223           }
00224           cur++;
00225         }
00226       }
00227     }
00228   }
00229   catch (std::exception err) {
00230     *outStream << err.what() << "\n\n";
00231     errorFlag = -1000;
00232   }
00233 
00234   try {
00235     const int deg = 2;
00236     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00237     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00238     FieldContainer<double> lattice( np_lattice , 2 );
00239     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00240                                                             myBasis.getBaseCellTopology() , 
00241                                                             deg , 
00242                                                             0 , 
00243                                                             POINTTYPE_EQUISPACED );        
00244     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice );
00245 
00246     
00247     myBasis.getValues( myBasisValues, lattice , OPERATOR_CURL );
00248 
00249     const double fiat_curls[] = {
00250 -7.000000000000000e+00,
00251 -2.500000000000000e+00,
00252 2.000000000000000e+00,
00253 -2.500000000000000e+00,
00254 2.000000000000000e+00,
00255 2.000000000000000e+00,
00256 2.000000000000000e+00,
00257 -2.500000000000000e+00,
00258 -7.000000000000000e+00,
00259 2.000000000000000e+00,
00260 -2.500000000000000e+00,
00261 2.000000000000000e+00,
00262 2.000000000000000e+00,
00263 -2.500000000000000e+00,
00264 -7.000000000000000e+00,
00265 2.000000000000000e+00,
00266 -2.500000000000000e+00,
00267 2.000000000000000e+00,
00268 2.000000000000000e+00,
00269 2.000000000000000e+00,
00270 2.000000000000000e+00,
00271 -2.500000000000000e+00,
00272 -2.500000000000000e+00,
00273 -7.000000000000000e+00,
00274 2.000000000000000e+00,
00275 2.000000000000000e+00,
00276 2.000000000000000e+00,
00277 -2.500000000000000e+00,
00278 -2.500000000000000e+00,
00279 -7.000000000000000e+00,
00280 -7.000000000000000e+00,
00281 -2.500000000000000e+00,
00282 2.000000000000000e+00,
00283 -2.500000000000000e+00,
00284 2.000000000000000e+00,
00285 2.000000000000000e+00,
00286 9.000000000000000e+00,
00287 4.500000000000000e+00,
00288 0.000000000000000e+00,
00289 0.000000000000000e+00,
00290 -4.500000000000000e+00,
00291 -9.000000000000000e+00,
00292 -9.000000000000000e+00,
00293 0.000000000000000e+00,
00294 9.000000000000000e+00,
00295 -4.500000000000000e+00,
00296 4.500000000000000e+00,
00297 0.000000000000000e+00
00298     };
00299 
00300     int cur=0;
00301     for (int i=0;i<myBasisValues.dimension(0);i++) {
00302       for (int j=0;j<myBasisValues.dimension(1);j++) {
00303         if (std::abs( myBasisValues(i,j) - fiat_curls[cur] ) > INTREPID_TOL ) {
00304           errorFlag++;
00305           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00306           
00307           // Output the multi-index of the value where the error is:
00308           *outStream << " At multi-index { ";
00309           *outStream << i << " " << j;
00310           *outStream << "}  computed value: " << myBasisValues(i,j)
00311                     << " but correct value: " << fiat_curls[cur] << "\n";
00312           *outStream << "Difference: " << std::abs( myBasisValues(i,j) - fiat_curls[cur] ) << "\n";
00313         }
00314         cur++;
00315       }
00316     }
00317   }
00318   catch (std::exception err) {
00319     *outStream << err.what() << "\n\n";
00320     errorFlag = -1000;
00321   }  
00322 
00323   if (errorFlag != 0)
00324     std::cout << "End Result: TEST FAILED\n";
00325   else
00326     std::cout << "End Result: TEST PASSED\n";
00327   
00328   // reset format state of std::cout
00329   std::cout.copyfmt(oldFormatState);
00330   
00331   return errorFlag;
00332 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HCURL_TRI_In_FEM/test_01.cpp
