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