MOAB: Mesh Oriented datABase  (version 5.3.1)
MsqHessianTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
00028 // -*-
00029 //
00030 //   SUMMARY:
00031 //     USAGE:
00032 //
00033 //    AUTHOR: Thomas Leurent <tleurent@mcs.anl.gov>
00034 //       ORG: Argonne National Laboratory
00035 //    E-MAIL: tleurent@mcs.anl.gov
00036 //
00037 // ORIG-DATE: 13-Jan-03 at 09:05:56
00038 //  LAST-MOD:  9-Apr-03 at 17:17:26 by Thomas Leurent
00039 //
00040 // DESCRIPTION:
00041 // ============
00042 /*! \file MsqHessianTest.cpp
00043 
00044 Unit testing of the MsqHessian class.
00045 
00046 */
00047 // DESCRIP-END.
00048 //
00049 
00050 #include "PatchDataInstances.hpp"
00051 #include "MsqHessian.hpp"
00052 
00053 #include <cmath>
00054 
00055 #include "cppunit/extensions/HelperMacros.h"
00056 
00057 using namespace MBMesquite;
00058 using std::cerr;
00059 using std::cout;
00060 using std::endl;
00061 
00062 class MsqHessianTest : public CppUnit::TestFixture, public MBMesquite::MsqHessian
00063 {
00064 
00065   private:
00066     CPPUNIT_TEST_SUITE( MsqHessianTest );
00067     CPPUNIT_TEST( test_initialize );
00068     CPPUNIT_TEST( test_axpy );
00069     CPPUNIT_TEST( test_cg_solver );
00070     CPPUNIT_TEST( test_cholesky_preconditioner );
00071     CPPUNIT_TEST_SUITE_END();
00072 
00073   private:
00074     PatchData twoTriangles;
00075 
00076   public:
00077     void setUp()
00078     {
00079         MsqPrintError err( cout );
00080         create_two_tri_patch( twoTriangles, err );
00081         CPPUNIT_ASSERT( !err );
00082     }
00083 
00084     void tearDown()
00085     {
00086         destroy_patch_with_domain( twoTriangles );
00087     }
00088 
00089   public:
00090     MsqHessianTest() {}
00091 
00092     void accumulate_entries( const PatchData& pd, size_t elem_index, const Matrix3D* mat3d_array, MsqError& err );
00093 
00094     void test_initialize()
00095     {
00096         MsqPrintError err( cout );
00097 
00098         MsqHessian::initialize( twoTriangles, err );
00099         CPPUNIT_ASSERT( !err );
00100 
00101         // Checks values of mRowStart are correct.
00102         size_t i, row_start[] = { 0, 4, 7, 8, 9 };
00103         for( i = 0; i < 5; ++i )
00104             CPPUNIT_ASSERT( mRowStart[i] == row_start[i] );
00105 
00106         // Checks values of mColIndex are correct.
00107         size_t col_index[] = { 0, 1, 2, 3, 1, 2, 3, 2, 3 };
00108         for( i = 0; i < 9; ++i )
00109             CPPUNIT_ASSERT( mColIndex[i] == col_index[i] );
00110     }
00111 
00112     void test_axpy()
00113     {
00114         size_t i;
00115         MsqPrintError err( cout );
00116 
00117         MsqHessian::initialize( twoTriangles, err );
00118         CPPUNIT_ASSERT( !err );
00119 
00120         size_t hs = MsqHessian::size();
00121 
00122         Vector3D* res = new Vector3D[hs];
00123         Vector3D* x   = new Vector3D[hs];
00124         Vector3D* y   = new Vector3D[hs];
00125         Vector3D* ans = new Vector3D[hs];
00126 
00127         Matrix3D blocks[6];  // 6 blocks correspond to a triangular element (n+1)n/2 .
00128 
00129         blocks[0] = "4 4 7   4 5 7   7 7 3 ";
00130         blocks[1] = "4 8 7   3 5 7   1 2 3 ";
00131         blocks[2] = "4 4 7   6 5 9   1 8 5 ";
00132         blocks[3] = "4 4 2   4 5 3   2 3 3 ";
00133         blocks[4] = "2 4 7   3 2 7   1 4 3 ";
00134         blocks[5] = "8 4 9   4 5 7   9 7 3 ";
00135 
00136         accumulate_entries( twoTriangles, 0, blocks, err );
00137         CPPUNIT_ASSERT( !err );
00138 
00139         blocks[2] += blocks[5];
00140         blocks[5] = blocks[3];
00141 
00142         accumulate_entries( twoTriangles, 1, blocks, err );
00143         CPPUNIT_ASSERT( !err );
00144 
00145         Matrix3D entries_6_ans( "2 3 1   4 2 4   7 7 3" );
00146         CPPUNIT_ASSERT( mEntries[6] == entries_6_ans );
00147 
00148         x[0].set( 4, 5, 6 );
00149         x[1].set( 2, 5, 9 );
00150         x[2].set( 1, 2, 6 );
00151         x[3].set( 1, 5, 9 );
00152 
00153         y[0].set( 0, 0, 0 );
00154         y[1].set( 0, 0, 0 );
00155         y[2].set( 0, 0, 0 );
00156         y[3].set( 0, 0, 0 );
00157 
00158         axpy( res, hs, *this, x, hs, y, hs, err );
00159         CPPUNIT_ASSERT( !err );
00160 
00161         ans[0].set( 636, 635, 453 );
00162         ans[1].set( 365, 460, 461 );
00163         ans[2].set( 150, 199, 220 );
00164         ans[3].set( 166, 204, 174 );
00165 
00166         for( i = 0; i < hs; ++i )
00167             CPPUNIT_ASSERT( res[i] == ans[i] );
00168 
00169         y[0].set( 3, 2, 6 );
00170         y[1].set( 1, 2, 4 );
00171         y[2].set( 3, 6, 9 );
00172         y[3].set( 2, 4, 4 );
00173 
00174         ans[0].set( 639, 637, 459 );
00175         ans[1].set( 366, 462, 465 );
00176         ans[2].set( 153, 205, 229 );
00177         ans[3].set( 168, 208, 178 );
00178 
00179         axpy( res, hs, *this, x, hs, y, hs, err );
00180         CPPUNIT_ASSERT( !err );
00181 
00182         for( i = 0; i < hs; ++i )
00183             CPPUNIT_ASSERT( res[i] == ans[i] );
00184 
00185         delete[] res;
00186         delete[] x;
00187         delete[] y;
00188         delete[] ans;
00189     }
00190 
00191     void test_cg_solver()
00192     {
00193         size_t i;
00194         MsqPrintError err( cout );
00195 
00196         MsqHessian::initialize( twoTriangles, err );
00197         CPPUNIT_ASSERT( !err );
00198 
00199         size_t hs = MsqHessian::size();
00200 
00201         Vector3D* res = new Vector3D[hs];
00202         Vector3D* x   = new Vector3D[hs];
00203         Vector3D* y   = new Vector3D[hs];
00204         Vector3D* ans = new Vector3D[hs];
00205 
00206         Matrix3D blocks[6];  // 6 blocks correspond to a triangular element (n+!)n/2 .
00207 
00208         blocks[0] = "14 4 7   4 15 7   7 7 13 ";
00209         blocks[1] = "4 8 7   3 5 7   1 2 3 ";
00210         blocks[2] = "4 4 7   6 5 9   1 8 5 ";
00211         blocks[3] = "14 4 2   4 15 3   2 3 13 ";
00212         blocks[4] = "2 4 7   3 2 7   1 4 3 ";
00213         blocks[5] = "18 4 9   4 15 7   9 7 13 ";
00214 
00215         accumulate_entries( twoTriangles, 0, blocks, err );
00216         CPPUNIT_ASSERT( !err );
00217 
00218         blocks[2] -= blocks[5];
00219         blocks[5] = blocks[3];
00220 
00221         accumulate_entries( twoTriangles, 1, blocks, err );
00222         CPPUNIT_ASSERT( !err );
00223 
00224         y[0].set( 3, 2, 6 );
00225         y[1].set( 1, 2, 4 );
00226         y[2].set( 3, 6, 9 );
00227         y[3].set( 2, 4, 4 );
00228 
00229         cg_solver( x, y, err );
00230         CPPUNIT_ASSERT( !err );
00231 
00232         //     for (int i=0; i<4; ++i)
00233         //       cout << x[i];
00234 
00235         ans[0].set( 3.2338, 7.6431, -7.0735 );
00236         ans[1].set( -1.0068, 4.5520, -4.6628 );
00237         ans[2].set( -0.4361, 4.7640, -8.1006 );
00238         ans[3].set( -0.1218, -0.8817, -4.5571 );
00239 
00240         for( i = 0; i < hs; ++i )
00241             for( short j = 0; j < 3; ++j )
00242                 CPPUNIT_ASSERT_DOUBLES_EQUAL( x[i][j], ans[i][j], 10e-1 );
00243 
00244         delete[] res;
00245         delete[] x;
00246         delete[] y;
00247         delete[] ans;
00248     }
00249 
00250     void test_cholesky_preconditioner()
00251     {
00252         MsqPrintError err( cout );
00253 
00254         MsqHessian::initialize( twoTriangles, err );
00255         CPPUNIT_ASSERT( !err );
00256         MsqHessian::zero_out();
00257 
00258         Matrix3D blocks[6];  // 6 blocks correspond to a triangular element (n+1)n/2 .
00259 
00260         blocks[0] = " 2  1    1 "
00261                     " 1  2.5  0.5 "
00262                     " 1  0.5  2.5 ";
00263         blocks[1] = "0 0 0   0 0 0   0 0 0 ";
00264         blocks[2] = "0 0 0   0 0 0   0 0 0 ";
00265         blocks[3] = "0 0 0   0 0 0   0 0 0 ";
00266         blocks[4] = "0 0 0   0 0 0   0 0 0 ";
00267         blocks[5] = "0 0 0   0 0 0   0 0 0 ";
00268 
00269         accumulate_entries( twoTriangles, 0, blocks, err );
00270         CPPUNIT_ASSERT( !err );
00271 
00272         MsqHessian::compute_preconditioner( err );
00273         CPPUNIT_ASSERT( !err );
00274         Matrix3D block_0 = mPreconditioner[0];
00275 
00276         Matrix3D correct( " 0.5  0.5  0.5 "
00277                           " 0    0.5  0  "
00278                           " 0    0    0.5 " );
00279 
00280         for( short i = 0; i < 3; ++i )
00281             for( short j = 0; j < 3; ++j )
00282                 CPPUNIT_ASSERT_DOUBLES_EQUAL( block_0[i][j], correct[i][j], 10e-10 );
00283 
00284         //     cout << "block 0: \n" << block_0 << endl;
00285         //     cout << "correct: \n" << correct << endl;
00286     }
00287 };
00288 
00289 void MsqHessianTest::accumulate_entries( const PatchData& pd, size_t elem_index, const Matrix3D* mat3d_array,
00290                                          MsqError& err )
00291 {
00292     const MsqMeshEntity& elem = pd.element_by_index( elem_index );
00293     const size_t nv           = elem.vertex_count();
00294     const size_t* v           = elem.get_vertex_index_array();
00295     for( size_t r = 0; r < nv; ++r )
00296     {
00297         for( size_t c = r; c < nv; ++c )
00298         {
00299             add( v[r], v[c], *mat3d_array, err );
00300             CPPUNIT_ASSERT( !err );
00301             ++mat3d_array;
00302         }
00303     }
00304 }
00305 
00306 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqHessianTest, "MsqHessianTest" );
00307 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqHessianTest, "Unit" );
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines