MOAB: Mesh Oriented datABase  (version 5.2.1)
NonSmoothDescent.hpp
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 /*!
00028   \file   NonSmoothDescent.hpp
00029   \brief
00030 
00031   The NonSmoothDescent Class implements the steepest descent algorythm in
00032   order to move a free vertex to an optimal position given an
00033   ObjectiveFunction object and a QaulityMetric object.
00034 
00035   \author Thomas Leurent
00036   \date   2002-06-13
00037 */
00038 
00039 #ifndef Mesquite_NonSmoothDescent_hpp
00040 #define Mesquite_NonSmoothDescent_hpp
00041 
00042 #include "Mesquite.hpp"
00043 #include "VertexMover.hpp"
00044 #include "ObjectiveFunction.hpp"
00045 #include "MsqFreeVertexIndexIterator.hpp"
00046 #include "MsqDebug.hpp"
00047 #include "QualityMetric.hpp"
00048 #include "VertexPatches.hpp"
00049 #include <vector>
00050 #include "Vector3D.hpp"
00051 
00052 namespace MBMesquite
00053 {
00054 
00055 class NonSmoothDescent : public VertexMover
00056 {
00057   public:
00058     MESQUITE_EXPORT
00059     NonSmoothDescent( QualityMetric* qm );
00060 
00061     MESQUITE_EXPORT virtual ~NonSmoothDescent() {}
00062 
00063     MESQUITE_EXPORT virtual std::string get_name() const;
00064 
00065     MESQUITE_EXPORT virtual PatchSet* get_patch_set();
00066 
00067   protected:
00068     struct ActiveSet
00069     {
00070         int num_equal;
00071         double true_active_value;
00072         std::vector< int > active_ind;
00073 
00074         void set( int index )
00075         {
00076             active_ind.clear();
00077             active_ind.push_back( index );
00078             num_equal = 0;
00079         }
00080         void add( int index, bool equal )
00081         {
00082             active_ind.push_back( index );
00083             num_equal += equal;
00084         }
00085     };
00086 
00087     class SymmetricMatrix
00088     {
00089         double* storage;
00090         size_t size;
00091         size_t index( size_t r, size_t c ) const
00092         {
00093             return ( r <= c ) ? size * r - r * ( r + 1 ) / 2 + c : size * c - c * ( c + 1 ) / 2 + r;
00094         }
00095 
00096       public:
00097         SymmetricMatrix() : storage( 0 ), size( 0 ) {}
00098 
00099         ~SymmetricMatrix()
00100         {
00101             free( storage );
00102         }
00103 
00104         void resize( size_t new_size )
00105         {
00106             size    = new_size;
00107             storage = (double*)realloc( storage, sizeof( double ) * size * ( size + 1 ) / 2 );
00108         }
00109 
00110         void fill( double value )
00111         {
00112             std::fill( storage, storage + ( size * ( size + 1 ) / 2 ), value );
00113         }
00114 
00115         double& operator()( size_t r, size_t c )
00116         {
00117             return storage[index( r, c )];
00118         }
00119 
00120         double operator()( size_t r, size_t c ) const
00121         {
00122             return storage[index( r, c )];
00123         }
00124 
00125         double condition3x3() const;
00126     };
00127 
00128     MESQUITE_EXPORT virtual void initialize( PatchData& pd, MsqError& err );
00129 
00130     MESQUITE_EXPORT virtual void optimize_vertex_positions( PatchData& pd, MsqError& err );
00131 
00132     MESQUITE_EXPORT virtual void initialize_mesh_iteration( PatchData& pd, MsqError& err );
00133 
00134     MESQUITE_EXPORT virtual void terminate_mesh_iteration( PatchData& pd, MsqError& err );
00135 
00136     MESQUITE_EXPORT virtual void cleanup();
00137 
00138   private:
00139     enum OptStatus
00140     {
00141         MSQ_NO_STATUS     = 0,
00142         MSQ_STEP_ACCEPTED = 100,
00143         MSQ_IMP_TOO_SMALL,
00144         MSQ_FLAT_NO_IMP,
00145         MSQ_STEP_TOO_SMALL,
00146         MSQ_EQUILIBRIUM,
00147         MSQ_ZERO_SEARCH,
00148         MSQ_MAX_ITER_EXCEEDED
00149     };
00150 
00151     enum Status
00152     {
00153         MSQ_NO_EQUIL = 101,
00154         MSQ_CHECK_TOP_DOWN,
00155         MSQ_CHECK_BOTTOM_UP,
00156         MSQ_TWO_PT_PLANE,
00157         MSQ_THREE_PT_PLANE,
00158         MSQ_CHECK_Y_COORD_DIRECTION,
00159         MSQ_CHECK_X_COORD_DIRECTION,
00160         MSQ_CHECK_Z_COORD_DIRECTION,
00161         MSQ_EQUIL,
00162         MSQ_HULL_TEST_ERROR
00163     };
00164 
00165     enum Direction
00166     {
00167         MSQ_XDIR = 0,
00168         MSQ_YDIR = 1,
00169         MSQ_ZDIR = 2
00170     };
00171 
00172     /* local copy of patch data */
00173     //    PatchData patch_data;
00174     size_t freeVertexIndex;
00175 
00176     /* smoothing parameters */
00177     double minStepSize;
00178 
00179     /* optimization data */
00180     VertexPatches patchSet;
00181     QualityMetric* currentQM;
00182     double originalValue;
00183     std::vector< double > mFunction, testFunction, originalFunction;
00184     std::vector< double > mGS;
00185     SymmetricMatrix mG;
00186     Matrix3D mPDG;
00187     std::vector< double > prevActiveValues;
00188     std::vector< Vector3D > mGradient, tmpGrad;
00189     std::vector< size_t > tmpIdx, qmHandles;
00190     OptStatus optStatus;
00191     size_t mSteepest;
00192     double mAlpha;
00193     double maxAlpha;
00194     int pdgInd[3];
00195     ActiveSet mActive, testActive, originalActive;
00196 
00197     /* functions */
00198     void init_opt( PatchData& pd, MsqError& err );
00199     void init_max_step_length( PatchData& pd, MsqError& err );
00200 
00201     /* optimize */
00202     void minmax_opt( PatchData& pd, MsqError& err );
00203     void step_acceptance( PatchData& pd, int iter_count, const Vector3D& search_dir, MsqError& err );
00204     void get_min_estimate( double* final_est, MsqError& err );
00205     void get_gradient_projections( const Vector3D& mSearch, MsqError& err );
00206     void compute_alpha( MsqError& err );
00207 
00208     /* function/gradient/active set computations */
00209     bool compute_function( PatchData* pd, std::vector< double >& function, MsqError& err );
00210     bool compute_gradient( PatchData* pd, std::vector< Vector3D >& grad_out, MsqError& err );
00211     void find_active_set( const std::vector< double >& function, ActiveSet& active_set );
00212     void print_active_set( const ActiveSet& active_set, const std::vector< double >& func );
00213 
00214     /* checking validity/improvement */
00215     bool improvement_check( MsqError& err );
00216     bool validity_check( PatchData& pd, MsqError& err );
00217 
00218     /* checking equilibrium routines */
00219     bool check_equilibrium( OptStatus& opt_status, MsqError& err );
00220     bool convex_hull_test( const std::vector< Vector3D >& vec, MsqError& err );
00221     bool check_vector_dots( const std::vector< Vector3D >& vec, const Vector3D& normal, MsqError& err );
00222     void find_plane_normal( const Vector3D& pt1, const Vector3D& pt2, const Vector3D& pt3, Vector3D& cross,
00223                             MsqError& /*err*/ );
00224     void find_plane_points( Direction dir1, Direction dir2, const std::vector< Vector3D >& vec, Vector3D& pt1,
00225                             Vector3D& pt2, Vector3D& pt3, Status& status, MsqError& );
00226 
00227     /* from the matrix file */
00228     void form_grammian( const std::vector< Vector3D >& vec, MsqError& err );
00229     void form_PD_grammian( MsqError& err );
00230     void singular_test( int n, const Matrix3D& A, bool& singular, MsqError& err );
00231     bool solve2x2( double a11, double a12, double a21, double a22, double b1, double b2, double x[2], MsqError& err );
00232     void form_reduced_matrix( SymmetricMatrix& P );
00233 
00234     /* search direction */
00235     void search_direction( PatchData& pd, Vector3D& search_dir_out, MsqError& err );
00236     void search_edges_faces( const Vector3D* dir, Vector3D& search, MsqError& err );
00237     void get_active_directions( const std::vector< Vector3D >& gradient, std::vector< Vector3D >& dir, MsqError& err );
00238     NonSmoothDescent( const NonSmoothDescent& pd );             // disable copying
00239     NonSmoothDescent& operator=( const NonSmoothDescent& pd );  // disable assignment
00240 };
00241 
00242 inline void NonSmoothDescent::find_plane_normal( const Vector3D& pt1, const Vector3D& pt2, const Vector3D& pt3,
00243                                                  Vector3D& cross, MsqError& /*err*/ )
00244 {
00245     Vector3D vecA = pt2 - pt1;
00246     Vector3D vecB = pt3 - pt1;
00247     cross         = vecA * vecB;
00248     cross /= cross.length();
00249 }
00250 
00251 inline bool NonSmoothDescent::improvement_check( MsqError& /*err*/ )
00252 {
00253     /* check to see that the mesh didn't get worse */
00254     if( originalValue < mActive.true_active_value )
00255     {
00256         MSQ_PRINT( 2 )
00257         ( "The local mesh got worse; initial value %f; final value %f\n", originalValue, mActive.true_active_value );
00258         return false;
00259     }
00260 
00261     return true;
00262 }
00263 
00264 }  // namespace MBMesquite
00265 
00266 #endif  // Mesquite_NonSmoothDescent_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines