Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/FunctionSpaceTools/test_02.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_HCURL_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 HCURL              |\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       shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >();    // cell type: hex
00110 
00111       /* Related to cubature. */
00112       DefaultCubatureFactory<double> cubFactory;                                                // create cubature factory
00113       int cubDegree = 20;                                                                       // cubature degree
00114       Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);           // create default cubature
00115       int spaceDim = myCub->getDimension();                                                     // get spatial dimension 
00116       int numCubPoints = myCub->getNumPoints();                                                 // get number of cubature points
00117 
00118       /* Related to basis. */
00119       Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexBasis;                         // create H-curl basis on a hex
00120       int numFields = hexBasis.getCardinality();                                                // get basis cardinality
00121  
00122       /* Cell geometries and orientations. */
00123       int numCells    = 4;
00124       int numNodes    = 8;
00125       int numCellData = numCells*numNodes*spaceDim;
00126       int numSignData = numCells*numFields;
00127 
00128       double hexnodes[] = {
00129         // hex 0  -- affine
00130         -1.0, -1.0, -1.0,
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         // hex 1  -- affine
00139         -3.0, -3.0, 1.0,
00140         6.0, 3.0, 1.0,
00141         7.0, 8.0, 0.0,
00142         -2.0, 2.0, 0.0,
00143         -3.0, -3.0, 4.0,
00144         6.0, 3.0, 4.0,
00145         7.0, 8.0, 3.0,
00146         -2.0, 2.0, 3.0,
00147         // hex 2  -- affine
00148         -3.0, -3.0, 0.0,
00149         9.0, 3.0, 0.0,
00150         15.0, 6.1, 0.0,
00151         3.0, 0.1, 0.0,
00152         9.0, 3.0, 0.1,
00153         21.0, 9.0, 0.1,
00154         27.0, 12.1, 0.1,
00155         15.0, 6.1, 0.1,
00156         // hex 3  -- nonaffine
00157         -2.0, -2.0, 0.0,
00158         2.0, -1.0, 0.0,
00159         1.0, 6.0, 0.0,
00160         -1.0, 1.0, 0.0,
00161         0.0, 0.0, 1.0,
00162         1.0, 0.0, 1.0,
00163         1.0, 1.0, 1.0,
00164         0.0, 1.0, 1.0
00165       };
00166 
00167       short edgesigns[] = {
00168         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00169         1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1,
00170         -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1,
00171         1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1
00172       };
00173 
00174       /* Computational arrays. */
00175       FieldContainer<double> cub_points(numCubPoints, spaceDim);
00176       FieldContainer<double> cub_weights(numCubPoints);
00177       FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
00178       FieldContainer<short>  field_signs(numCells, numFields);
00179       FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
00180       FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
00181       FieldContainer<double> jacobian_det(numCells, numCubPoints);
00182       FieldContainer<double> weighted_measure(numCells, numCubPoints);
00183 
00184       FieldContainer<double> curl_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
00185       FieldContainer<double> transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00186       FieldContainer<double> weighted_transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00187       FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
00188 
00189       FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
00190       FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00191       FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
00192       FieldContainer<double> mass_matrices(numCells, numFields, numFields);
00193 
00194       /******************* START COMPUTATION ***********************/
00195 
00196       // get cubature points and weights
00197       myCub->getCubature(cub_points, cub_weights);
00198 
00199       // fill cell vertex array
00200       cell_nodes.setValues(hexnodes, numCellData);
00201 
00202       // set basis function signs, for each cell
00203       field_signs.setValues(edgesigns, numSignData);
00204 
00205       // compute geometric cell information
00206       CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
00207       CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
00208       CellTools<double>::setJacobianDet(jacobian_det, jacobian);
00209 
00210       // compute weighted measure
00211       fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
00212 
00213       // Computing stiffness matrices:
00214       // tabulate curls of basis functions at (reference) cubature points
00215       hexBasis.getValues(curl_of_basis_at_cub_points, cub_points, OPERATOR_CURL);
00216 
00217       // transform curls of basis functions 
00218       fst::HCURLtransformCURL<double>(transformed_curl_of_basis_at_cub_points,
00219                                       jacobian,
00220                                       jacobian_det,
00221                                       curl_of_basis_at_cub_points);
00222 
00223       // multiply with weighted measure
00224       fst::multiplyMeasure<double>(weighted_transformed_curl_of_basis_at_cub_points,
00225                                    weighted_measure,
00226                                    transformed_curl_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_curl_of_basis_at_cub_points, field_signs);
00230       fst::applyFieldSigns<double>(weighted_transformed_curl_of_basis_at_cub_points, field_signs);
00231 
00232       // compute stiffness matrices
00233       fst::integrate<double>(stiffness_matrices,
00234                              transformed_curl_of_basis_at_cub_points,
00235                              weighted_transformed_curl_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::HCURLtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
00244                                        jacobian_inv,
00245                                        value_of_basis_at_cub_points);
00246 
00247       // multiply with weighted measure
00248       fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
00249                                    weighted_measure,
00250                                    transformed_value_of_basis_at_cub_points);
00251 
00252       // compute mass matrices
00253       fst::integrate<double>(mass_matrices,
00254                              transformed_value_of_basis_at_cub_points,
00255                              weighted_transformed_value_of_basis_at_cub_points,
00256                              COMP_CPP);
00257 
00258       // apply field signs (after the fact, as a post-processing step)
00259       fst::applyLeftFieldSigns<double>(mass_matrices, field_signs);
00260       fst::applyRightFieldSigns<double>(mass_matrices, field_signs);
00261 
00262       /*******************  STOP COMPUTATION ***********************/
00263 
00264 
00265       /******************* START COMPARISON ***********************/
00266       string basedir = "./testdata";
00267       for (int cell_id = 0; cell_id < numCells-1; cell_id++) {
00268 
00269         stringstream namestream;
00270         string filename;
00271         namestream <<  basedir << "/mass_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00272         namestream >> filename;
00273 
00274         ifstream massfile(&filename[0]);
00275         if (massfile.is_open()) {
00276           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
00277             errorFlag++;
00278           massfile.close();
00279         }
00280         else {
00281           errorFlag = -1;
00282           std::cout << "End Result: TEST FAILED\n";
00283           return errorFlag;
00284         }
00285 
00286         namestream.clear();
00287         namestream << basedir << "/stiff_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00288         namestream >> filename;
00289 
00290         ifstream stifffile(&filename[0]);
00291         if (stifffile.is_open())
00292         {
00293           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
00294             errorFlag++;
00295           stifffile.close();
00296         }
00297         else {
00298           errorFlag = -1;
00299           std::cout << "End Result: TEST FAILED\n";
00300           return errorFlag;
00301         }
00302 
00303       }
00304       for (int cell_id = 3; cell_id < numCells; cell_id++) {
00305 
00306         stringstream namestream;
00307         string filename;
00308         namestream <<  basedir << "/mass_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00309         namestream >> filename;
00310 
00311         ifstream massfile(&filename[0]);
00312         if (massfile.is_open()) {
00313           if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
00314             errorFlag++;
00315           massfile.close();
00316         }
00317         else {
00318           errorFlag = -1;
00319           std::cout << "End Result: TEST FAILED\n";
00320           return errorFlag;
00321         }
00322 
00323         namestream.clear();
00324         namestream << basedir << "/stiff_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
00325         namestream >> filename;
00326 
00327         ifstream stifffile(&filename[0]);
00328         if (stifffile.is_open())
00329         {
00330           if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
00331             errorFlag++;
00332           stifffile.close();
00333         }
00334         else {
00335           errorFlag = -1;
00336           std::cout << "End Result: TEST FAILED\n";
00337           return errorFlag;
00338         }
00339 
00340       }
00341       /******************* STOP COMPARISON ***********************/
00342 
00343       *outStream << "\n";
00344   }
00345   catch (std::logic_error err) {
00346     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00347     *outStream << err.what() << '\n';
00348     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00349     errorFlag = -1000;
00350   };
00351 
00352 
00353   if (errorFlag != 0)
00354     std::cout << "End Result: TEST FAILED\n";
00355   else
00356     std::cout << "End Result: TEST PASSED\n";
00357 
00358   // reset format state of std::cout
00359   std::cout.copyfmt(oldFormatState);
00360 
00361   return errorFlag;
00362 }