Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HCURL_TRI_In_FEM/test_01.cpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00044 
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Teuchos_oblackholestream.hpp"
00052 #include "Teuchos_RCP.hpp"
00053 #include "Teuchos_GlobalMPISession.hpp"
00054 #include "Intrepid_PointTools.hpp"
00055 #include "Intrepid_HDIV_TRI_In_FEM.hpp"
00056 #include "Intrepid_HCURL_TRI_In_FEM.hpp"
00057 #include "Shards_CellTopology.hpp"
00058 
00059 #include <iostream>
00060 using namespace Intrepid;
00061 
00066 int main(int argc, char *argv[]) {
00067 
00068   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00069   
00070   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
00071   int iprint     = argc - 1;
00072   
00073   Teuchos::RCP<std::ostream> outStream;
00074   Teuchos::oblackholestream bhs; // outputs nothing
00075   
00076   if (iprint > 0)
00077     outStream = Teuchos::rcp(&std::cout, false);
00078   else
00079     outStream = Teuchos::rcp(&bhs, false);
00080   
00081   // Save the format state of the original std::cout.
00082   Teuchos::oblackholestream oldFormatState;
00083   oldFormatState.copyfmt(std::cout);
00084   
00085   *outStream \
00086     << "===============================================================================\n" \
00087     << "|                                                                             |\n" \
00088     << "|                           Unit Test HCURL_TRI_In_FEM                        |\n" \
00089     << "|                                                                             |\n" \
00090     << "|     1) Tests triangular Nedelec-Thomas basis                                |\n" \
00091     << "|                                                                             |\n" \
00092     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00093     << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00094     << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00095     << "|                                                                             |\n" \
00096     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00097     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00098     << "|                                                                             |\n" \
00099     << "===============================================================================\n";
00100   
00101   int errorFlag  = 0;
00102 
00103   // test for basis values, compare against H(div) basis, for they should be related by rotation
00104   // in the lowest order case
00105   try {
00106     const int deg = 1;
00107     Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >  myBasisDIV( deg , POINTTYPE_EQUISPACED );
00108     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasisCURL( deg , POINTTYPE_EQUISPACED );
00109 
00110     // Get a lattice
00111     const int np_lattice = PointTools::getLatticeSize( myBasisDIV.getBaseCellTopology() , deg , 0 );
00112     FieldContainer<double> lattice( np_lattice , 2 );
00113 
00114     FieldContainer<double> myBasisValuesDIV( myBasisDIV.getCardinality() , np_lattice , 2 );
00115     FieldContainer<double> myBasisValuesCURL( myBasisDIV.getCardinality() , np_lattice , 2 );
00116     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00117                                                             myBasisDIV.getBaseCellTopology() , 
00118                                                             deg , 
00119                                                             0 , 
00120                                                             POINTTYPE_EQUISPACED );    
00121 
00122     myBasisDIV.getValues( myBasisValuesDIV , lattice , OPERATOR_VALUE );
00123     myBasisCURL.getValues( myBasisValuesCURL, lattice , OPERATOR_VALUE );
00124 
00125     for (int i=0;i<myBasisValuesDIV.dimension(0);i++) {
00126       for (int j=0;j<myBasisValuesDIV.dimension(1);j++) {
00127         if (std::abs( myBasisValuesDIV(i,j,1) + myBasisValuesCURL(i,j,0) ) > INTREPID_TOL ) {
00128           errorFlag++;
00129           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00130           // Output the multi-index of the value where the error is:
00131           *outStream << " At multi-index { ";
00132           *outStream << i << " " << j << " and component 0";
00133           *outStream << "}  computed value: " << myBasisValuesCURL(i,j,0)
00134                     << " but correct value: " << -myBasisValuesDIV(i,j,1) << "\n";
00135           *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,0) + myBasisValuesDIV(i,j,1) ) << "\n";
00136         }
00137         if (std::abs( myBasisValuesDIV(i,j,0) - myBasisValuesCURL(i,j,1) ) > INTREPID_TOL ) {
00138           errorFlag++;
00139           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00140           // Output the multi-index of the value where the error is:
00141           *outStream << " At multi-index { ";
00142           *outStream << i << " " << j << " and component 1";
00143           *outStream << "}  computed value: " << myBasisValuesCURL(i,j,1)
00144                     << " but correct value: " << myBasisValuesDIV(i,j,0) << "\n";
00145           *outStream << "Difference: " << std::abs( myBasisValuesCURL(i,j,1) - myBasisValuesDIV(i,j,1) ) << "\n";
00146         }
00147       }
00148     }
00149   }
00150   catch (std::exception err) {
00151     *outStream << err.what() << "\n\n";
00152     errorFlag = -1000;
00153   }
00154 
00155   // next test: code-code comparison with FIAT 
00156   try {
00157     const int deg = 2;
00158     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00159     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00160     FieldContainer<double> lattice( np_lattice , 2 );
00161     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00162                                                             myBasis.getBaseCellTopology() , 
00163                                                             deg , 
00164                                                             0 , 
00165                                                             POINTTYPE_EQUISPACED );        
00166     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice , 2 );
00167 
00168     
00169     myBasis.getValues( myBasisValues, lattice , OPERATOR_VALUE );
00170 
00171     const double fiat_vals[] = {
00172 2.000000000000000e+00, 0.000000000000000e+00,
00173 5.000000000000000e-01, 2.500000000000000e-01,
00174 -1.000000000000000e+00, -1.000000000000000e+00,
00175 2.500000000000000e-01, 0.000000000000000e+00,
00176 -5.000000000000000e-01, -5.000000000000000e-01,
00177 0.000000000000000e+00, 0.000000000000000e+00,
00178 -1.000000000000000e+00, 0.000000000000000e+00,
00179 5.000000000000000e-01, 2.500000000000000e-01,
00180 2.000000000000000e+00, 2.000000000000000e+00,
00181 -5.000000000000000e-01, 0.000000000000000e+00,
00182 2.500000000000000e-01, 2.500000000000000e-01,
00183 0.000000000000000e+00, 0.000000000000000e+00,
00184 0.000000000000000e+00, 0.000000000000000e+00,
00185 0.000000000000000e+00, 2.500000000000000e-01,
00186 0.000000000000000e+00, 2.000000000000000e+00,
00187 5.000000000000000e-01, 0.000000000000000e+00,
00188 -2.500000000000000e-01, 2.500000000000000e-01,
00189 1.000000000000000e+00, 0.000000000000000e+00,
00190 0.000000000000000e+00, 0.000000000000000e+00,
00191 0.000000000000000e+00, -5.000000000000000e-01,
00192 0.000000000000000e+00, -1.000000000000000e+00,
00193 -2.500000000000000e-01, 0.000000000000000e+00,
00194 -2.500000000000000e-01, 2.500000000000000e-01,
00195 -2.000000000000000e+00, 0.000000000000000e+00,
00196 0.000000000000000e+00, 1.000000000000000e+00,
00197 0.000000000000000e+00, 5.000000000000000e-01,
00198 0.000000000000000e+00, 0.000000000000000e+00,
00199 -2.500000000000000e-01, -5.000000000000000e-01,
00200 -2.500000000000000e-01, -2.500000000000000e-01,
00201 -2.000000000000000e+00, -2.000000000000000e+00,
00202 0.000000000000000e+00, -2.000000000000000e+00,
00203 0.000000000000000e+00, -2.500000000000000e-01,
00204 0.000000000000000e+00, 0.000000000000000e+00,
00205 -2.500000000000000e-01, -5.000000000000000e-01,
00206 5.000000000000000e-01, 5.000000000000000e-01,
00207 1.000000000000000e+00, 1.000000000000000e+00,
00208 0.000000000000000e+00, 0.000000000000000e+00,
00209 0.000000000000000e+00, -7.500000000000000e-01,
00210 0.000000000000000e+00, 0.000000000000000e+00,
00211 1.500000000000000e+00, 0.000000000000000e+00,
00212 7.500000000000000e-01, 7.500000000000000e-01,
00213 0.000000000000000e+00, 0.000000000000000e+00,
00214 0.000000000000000e+00, 0.000000000000000e+00,
00215 0.000000000000000e+00, 1.500000000000000e+00,
00216 0.000000000000000e+00, 0.000000000000000e+00,
00217 -7.500000000000000e-01, 0.000000000000000e+00,
00218 7.500000000000000e-01, 7.500000000000000e-01,
00219 0.000000000000000e+00, 0.000000000000000e+00
00220     };
00221 
00222     int cur=0;
00223     for (int i=0;i<myBasisValues.dimension(0);i++) {
00224       for (int j=0;j<myBasisValues.dimension(1);j++) {
00225         for (int k=0;k<myBasisValues.dimension(2);k++) {
00226           if (std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) > INTREPID_TOL ) {
00227             errorFlag++;
00228             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00229             
00230             // Output the multi-index of the value where the error is:
00231             *outStream << " At multi-index { ";
00232             *outStream << i << " " << j << " " << k;
00233             *outStream << "}  computed value: " << myBasisValues(i,j,k)
00234                       << " but correct value: " << fiat_vals[cur] << "\n";
00235           *outStream << "Difference: " << std::abs( myBasisValues(i,j,k) - fiat_vals[cur] ) << "\n";
00236           }
00237           cur++;
00238         }
00239       }
00240     }
00241   }
00242   catch (std::exception err) {
00243     *outStream << err.what() << "\n\n";
00244     errorFlag = -1000;
00245   }
00246 
00247   try {
00248     const int deg = 2;
00249     Basis_HCURL_TRI_In_FEM<double,FieldContainer<double> > myBasis( deg , POINTTYPE_EQUISPACED );
00250     const int np_lattice = PointTools::getLatticeSize( myBasis.getBaseCellTopology() , deg , 0 );
00251     FieldContainer<double> lattice( np_lattice , 2 );
00252     PointTools::getLattice<double,FieldContainer<double> >( lattice , 
00253                                                             myBasis.getBaseCellTopology() , 
00254                                                             deg , 
00255                                                             0 , 
00256                                                             POINTTYPE_EQUISPACED );        
00257     FieldContainer<double> myBasisValues( myBasis.getCardinality() , np_lattice );
00258 
00259     
00260     myBasis.getValues( myBasisValues, lattice , OPERATOR_CURL );
00261 
00262     const double fiat_curls[] = {
00263 7.000000000000000e+00,
00264 2.500000000000000e+00,
00265 -2.000000000000000e+00,
00266 2.500000000000000e+00,
00267 -2.000000000000000e+00,
00268 -2.000000000000000e+00,
00269 -2.000000000000000e+00,
00270 2.500000000000000e+00,
00271 7.000000000000000e+00,
00272 -2.000000000000000e+00,
00273 2.500000000000000e+00,
00274 -2.000000000000000e+00,
00275 -2.000000000000000e+00,
00276 2.500000000000000e+00,
00277 7.000000000000000e+00,
00278 -2.000000000000000e+00,
00279 2.500000000000000e+00,
00280 -2.000000000000000e+00,
00281 -2.000000000000000e+00,
00282 -2.000000000000000e+00,
00283 -2.000000000000000e+00,
00284 2.500000000000000e+00,
00285 2.500000000000000e+00,
00286 7.000000000000000e+00,
00287 -2.000000000000000e+00,
00288 -2.000000000000000e+00,
00289 -2.000000000000000e+00,
00290 2.500000000000000e+00,
00291 2.500000000000000e+00,
00292 7.000000000000000e+00,
00293 7.000000000000000e+00,
00294 2.500000000000000e+00,
00295 -2.000000000000000e+00,
00296 2.500000000000000e+00,
00297 -2.000000000000000e+00,
00298 -2.000000000000000e+00,
00299 -9.000000000000000e+00,
00300 -4.500000000000000e+00,
00301 0.000000000000000e+00,
00302 0.000000000000000e+00,
00303 4.500000000000000e+00,
00304 9.000000000000000e+00,
00305 9.000000000000000e+00,
00306 0.000000000000000e+00,
00307 -9.000000000000000e+00,
00308 4.500000000000000e+00,
00309 -4.500000000000000e+00,
00310 0.000000000000000e+00
00311     };
00312 
00313     int cur=0;
00314     for (int i=0;i<myBasisValues.dimension(0);i++) {
00315       for (int j=0;j<myBasisValues.dimension(1);j++) {
00316         if (std::abs( myBasisValues(i,j) - fiat_curls[cur] ) > INTREPID_TOL ) {
00317           errorFlag++;
00318           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00319           
00320           // Output the multi-index of the value where the error is:
00321           *outStream << " At multi-index { ";
00322           *outStream << i << " " << j;
00323           *outStream << "}  computed value: " << myBasisValues(i,j)
00324                     << " but correct value: " << fiat_curls[cur] << "\n";
00325           *outStream << "Difference: " << std::abs( myBasisValues(i,j) - fiat_curls[cur] ) << "\n";
00326         }
00327         cur++;
00328       }
00329     }
00330   }
00331   catch (std::exception err) {
00332     *outStream << err.what() << "\n\n";
00333     errorFlag = -1000;
00334   }  
00335 
00336   if (errorFlag != 0)
00337     std::cout << "End Result: TEST FAILED\n";
00338   else
00339     std::cout << "End Result: TEST PASSED\n";
00340   
00341   // reset format state of std::cout
00342   std::cout.copyfmt(oldFormatState);
00343   
00344   return errorFlag;
00345 }