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