Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : /*!
28 : : \file NonSmoothDescent.hpp
29 : : \brief
30 : :
31 : : The NonSmoothDescent Class implements the steepest descent algorythm in
32 : : order to move a free vertex to an optimal position given an
33 : : ObjectiveFunction object and a QaulityMetric object.
34 : :
35 : : \author Thomas Leurent
36 : : \date 2002-06-13
37 : : */
38 : :
39 : : #ifndef Mesquite_NonSmoothDescent_hpp
40 : : #define Mesquite_NonSmoothDescent_hpp
41 : :
42 : : #include "Mesquite.hpp"
43 : : #include "VertexMover.hpp"
44 : : #include "ObjectiveFunction.hpp"
45 : : #include "MsqFreeVertexIndexIterator.hpp"
46 : : #include "MsqDebug.hpp"
47 : : #include "QualityMetric.hpp"
48 : : #include "VertexPatches.hpp"
49 : : #include <vector>
50 : : #include "Vector3D.hpp"
51 : :
52 : : namespace MBMesquite
53 : : {
54 : :
55 : : class NonSmoothDescent : public VertexMover
56 : : {
57 : : public:
58 : : MESQUITE_EXPORT
59 : : NonSmoothDescent( QualityMetric* qm );
60 : :
61 [ - + ]: 2 : MESQUITE_EXPORT virtual ~NonSmoothDescent() {}
62 : :
63 : : MESQUITE_EXPORT virtual std::string get_name() const;
64 : :
65 : : MESQUITE_EXPORT virtual PatchSet* get_patch_set();
66 : :
67 : : protected:
68 : 70218 : struct ActiveSet
69 : : {
70 : : int num_equal;
71 : : double true_active_value;
72 : : std::vector< int > active_ind;
73 : :
74 : 100495 : void set( int index )
75 : : {
76 : 100495 : active_ind.clear();
77 : 100495 : active_ind.push_back( index );
78 : 100495 : num_equal = 0;
79 : 100495 : }
80 : 23171 : void add( int index, bool equal )
81 : : {
82 : 23171 : active_ind.push_back( index );
83 : 23171 : num_equal += equal;
84 : 23171 : }
85 : : };
86 : :
87 : : class SymmetricMatrix
88 : : {
89 : : double* storage;
90 : : size_t size;
91 : 221646 : size_t index( size_t r, size_t c ) const
92 : : {
93 [ + + ]: 221646 : return ( r <= c ) ? size * r - r * ( r + 1 ) / 2 + c : size * c - c * ( c + 1 ) / 2 + r;
94 : : }
95 : :
96 : : public:
97 : 13464 : SymmetricMatrix() : storage( 0 ), size( 0 ) {}
98 : :
99 : 13464 : ~SymmetricMatrix()
100 : : {
101 : 13464 : free( storage );
102 : 13464 : }
103 : :
104 : 11789 : void resize( size_t new_size )
105 : : {
106 : 11789 : size = new_size;
107 : 11789 : storage = (double*)realloc( storage, sizeof( double ) * size * ( size + 1 ) / 2 );
108 : 11789 : }
109 : :
110 : 1322 : void fill( double value )
111 : : {
112 : 1322 : std::fill( storage, storage + ( size * ( size + 1 ) / 2 ), value );
113 : 1322 : }
114 : :
115 : 194727 : double& operator()( size_t r, size_t c )
116 : : {
117 : 194727 : return storage[index( r, c )];
118 : : }
119 : :
120 : 26919 : double operator()( size_t r, size_t c ) const
121 : : {
122 : 26919 : return storage[index( r, c )];
123 : : }
124 : :
125 : : double condition3x3() const;
126 : : };
127 : :
128 : : MESQUITE_EXPORT virtual void initialize( PatchData& pd, MsqError& err );
129 : :
130 : : MESQUITE_EXPORT virtual void optimize_vertex_positions( PatchData& pd, MsqError& err );
131 : :
132 : : MESQUITE_EXPORT virtual void initialize_mesh_iteration( PatchData& pd, MsqError& err );
133 : :
134 : : MESQUITE_EXPORT virtual void terminate_mesh_iteration( PatchData& pd, MsqError& err );
135 : :
136 : : MESQUITE_EXPORT virtual void cleanup();
137 : :
138 : : private:
139 : : enum OptStatus
140 : : {
141 : : MSQ_NO_STATUS = 0,
142 : : MSQ_STEP_ACCEPTED = 100,
143 : : MSQ_IMP_TOO_SMALL,
144 : : MSQ_FLAT_NO_IMP,
145 : : MSQ_STEP_TOO_SMALL,
146 : : MSQ_EQUILIBRIUM,
147 : : MSQ_ZERO_SEARCH,
148 : : MSQ_MAX_ITER_EXCEEDED
149 : : };
150 : :
151 : : enum Status
152 : : {
153 : : MSQ_NO_EQUIL = 101,
154 : : MSQ_CHECK_TOP_DOWN,
155 : : MSQ_CHECK_BOTTOM_UP,
156 : : MSQ_TWO_PT_PLANE,
157 : : MSQ_THREE_PT_PLANE,
158 : : MSQ_CHECK_Y_COORD_DIRECTION,
159 : : MSQ_CHECK_X_COORD_DIRECTION,
160 : : MSQ_CHECK_Z_COORD_DIRECTION,
161 : : MSQ_EQUIL,
162 : : MSQ_HULL_TEST_ERROR
163 : : };
164 : :
165 : : enum Direction
166 : : {
167 : : MSQ_XDIR = 0,
168 : : MSQ_YDIR = 1,
169 : : MSQ_ZDIR = 2
170 : : };
171 : :
172 : : /* local copy of patch data */
173 : : // PatchData patch_data;
174 : : size_t freeVertexIndex;
175 : :
176 : : /* smoothing parameters */
177 : : double minStepSize;
178 : :
179 : : /* optimization data */
180 : : VertexPatches patchSet;
181 : : QualityMetric* currentQM;
182 : : double originalValue;
183 : : std::vector< double > mFunction, testFunction, originalFunction;
184 : : std::vector< double > mGS;
185 : : SymmetricMatrix mG;
186 : : Matrix3D mPDG;
187 : : std::vector< double > prevActiveValues;
188 : : std::vector< Vector3D > mGradient, tmpGrad;
189 : : std::vector< size_t > tmpIdx, qmHandles;
190 : : OptStatus optStatus;
191 : : size_t mSteepest;
192 : : double mAlpha;
193 : : double maxAlpha;
194 : : int pdgInd[3];
195 : : ActiveSet mActive, testActive, originalActive;
196 : :
197 : : /* functions */
198 : : void init_opt( PatchData& pd, MsqError& err );
199 : : void init_max_step_length( PatchData& pd, MsqError& err );
200 : :
201 : : /* optimize */
202 : : void minmax_opt( PatchData& pd, MsqError& err );
203 : : void step_acceptance( PatchData& pd, int iter_count, const Vector3D& search_dir, MsqError& err );
204 : : void get_min_estimate( double* final_est, MsqError& err );
205 : : void get_gradient_projections( const Vector3D& mSearch, MsqError& err );
206 : : void compute_alpha( MsqError& err );
207 : :
208 : : /* function/gradient/active set computations */
209 : : bool compute_function( PatchData* pd, std::vector< double >& function, MsqError& err );
210 : : bool compute_gradient( PatchData* pd, std::vector< Vector3D >& grad_out, MsqError& err );
211 : : void find_active_set( const std::vector< double >& function, ActiveSet& active_set );
212 : : void print_active_set( const ActiveSet& active_set, const std::vector< double >& func );
213 : :
214 : : /* checking validity/improvement */
215 : : bool improvement_check( MsqError& err );
216 : : bool validity_check( PatchData& pd, MsqError& err );
217 : :
218 : : /* checking equilibrium routines */
219 : : bool check_equilibrium( OptStatus& opt_status, MsqError& err );
220 : : bool convex_hull_test( const std::vector< Vector3D >& vec, MsqError& err );
221 : : bool check_vector_dots( const std::vector< Vector3D >& vec, const Vector3D& normal, MsqError& err );
222 : : void find_plane_normal( const Vector3D& pt1, const Vector3D& pt2, const Vector3D& pt3, Vector3D& cross,
223 : : MsqError& /*err*/ );
224 : : void find_plane_points( Direction dir1, Direction dir2, const std::vector< Vector3D >& vec, Vector3D& pt1,
225 : : Vector3D& pt2, Vector3D& pt3, Status& status, MsqError& );
226 : :
227 : : /* from the matrix file */
228 : : void form_grammian( const std::vector< Vector3D >& vec, MsqError& err );
229 : : void form_PD_grammian( MsqError& err );
230 : : void singular_test( int n, const Matrix3D& A, bool& singular, MsqError& err );
231 : : bool solve2x2( double a11, double a12, double a21, double a22, double b1, double b2, double x[2], MsqError& err );
232 : : void form_reduced_matrix( SymmetricMatrix& P );
233 : :
234 : : /* search direction */
235 : : void search_direction( PatchData& pd, Vector3D& search_dir_out, MsqError& err );
236 : : void search_edges_faces( const Vector3D* dir, Vector3D& search, MsqError& err );
237 : : void get_active_directions( const std::vector< Vector3D >& gradient, std::vector< Vector3D >& dir, MsqError& err );
238 : : NonSmoothDescent( const NonSmoothDescent& pd ); // disable copying
239 : : NonSmoothDescent& operator=( const NonSmoothDescent& pd ); // disable assignment
240 : : };
241 : :
242 : 5694 : inline void NonSmoothDescent::find_plane_normal( const Vector3D& pt1, const Vector3D& pt2, const Vector3D& pt3,
243 : : Vector3D& cross, MsqError& /*err*/ )
244 : : {
245 [ + - ]: 5694 : Vector3D vecA = pt2 - pt1;
246 [ + - ]: 5694 : Vector3D vecB = pt3 - pt1;
247 [ + - ][ + - ]: 5694 : cross = vecA * vecB;
248 [ + - ][ + - ]: 5694 : cross /= cross.length();
249 : 5694 : }
250 : :
251 : 40309 : inline bool NonSmoothDescent::improvement_check( MsqError& /*err*/ )
252 : : {
253 : : /* check to see that the mesh didn't get worse */
254 [ + + ]: 40309 : if( originalValue < mActive.true_active_value )
255 : : {
256 : : MSQ_PRINT( 2 )
257 : : ( "The local mesh got worse; initial value %f; final value %f\n", originalValue, mActive.true_active_value );
258 : 5954 : return false;
259 : : }
260 : :
261 : 34355 : return true;
262 : : }
263 : :
264 : : } // namespace MBMesquite
265 : :
266 : : #endif // Mesquite_NonSmoothDescent_hpp
|
