Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/FunctionSpaceTools/test_03.cpp
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 
00049 #include "Intrepid_FunctionSpaceTools.hpp"
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Intrepid_CellTools.hpp"
00052 #include "Intrepid_HDIV_HEX_I1_FEM.hpp"
00053 #include "Intrepid_DefaultCubatureFactory.hpp"
00054 #include "Intrepid_Utils.hpp"
00055 #include "Teuchos_oblackholestream.hpp"
00056 #include "Teuchos_RCP.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using namespace std;
00060 using namespace Intrepid;
00061 
00062 int main(int argc, char *argv[]) {
00063 
00064   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00065 
00066   // This little trick lets us print to std::cout only if
00067   // a (dummy) command-line argument is provided.
00068   int iprint     = argc - 1;
00069   Teuchos::RCP<std::ostream> outStream;
00070   Teuchos::oblackholestream bhs; // outputs nothing
00071   if (iprint > 0)
00072     outStream = Teuchos::rcp(&std::cout, false);
00073   else
00074     outStream = Teuchos::rcp(&bhs, false);
00075 
00076   // Save the format state of the original std::cout.
00077   Teuchos::oblackholestream oldFormatState;
00078   oldFormatState.copyfmt(std::cout);
00079 
00080   *outStream \
00081     << "===============================================================================\n" \
00082     << "|                                                                             |\n" \
00083     << "|                      Unit Test (FunctionSpaceTools)                         |\n" \
00084     << "|                                                                             |\n" \
00085     << "|     1) Basic operator transformations and integration in HDIV:              |\n" \
00086     << "|                                                                             |\n" \
00087     << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00088     << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00089     << "|                                                                             |\n" \
00090     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00091     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00092     << "|                                                                             |\n" \
00093     << "===============================================================================\n";
00094   
00095 
00096   int errorFlag = 0;
00097 
00098   typedef FunctionSpaceTools fst; 
00099 
00100   *outStream \
00101   << "\n"
00102   << "===============================================================================\n"\
00103   << "| TEST 1: correctness of math operations                                      |\n"\
00104   << "===============================================================================\n";
00105 
00106   outStream->precision(20);
00107 
00108   try {
00109     
00110       shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >();    // cell type: hex
00111 
00112       /* Related to cubature. */
00113       DefaultCubatureFactory<double> cubFactory;                                                // create cubature factory
00114       int cubDegree = 20;                                                                       // cubature degree
00115       Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);           // create default cubature
00116       int spaceDim = myCub->getDimension();                                                     // get spatial dimension 
00117       int numCubPoints = myCub->getNumPoints();                                                 // get number of cubature points
00118 
00119       /* Related to basis. */
00120       Basis_HDIV_HEX_I1_FEM<double, FieldContainer<double> > hexBasis;                          // create H-div basis on a hex
00121       int numFields = hexBasis.getCardinality();                                                // get basis cardinality
00122  
00123       /* Cell geometries and orientations. */
00124       int numCells    = 4;
00125       int numNodes    = 8;
00126       int numCellData = numCells*numNodes*spaceDim;
00127       int numSignData = numCells*numFields;
00128 
00129       double hexnodes[] = {
00130         // hex 0  -- affine
00131         -1.0, -1.0, -1.0,
00132         1.0, -1.0, -1.0,
00133         1.0, 1.0, -1.0,
00134         -1.0, 1.0, -1.0,
00135         -1.0, -1.0, 1.0,
00136         1.0, -1.0, 1.0,
00137         1.0, 1.0, 1.0,
00138         -1.0, 1.0, 1.0,
00139         // hex 1  -- affine
00140         -3.0, -3.0, 1.0,
00141         6.0, 3.0, 1.0,
00142         7.0, 8.0, 0.0,
00143         -2.0, 2.0, 0.0,
00144         -3.0, -3.0, 4.0,
00145         6.0, 3.0, 4.0,
00146         7.0, 8.0, 3.0,
00147         -2.0, 2.0, 3.0,
00148         // hex 2  -- affine
00149         -3.0, -3.0, 0.0,
00150         9.0, 3.0, 0.0,
00151         15.0, 6.1, 0.0,
00152         3.0, 0.1, 0.0,
00153         9.0, 3.0, 0.1,
00154         21.0, 9.0, 0.1,
00155         27.0, 12.1, 0.1,
00156         15.0, 6.1, 0.1,
00157         // hex 3  -- nonaffine
00158         -2.0, -2.0, 0.0,
00159         2.0, -1.0, 0.0,
00160         1.0, 6.0, 0.0,
00161         -1.0, 1.0, 0.0,
00162         0.0, 0.0, 1.0,
00163         1.0, 0.0, 1.0,
00164         1.0, 1.0, 1.0,
00165         0.0, 1.0, 1.0
00166       };
00167 
00168       short facesigns[] = {
00169         1, 1, 1, 1, 1, 1,
00170         1, -1, 1, -1, 1, -1,
00171         -1, -1, 1, 1, -1, 1,
00172         -1, -1, 1, 1, -1, -1
00173       };
00174 
00175       /* Computational arrays. */
00176       FieldContainer<double> cub_points(numCubPoints, spaceDim);
00177       FieldContainer<double> cub_weights(numCubPoints);
00178       FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
00179       FieldContainer<short>  field_signs(numCells, numFields);
00180       FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
00181       //FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
00182       FieldContainer<double> jacobian_det(numCells, numCubPoints);
00183       FieldContainer<double> weighted_measure(numCells, numCubPoints);
00184 
00185       FieldContainer<double> div_of_basis_at_cub_points(numFields, numCubPoints);
00186       FieldContainer<double> transformed_div_of_basis_at_cub_points(numCells, numFields, numCubPoints);
00187       FieldContainer<double> weighted_transformed_div_of_basis_at_cub_points(numCells, numFields, numCubPoints);
00188       FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
00189 
00190       FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
00191       FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00192       FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00193       FieldContainer<double> mass_matrices(numCells, numFields, numFields);
00194 
00195       /******************* START COMPUTATION ***********************/
00196 
00197       // get cubature points and weights
00198       myCub->getCubature(cub_points, cub_weights);
00199 
00200       // fill cell vertex array
00201       cell_nodes.setValues(hexnodes, numCellData);
00202 
00203       // set basis function signs, for each cell
00204       field_signs.setValues(facesigns, numSignData);
00205 
00206       // compute geometric cell information
00207       CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
00208       //CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
00209       CellTools<double>::setJacobianDet(jacobian_det, jacobian);
00210 
00211       // compute weighted measure
00212       fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
00213 
00214       // Computing stiffness matrices:
00215       // tabulate divergences of basis functions at (reference) cubature points
00216       hexBasis.getValues(div_of_basis_at_cub_points, cub_points, OPERATOR_DIV);
00217 
00218       // transform divergences of basis functions
00219       fst::HDIVtransformDIV<double>(transformed_div_of_basis_at_cub_points,
00220                                     jacobian_det,
00221                                     div_of_basis_at_cub_points);
00222 
00223       // multiply with weighted measure
00224       fst::multiplyMeasure<double>(weighted_transformed_div_of_basis_at_cub_points,
00225                                    weighted_measure,
00226                                    transformed_div_of_basis_at_cub_points);
00227 
00228       // we can apply the field signs to the basis function arrays, or after the fact, see below
00229       fst::applyFieldSigns<double>(transformed_div_of_basis_at_cub_points, field_signs);
00230       fst::applyFieldSigns<double>(weighted_transformed_div_of_basis_at_cub_points, field_signs);
00231 
00232       // compute stiffness matrices
00233       fst::integrate<double>(stiffness_matrices,
00234                              transformed_div_of_basis_at_cub_points,
00235                              weighted_transformed_div_of_basis_at_cub_points,
00236                              COMP_CPP);
00237 
00238       // Computing mass matrices:
00239       // tabulate values of basis functions at (reference) cubature points
00240       hexBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
00241 
00242       // transform values of basis functions
00243       fst::HDIVtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
00244                                       jacobian,
00245                                       jacobian_det,
00246                                       value_of_basis_at_cub_points);
00247 
00248       // multiply with weighted measure
00249       fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
00250                                    weighted_measure,
00251                                    transformed_value_of_basis_at_cub_points);
00252 
00253       // compute mass matrices
00254       fst::integrate<double>(mass_matrices,
00255                              transformed_value_of_basis_at_cub_points,
00256                              weighted_transformed_value_of_basis_at_cub_points,
00257                              COMP_CPP);
00258 
00259       // apply field signs
00260       fst::applyLeftFieldSigns<double>(mass_matrices, field_signs);
00261       fst::applyRightFieldSigns<double>(mass_matrices, field_signs);
00262 
00263       /*******************  STOP COMPUTATION ***********************/
00264 
00265 
00266       /******************* START COMPARISON ***********************/
00267       string basedir = "./testdata";
00268       for (int cell_id = 0; cell_id < numCells-1; cell_id++) {
00269 
00270         stringstream namestream;
00271         string filename;
00272         namestream <<  basedir << "/mass_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00273         namestream >> filename;
00274 
00275         ifstream massfile(&filename[0]);
00276         if (massfile.is_open()) {
00277           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
00278             errorFlag++;
00279           massfile.close();
00280         }
00281         else {
00282           errorFlag = -1;
00283           std::cout << "End Result: TEST FAILED\n";
00284           return errorFlag;
00285         }
00286 
00287         namestream.clear();
00288         namestream << basedir << "/stiff_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00289         namestream >> filename;
00290 
00291         ifstream stifffile(&filename[0]);
00292         if (stifffile.is_open())
00293         {
00294           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
00295             errorFlag++;
00296           stifffile.close();
00297         }
00298         else {
00299           errorFlag = -1;
00300           std::cout << "End Result: TEST FAILED\n";
00301           return errorFlag;
00302         }
00303 
00304       }
00305       for (int cell_id = 3; cell_id < numCells; cell_id++) {
00306 
00307         stringstream namestream;
00308         string filename;
00309         namestream <<  basedir << "/mass_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00310         namestream >> filename;
00311 
00312         ifstream massfile(&filename[0]);
00313         if (massfile.is_open()) {
00314           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
00315             errorFlag++;
00316           massfile.close();
00317         }
00318         else {
00319           errorFlag = -1;
00320           std::cout << "End Result: TEST FAILED\n";
00321           return errorFlag;
00322         }
00323 
00324         namestream.clear();
00325         namestream << basedir << "/stiff_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00326         namestream >> filename;
00327 
00328         ifstream stifffile(&filename[0]);
00329         if (stifffile.is_open())
00330         {
00331           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
00332             errorFlag++;
00333           stifffile.close();
00334         }
00335         else {
00336           errorFlag = -1;
00337           std::cout << "End Result: TEST FAILED\n";
00338           return errorFlag;
00339         }
00340 
00341       }
00342       /******************* STOP COMPARISON ***********************/
00343 
00344       *outStream << "\n";
00345   }// try Basis_HDIV_HEX_I1
00346   catch (std::logic_error err) {
00347     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00348     *outStream << err.what() << '\n';
00349     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00350     errorFlag = -1000;
00351   };
00352 
00353 
00354   if (errorFlag != 0)
00355     std::cout << "End Result: TEST FAILED\n";
00356   else
00357     std::cout << "End Result: TEST PASSED\n";
00358 
00359   // reset format state of std::cout
00360   std::cout.copyfmt(oldFormatState);
00361 
00362   return errorFlag;
00363 }