Branch data Line data Source code
1 : : #include "moab/DiscreteGeometry/HiReconstruction.hpp"
2 : : #include "moab/DiscreteGeometry/DGMSolver.hpp"
3 : : #include "moab/HalfFacetRep.hpp"
4 : :
5 : : #ifdef MOAB_HAVE_MPI
6 : : #include "moab/ParallelComm.hpp"
7 : : #endif
8 : : #include "MBTagConventions.hpp"
9 : :
10 : : #define HIREC_USE_AHF
11 : :
12 : : #include <math.h>
13 : : #include <deque>
14 : : #include <iostream>
15 : :
16 : : namespace moab
17 : : {
18 : 21 : HiReconstruction::HiReconstruction( Core* impl, ParallelComm* comm, EntityHandle meshIn, int minpnts, bool recwhole )
19 [ + - ][ + - ]: 21 : : mbImpl( impl ), pcomm( comm ), _mesh2rec( meshIn ), _MINPNTS( minpnts )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
20 : : {
21 [ - + ]: 21 : assert( NULL != impl );
22 : : ErrorCode error;
23 : 21 : _MINEPS = 1e-12;
24 : 21 : _dim = 0;
25 : 21 : _hasfittings = false;
26 : 21 : _hasderiv = false;
27 : : #ifdef MOAB_HAVE_MPI
28 [ + + ][ + - ]: 21 : if( !pcomm ) { pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 ); }
29 : : #endif
30 : :
31 [ + - ]: 21 : error = initialize( recwhole );
32 [ - + ]: 21 : if( MB_SUCCESS != error )
33 : : {
34 [ # # ][ # # ]: 0 : std::cout << "Error initializing HiReconstruction\n" << std::endl;
35 : 0 : exit( 1 );
36 : : }
37 : 21 : }
38 : :
39 : 42 : HiReconstruction::~HiReconstruction()
40 : : {
41 : : #ifdef MOAB_HAVE_AHF
42 : : ahf = NULL;
43 : : #else
44 [ + - ]: 21 : delete ahf;
45 : : #endif
46 : 21 : }
47 : :
48 : 21 : ErrorCode HiReconstruction::initialize( bool recwhole )
49 : : {
50 : : ErrorCode error;
51 : :
52 : : #ifdef HIREC_USE_AHF
53 : 21 : std::cout << "HIREC_USE_AHF: Initializing" << std::endl;
54 [ + - ]: 21 : ahf = new HalfFacetRep( mbImpl, pcomm, _mesh2rec, false );
55 [ - + ]: 21 : if( !ahf ) { return MB_MEMORY_ALLOCATION_FAILED; }
56 [ - + ][ # # ]: 21 : error = ahf->initialize();MB_CHK_ERR( error );
57 : : #else
58 : : ahf = NULL;
59 : : #endif
60 : :
61 : : // error = ahf->get_entity_ranges(_inverts,_inedges,_infaces,_incells); MB_CHK_ERR(error);
62 [ - + ][ # # ]: 21 : error = mbImpl->get_entities_by_dimension( _mesh2rec, 0, _inverts );MB_CHK_ERR( error );
63 [ - + ][ # # ]: 21 : error = mbImpl->get_entities_by_dimension( _mesh2rec, 1, _inedges );MB_CHK_ERR( error );
64 [ - + ][ # # ]: 21 : error = mbImpl->get_entities_by_dimension( _mesh2rec, 2, _infaces );MB_CHK_ERR( error );
65 [ - + ][ # # ]: 21 : error = mbImpl->get_entities_by_dimension( _mesh2rec, 3, _incells );MB_CHK_ERR( error );
66 [ + + ][ + - ]: 21 : if( _inedges.size() && _infaces.empty() && _incells.empty() )
[ + - ][ + + ]
67 : : {
68 : 4 : _dim = 1;
69 : 4 : _MAXPNTS = 13;
70 : : }
71 [ + - ][ + - ]: 17 : else if( _infaces.size() && _incells.empty() )
[ + - ]
72 : : {
73 : 17 : _dim = 2;
74 : 17 : _MAXPNTS = 128;
75 : : }
76 : : else
77 : : {
78 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Encountered a non-manifold mesh or a mesh with volume elements" );
[ # # ][ # # ]
[ # # ]
79 : : }
80 : :
81 : : // get locally hosted vertices by filtering pstatus
82 : : #ifdef MOAB_HAVE_MPI
83 [ + + ]: 21 : if( pcomm )
84 : : {
85 [ - + ][ # # ]: 5 : error = pcomm->filter_pstatus( _inverts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts2rec );MB_CHK_ERR( error );
86 : : }
87 : : else
88 : : {
89 : 16 : _verts2rec = _inverts;
90 : : }
91 : : #else
92 : : _verts2rec = _inverts;
93 : : #endif
94 : 21 : _nv2rec = _verts2rec.size();
95 : :
96 [ + - ]: 21 : if( recwhole )
97 : : {
98 : : // compute normals(surface) or tangent vector(curve) for all locally hosted vertices
99 [ + + ]: 21 : if( 2 == _dim ) { compute_average_vertex_normals_surf(); }
100 [ + - ]: 4 : else if( 1 == _dim )
101 : : {
102 : 4 : compute_average_vertex_tangents_curve();
103 : : }
104 : : else
105 : : {
106 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "Unknow space dimension" );
[ # # ][ # # ]
[ # # ]
107 : : }
108 : 21 : _hasderiv = true;
109 : : }
110 : 21 : return error;
111 : : }
112 : :
113 : : /***************************************************
114 : : * User Interface for Reconstruction of Geometry *
115 : : ***************************************************/
116 : :
117 : 159 : ErrorCode HiReconstruction::reconstruct3D_surf_geom( int degree, bool interp, bool safeguard, bool reset )
118 : : {
119 [ - + ]: 159 : assert( 2 == _dim );
120 [ + + ][ - + ]: 159 : if( _hasfittings && !reset )
121 : : {
122 : : // This object has precomputed fitting results and user don't want to reset
123 : 0 : return MB_SUCCESS;
124 : : }
125 : : else
126 : : {
127 : 159 : _initfittings = _hasfittings = false;
128 : : }
129 : : // initialize for geometric information
130 : 159 : initialize_surf_geom( degree );
131 : : ErrorCode error;
132 : : double *coeffs, *coords;
133 : : int* degree_out;
134 : 159 : int ncoeffs = ( degree + 2 ) * ( degree + 1 ) / 2;
135 : :
136 : : // DBG
137 : 159 : int dcount = 0;
138 : :
139 [ + - ][ + - ]: 67505 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
[ + - ][ + - ]
[ + + ]
140 : : {
141 [ + - ][ + - ]: 67346 : int index = _verts2rec.index( *ivert );
142 : : // for debug
143 : : /*if(index==70){
144 : : EntityHandle vid = *ivert;
145 : : double vertcoords[3];
146 : : error = mbImpl->get_coords(&vid,1,vertcoords);
147 : : }*/
148 : :
149 [ + - ]: 67346 : size_t istr = _vertID2coeffID[index];
150 [ + - ]: 67346 : coords = &( _local_coords[9 * index] );
151 [ + - ]: 67346 : coeffs = &( _local_fit_coeffs[istr] );
152 [ + - ]: 67346 : degree_out = &( _degrees_out[index] );
153 [ + - ]: 67346 : _interps[index] = interp;
154 [ + - ]: 67346 : error = polyfit3d_walf_surf_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 9, coords, degree_out, ncoeffs,
155 [ + - ][ - + ]: 67346 : coeffs );MB_CHK_ERR( error );
[ # # ][ # # ]
156 : :
157 : : // DBG
158 [ + + ]: 67346 : if( degree_out[0] < degree ) dcount += 1;
159 : : }
160 : :
161 : : // DBG
162 : : // std::cout<<"Total #points ="<<_verts2rec.size()<<", #degraded points = "<<dcount<<std::endl;
163 : :
164 : 159 : _geom = HISURFACE;
165 : 159 : _hasfittings = true;
166 : 159 : return error;
167 : : }
168 : :
169 : 0 : ErrorCode HiReconstruction::reconstruct3D_surf_geom( size_t npts, int* degrees, bool* interps, bool safeguard,
170 : : bool reset )
171 : : {
172 [ # # ]: 0 : assert( _dim == 2 );
173 [ # # ][ # # ]: 0 : if( npts != _nv2rec ) { MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match number of vertices" ); }
[ # # ][ # # ]
[ # # ][ # # ]
174 [ # # ][ # # ]: 0 : if( _hasfittings && !reset ) { return MB_SUCCESS; }
175 : : else
176 : : {
177 : 0 : _initfittings = _hasfittings = false;
178 : : }
179 : : ErrorCode error;
180 : : // initialization for fitting
181 : 0 : initialize_surf_geom( npts, degrees );
182 : : double *coeffs, *coords;
183 : : int* degree_out;
184 : 0 : size_t i = 0;
185 [ # # ][ # # ]: 0 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
[ # # ][ # # ]
[ # # ]
186 : : {
187 [ # # ][ # # ]: 0 : int index = _verts2rec.index( *ivert );
188 [ # # ]: 0 : assert( -1 != index );
189 [ # # ]: 0 : size_t istr = _vertID2coeffID[index];
190 [ # # ]: 0 : coords = &( _local_coords[9 * index] );
191 [ # # ]: 0 : coeffs = &( _local_fit_coeffs[istr] );
192 [ # # ]: 0 : degree_out = &( _degrees_out[index] );
193 [ # # ]: 0 : _interps[index] = interps[i];
194 : 0 : int ncoeffs = ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
195 [ # # ]: 0 : error = polyfit3d_walf_surf_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 9, coords, degree_out,
196 [ # # ][ # # ]: 0 : ncoeffs, coeffs );MB_CHK_ERR( error );
[ # # ][ # # ]
197 : : }
198 : 0 : _geom = HISURFACE;
199 : 0 : _hasfittings = true;
200 : 0 : return error;
201 : : }
202 : :
203 : 48 : ErrorCode HiReconstruction::reconstruct3D_curve_geom( int degree, bool interp, bool safeguard, bool reset )
204 : : {
205 [ - + ]: 48 : assert( _dim == 1 );
206 [ + + ][ - + ]: 48 : if( _hasfittings && !reset ) { return MB_SUCCESS; }
207 : : else
208 : : {
209 : 48 : _initfittings = _hasfittings = false;
210 : : }
211 : 48 : initialize_3Dcurve_geom( degree );
212 : : ErrorCode error;
213 : 48 : double *coords = 0, *coeffs;
214 : : int* degree_out;
215 : 48 : int ncoeffs = 3 * ( degree + 1 );
216 [ + - ][ + - ]: 492 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert )
[ + - ][ + - ]
[ + + ]
217 : : {
218 [ + - ][ + - ]: 444 : int index = _verts2rec.index( *ivert );
219 [ - + ]: 444 : assert( index != -1 );
220 [ + - ]: 444 : size_t istr = _vertID2coeffID[index];
221 [ + - ]: 444 : coeffs = &( _local_fit_coeffs[istr] );
222 [ + - ]: 444 : degree_out = &( _degrees_out[index] );
223 [ + - ]: 444 : _interps[index] = interp;
224 [ + - ]: 444 : error = polyfit3d_walf_curve_vertex( *ivert, interp, degree, _MINPNTS, safeguard, 0, coords, degree_out,
225 [ + - ][ - + ]: 444 : ncoeffs, coeffs );MB_CHK_ERR( error );
[ # # ][ # # ]
226 : : }
227 : 48 : _geom = HI3DCURVE;
228 : 48 : _hasfittings = true;
229 : 48 : return error;
230 : : }
231 : :
232 : 0 : ErrorCode HiReconstruction::reconstruct3D_curve_geom( size_t npts, int* degrees, bool* interps, bool safeguard,
233 : : bool reset )
234 : : {
235 [ # # ]: 0 : assert( _dim == 1 );
236 : : ErrorCode error;
237 [ # # ][ # # ]: 0 : if( npts != _nv2rec ) { MB_SET_ERR( MB_FAILURE, "Input number of degrees doesn't match the number of vertices" ); }
[ # # ][ # # ]
[ # # ][ # # ]
238 [ # # ][ # # ]: 0 : if( _hasfittings && !reset ) { return MB_SUCCESS; }
239 : : else
240 : : {
241 : 0 : _initfittings = _hasfittings = false;
242 : : }
243 : : // initialize
244 : 0 : initialize_3Dcurve_geom( npts, degrees );
245 : 0 : double *coords = 0, *coeffs;
246 : : int* degree_out;
247 : 0 : size_t i = 0;
248 [ # # ][ # # ]: 0 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++i )
[ # # ][ # # ]
[ # # ]
249 : : {
250 [ # # ][ # # ]: 0 : int index = _verts2rec.index( *ivert );
251 [ # # ]: 0 : size_t istr = _vertID2coeffID[index];
252 [ # # ]: 0 : coeffs = &( _local_fit_coeffs[istr] );
253 [ # # ]: 0 : degree_out = &( _degrees_out[index] );
254 [ # # ]: 0 : _interps[index] = interps[i];
255 : 0 : int ncoeffs = 3 * ( degrees[i] + 1 );
256 [ # # ]: 0 : error = polyfit3d_walf_curve_vertex( *ivert, interps[i], degrees[i], _MINPNTS, safeguard, 0, coords, degree_out,
257 [ # # ][ # # ]: 0 : ncoeffs, coeffs );MB_CHK_ERR( error );
[ # # ][ # # ]
258 : : }
259 : 0 : _geom = HI3DCURVE;
260 : 0 : _hasfittings = true;
261 : 0 : return error;
262 : : }
263 : :
264 : 67346 : ErrorCode HiReconstruction::polyfit3d_walf_surf_vertex( const EntityHandle vid, const bool interp, int degree,
265 : : int minpnts, const bool safeguard, const int ncoords,
266 : : double* coords, int* degree_out, const int ncoeffs,
267 : : double* coeffs )
268 : : {
269 [ - + ]: 67346 : assert( _dim == 2 );
270 : : ErrorCode error;
271 [ + - ]: 67346 : int ring = estimate_num_rings( degree, interp );
272 : : // std::cout<<"ring = "<<ring<<std::endl;
273 : : // get n-ring neighbors
274 [ + - ]: 67346 : Range ngbvs;
275 [ + - ][ - + ]: 67346 : error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
[ # # ][ # # ]
276 : : // for debug
277 : : /*if(_verts2rec.index(vid)==70){
278 : : for(Range::iterator ingb=ngbvs.begin();ingb!=ngbvs.end();++ingb) std::cerr <<
279 : : _verts2rec.index(*ingb) << " "; std::cout << std::endl;
280 : : }*/
281 : :
282 : : // get coordinates;
283 [ + - ]: 67346 : size_t nverts = ngbvs.size();
284 [ - + ]: 67346 : assert( nverts );
285 [ + - ][ + - ]: 67346 : double* ngbcoords = new double[nverts * 3];
286 [ + - ][ - + ]: 67346 : error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
[ # # ][ # # ]
287 : : // get normals
288 [ + - ][ + - ]: 67346 : double* ngbnrms = new double[nverts * 3];
289 [ + - ][ - + ]: 67346 : error = get_normals_surf( ngbvs, ngbnrms );MB_CHK_ERR( error );
[ # # ][ # # ]
290 : : // switch vid to first one
291 [ + - ]: 67346 : int index = ngbvs.index( vid );
292 [ - + ]: 67346 : assert( index != -1 );
293 : 67346 : std::swap( ngbcoords[0], ngbcoords[3 * index] );
294 : 67346 : std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
295 : 67346 : std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
296 : 67346 : std::swap( ngbnrms[0], ngbnrms[3 * index] );
297 : 67346 : std::swap( ngbnrms[1], ngbnrms[3 * index + 1] );
298 : 67346 : std::swap( ngbnrms[2], ngbnrms[3 * index + 2] );
299 : : // local WLS fitting
300 : : int degree_pnt, degree_qr;
301 : : polyfit3d_surf_get_coeff( nverts, ngbcoords, ngbnrms, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
302 [ + - ]: 67346 : degree_out, °ree_pnt, °ree_qr );
303 [ + - ]: 67346 : delete[] ngbcoords;
304 [ + - ]: 67346 : delete[] ngbnrms;
305 : 67346 : return error;
306 : : }
307 : :
308 : 444 : ErrorCode HiReconstruction::polyfit3d_walf_curve_vertex( const EntityHandle vid, const bool interp, int degree,
309 : : int minpnts, const bool safeguard, const int ncoords,
310 : : double* coords, int* degree_out, const int ncoeffs,
311 : : double* coeffs )
312 : : {
313 : : ErrorCode error;
314 [ + - ]: 444 : int ring = estimate_num_rings( degree, interp );
315 : : // get n-ring neighbors
316 [ + - ]: 444 : Range ngbvs;
317 [ + - ][ - + ]: 444 : error = obtain_nring_ngbvs( vid, ring, minpnts, ngbvs );MB_CHK_ERR( error );
[ # # ][ # # ]
318 : : // get coordinates
319 [ + - ]: 444 : size_t nverts = ngbvs.size();
320 [ - + ]: 444 : assert( nverts );
321 [ + - ][ + - ]: 444 : double* ngbcoords = new double[nverts * 3];
322 [ + - ][ - + ]: 444 : error = mbImpl->get_coords( ngbvs, ngbcoords );MB_CHK_ERR( error );
[ # # ][ # # ]
323 : : // get tangent vectors
324 [ + - ][ + - ]: 444 : double* ngbtangs = new double[nverts * 3];
325 [ + - ][ - + ]: 444 : error = get_tangents_curve( ngbvs, ngbtangs );MB_CHK_ERR( error );
[ # # ][ # # ]
326 : : // switch vid to first one
327 [ + - ]: 444 : int index = ngbvs.index( vid );
328 [ - + ]: 444 : assert( index != -1 );
329 : 444 : std::swap( ngbcoords[0], ngbcoords[3 * index] );
330 : 444 : std::swap( ngbcoords[1], ngbcoords[3 * index + 1] );
331 : 444 : std::swap( ngbcoords[2], ngbcoords[3 * index + 2] );
332 : 444 : std::swap( ngbtangs[0], ngbtangs[3 * index] );
333 : 444 : std::swap( ngbtangs[1], ngbtangs[3 * index + 1] );
334 : 444 : std::swap( ngbtangs[2], ngbtangs[3 * index + 2] );
335 : : // local WLS fittings
336 : : polyfit3d_curve_get_coeff( nverts, ngbcoords, ngbtangs, degree, interp, safeguard, ncoords, coords, ncoeffs, coeffs,
337 [ + - ]: 444 : degree_out );
338 [ + - ]: 444 : delete[] ngbcoords;
339 [ + - ]: 444 : delete[] ngbtangs;
340 : 444 : return error;
341 : : }
342 : :
343 : : /**************************************************************
344 : : * User Interface for Evaluation via Reconstructed Geometry *
345 : : **************************************************************/
346 : :
347 : 104212 : ErrorCode HiReconstruction::hiproj_walf_in_element( EntityHandle elem, const int nvpe, const int npts2fit,
348 : : const double* naturalcoords2fit, double* newcoords )
349 : : {
350 [ - + ]: 104212 : assert( newcoords );
351 : : ErrorCode error;
352 : : // get connectivity table
353 [ + - ]: 104212 : std::vector< EntityHandle > elemconn;
354 [ + - ][ - + ]: 104212 : error = mbImpl->get_connectivity( &elem, 1, elemconn );MB_CHK_ERR( error );
[ # # ][ # # ]
355 [ - + ]: 104212 : if( nvpe != (int)elemconn.size() )
356 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "element connectivity table size doesn't match input size" ); }
[ # # ][ # # ]
[ # # ]
357 : :
358 [ - + ][ # # ]: 104212 : if( !_hasfittings ) { MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" ); }
[ # # ][ # # ]
[ # # ][ # # ]
359 : : else
360 : : {
361 [ + - ]: 104212 : std::ostringstream convert;
362 [ + - ]: 104212 : convert << elem;
363 [ + - ][ + - ]: 208424 : std::string ID = convert.str();
364 [ + + ][ + - ]: 448684 : for( int i = 0; i < nvpe; ++i )
365 : : {
366 [ + - ][ + - ]: 344472 : if( -1 == _verts2rec.index( elemconn[i] ) )
[ - + ]
367 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for element " + ID ); }
[ # # ][ # # ]
[ # # ][ # # ]
368 : 104212 : }
369 : : }
370 : : // check correctness of input
371 [ + + ]: 685784 : for( int i = 0; i < npts2fit; ++i )
372 : : {
373 [ + - ][ - + ]: 581572 : if( !check_barycentric_coords( nvpe, naturalcoords2fit + i * nvpe ) )
374 [ # # ][ # # ]: 0 : { MB_SET_ERR( MB_FAILURE, "Wrong barycentric coordinates" ); }
[ # # ][ # # ]
[ # # ]
375 : : }
376 : :
377 [ + - ][ + - ]: 104212 : double* elemcoords = new double[nvpe * 3];
378 [ + - ][ + - ]: 104212 : error = mbImpl->get_coords( &( elemconn[0] ), nvpe, elemcoords );MB_CHK_ERR( error );
[ - + ][ # # ]
[ # # ]
379 : :
380 [ + - ][ + - ]: 1848928 : double* coords2fit = new double[3 * npts2fit]();
[ + + ]
381 [ + + ]: 685784 : for( int i = 0; i < npts2fit; ++i )
382 : : {
383 [ + + ]: 2358124 : for( int j = 0; j < nvpe; ++j )
384 : : {
385 : 1776552 : coords2fit[3 * i] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j];
386 : 1776552 : coords2fit[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 1];
387 : 1776552 : coords2fit[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * elemcoords[3 * j + 2];
388 : : }
389 : : }
390 : :
391 [ + - ][ + - ]: 104212 : double* hiproj_new = new double[3 * npts2fit];
392 : : // initialize output
393 [ + + ]: 685784 : for( int i = 0; i < npts2fit; ++i )
394 : : {
395 : 581572 : newcoords[3 * i] = newcoords[3 * i + 1] = newcoords[3 * i + 2] = 0;
396 : : }
397 : : // for each input vertex, call nvpe fittings and take average
398 [ + + ]: 448684 : for( int j = 0; j < nvpe; ++j )
399 : : {
400 [ + - ][ + - ]: 344472 : error = hiproj_walf_around_vertex( elemconn[j], npts2fit, coords2fit, hiproj_new );MB_CHK_ERR( error );
[ - + ][ # # ]
[ # # ]
401 [ + + ]: 2121024 : for( int i = 0; i < npts2fit; ++i )
402 : : {
403 : 1776552 : newcoords[3 * i] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i];
404 : 1776552 : newcoords[3 * i + 1] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 1];
405 : 1776552 : newcoords[3 * i + 2] += naturalcoords2fit[i * nvpe + j] * hiproj_new[3 * i + 2];
406 : : }
407 : : }
408 [ + - ]: 104212 : delete[] elemcoords;
409 [ + - ]: 104212 : delete[] coords2fit;
410 [ + - ]: 104212 : delete[] hiproj_new;
411 : 104212 : return error;
412 : : }
413 : :
414 : 344472 : ErrorCode HiReconstruction::hiproj_walf_around_vertex( EntityHandle vid, const int npts2fit, const double* coords2fit,
415 : : double* hiproj_new )
416 : : {
417 [ - + ][ # # ]: 344472 : if( !_hasfittings ) { MB_SET_ERR( MB_FAILURE, "There is no existing fitting results" ); }
[ # # ][ # # ]
[ # # ][ # # ]
418 [ + - ][ - + ]: 344472 : else if( -1 == _verts2rec.index( vid ) )
419 : : {
420 [ # # ]: 0 : std::ostringstream convert;
421 [ # # ]: 0 : convert << vid;
422 [ # # ]: 0 : std::string VID = convert.str();
423 [ # # ][ # # ]: 0 : MB_SET_ERR( MB_FAILURE, "There is no existing fitting results for vertex " + VID );
[ # # ][ # # ]
[ # # ][ # # ]
424 : : }
425 : : ErrorCode error;
426 : : // get center of local coordinates system
427 : : double local_origin[3];
428 [ + - ][ - + ]: 344472 : error = mbImpl->get_coords( &vid, 1, local_origin );MB_CHK_ERR( error );
[ # # ][ # # ]
429 : : // get local fitting parameters
430 [ + - ]: 344472 : int index = _verts2rec.index( vid );
431 [ + - ]: 344472 : bool interp = _interps[index];
432 [ + - ]: 344472 : int local_deg = _degrees_out[index];
433 : : double *uvw_coords, *local_coeffs;
434 [ + + ]: 344472 : if( _geom == HISURFACE )
435 : : {
436 [ + - ]: 343584 : uvw_coords = &( _local_coords[9 * index] );
437 : : // int ncoeffs = (local_deg+2)*(local_deg+1)>>1;
438 [ + - ]: 343584 : size_t istr = _vertID2coeffID[index];
439 [ + - ]: 343584 : local_coeffs = &( _local_fit_coeffs[istr] );
440 : : walf3d_surf_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
441 [ + - ]: 343584 : hiproj_new );
442 : : }
443 [ + - ]: 888 : else if( _geom == HI3DCURVE )
444 : : {
445 [ + - ]: 888 : uvw_coords = &( _local_coords[3 * index] );
446 [ + - ]: 888 : size_t istr = _vertID2coeffID[index];
447 [ + - ]: 888 : local_coeffs = &( _local_fit_coeffs[istr] );
448 : : walf3d_curve_vertex_eval( local_origin, uvw_coords, local_deg, local_coeffs, interp, npts2fit, coords2fit,
449 [ + - ]: 888 : hiproj_new );
450 : : }
451 : 344472 : return error;
452 : : }
453 : :
454 : 343584 : void HiReconstruction::walf3d_surf_vertex_eval( const double* local_origin, const double* local_coords,
455 : : const int local_deg, const double* local_coeffs, const bool interp,
456 : : const int npts2fit, const double* coords2fit, double* hiproj_new )
457 : : {
458 : : double xaxis[3], yaxis[3], zaxis[3];
459 [ + + ]: 1374336 : for( int i = 0; i < 3; ++i )
460 : : {
461 : 1030752 : xaxis[i] = local_coords[i];
462 : 1030752 : yaxis[i] = local_coords[3 + i];
463 : 1030752 : zaxis[i] = local_coords[6 + i];
464 : : }
465 : : // double *basis = new double[(local_deg+2)*(local_deg+1)/2-1];
466 [ + - ]: 343584 : std::vector< double > basis( ( local_deg + 2 ) * ( local_deg + 1 ) / 2 - 1 );
467 [ + + ]: 2119248 : for( int i = 0; i < npts2fit; ++i )
468 : : {
469 : : double local_pos[3];
470 [ + + ]: 7102656 : for( int j = 0; j < 3; ++j )
471 : : {
472 : 5326992 : local_pos[j] = coords2fit[3 * i + j] - local_origin[j];
473 : : }
474 : 1775664 : double u, v, height = 0;
475 [ + - ]: 1775664 : u = DGMSolver::vec_innerprod( 3, local_pos, xaxis );
476 [ + - ]: 1775664 : v = DGMSolver::vec_innerprod( 3, local_pos, yaxis );
477 [ + - ]: 1775664 : basis[0] = u;
478 [ + - ]: 1775664 : basis[1] = v;
479 : 1775664 : int l = 1;
480 [ + + ]: 6194696 : for( int k = 2; k <= local_deg; ++k )
481 : : {
482 : 4419032 : ++l;
483 [ + - ][ + - ]: 4419032 : basis[l] = u * basis[l - k];
484 [ + + ]: 19126148 : for( int id = 0; id < k; ++id )
485 : : {
486 : 14707116 : ++l;
487 [ + - ][ + - ]: 14707116 : basis[l] = basis[l - k - 1] * v;
488 : : }
489 : : }
490 [ + + ]: 1775664 : if( !interp ) { height = local_coeffs[0]; }
491 [ + + ]: 24453140 : for( int p = 0; p <= l; ++p )
492 : : {
493 [ + - ]: 22677476 : height += local_coeffs[p + 1] * basis[p];
494 : : }
495 : 1775664 : hiproj_new[3 * i] = local_origin[0] + u * xaxis[0] + v * yaxis[0] + height * zaxis[0];
496 : 1775664 : hiproj_new[3 * i + 1] = local_origin[1] + u * xaxis[1] + v * yaxis[1] + height * zaxis[1];
497 : 1775664 : hiproj_new[3 * i + 2] = local_origin[2] + u * xaxis[2] + v * yaxis[2] + height * zaxis[2];
498 : 343584 : }
499 : : // delete [] basis;
500 : 343584 : }
501 : :
502 : 888 : void HiReconstruction::walf3d_curve_vertex_eval( const double* local_origin, const double* local_coords,
503 : : const int local_deg, const double* local_coeffs, const bool interp,
504 : : const int npts2fit, const double* coords2fit, double* hiproj_new )
505 : : {
506 [ + - ][ + - ]: 888 : assert( local_origin && local_coords && local_coeffs );
[ - + ]
507 : 888 : int ncoeffspvpd = local_deg + 1;
508 [ + + ]: 1776 : for( int i = 0; i < npts2fit; ++i )
509 : : {
510 : : // get the vector from center to current point, project to tangent line
511 : 888 : double vec[3], ans[3] = { 0, 0, 0 };
512 [ + - ]: 888 : DGMSolver::vec_linear_operation( 3, 1, coords2fit + 3 * i, -1, local_origin, vec );
513 [ + - ]: 888 : double u = DGMSolver::vec_innerprod( 3, local_coords, vec );
514 : : // evaluate polynomials
515 [ + + ]: 888 : if( !interp )
516 : : {
517 : 444 : ans[0] = local_coeffs[0];
518 : 444 : ans[1] = local_coeffs[ncoeffspvpd];
519 : 444 : ans[2] = local_coeffs[2 * ncoeffspvpd];
520 : : }
521 : 888 : double uk = 1; // degree_out and degree different, stored in columnwise contiguously
522 [ + + ]: 3288 : for( int j = 1; j < ncoeffspvpd; ++j )
523 : : {
524 : 2400 : uk *= u;
525 : 2400 : ans[0] += uk * local_coeffs[j];
526 : 2400 : ans[1] += uk * local_coeffs[j + ncoeffspvpd];
527 : 2400 : ans[2] += uk * local_coeffs[j + 2 * ncoeffspvpd];
528 : : }
529 : 888 : hiproj_new[3 * i] = ans[0] + local_origin[0];
530 : 888 : hiproj_new[3 * i + 1] = ans[1] + local_origin[1];
531 : 888 : hiproj_new[3 * i + 2] = ans[2] + local_origin[2];
532 : : }
533 : 888 : }
534 : :
535 : 0 : bool HiReconstruction::get_fittings_data( EntityHandle vid, GEOMTYPE& geomtype, std::vector< double >& coords,
536 : : int& degree_out, std::vector< double >& coeffs, bool& interp )
537 : : {
538 [ # # ]: 0 : if( !_hasfittings ) { return false; }
539 : : else
540 : : {
541 : 0 : int index = _verts2rec.index( vid );
542 [ # # ]: 0 : if( -1 == index )
543 : : {
544 : 0 : std::cout << "Input vertex is not locally hosted vertex in this mesh set" << std::endl;
545 : 0 : return false;
546 : : }
547 : 0 : geomtype = _geom;
548 [ # # ]: 0 : if( HISURFACE == _geom )
549 : : {
550 [ # # ][ # # ]: 0 : coords.insert( coords.end(), _local_coords.begin() + 9 * index, _local_coords.begin() + 9 * index + 9 );
[ # # ][ # # ]
551 : 0 : degree_out = _degrees_out[index];
552 : 0 : interp = _interps[index];
553 : 0 : int ncoeffs = ( degree_out + 2 ) * ( degree_out + 1 ) >> 1;
554 : 0 : size_t istr = _vertID2coeffID[index];
555 [ # # ][ # # ]: 0 : coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
[ # # ][ # # ]
556 : : }
557 [ # # ]: 0 : else if( HI3DCURVE == _geom )
558 : : {
559 [ # # ][ # # ]: 0 : coords.insert( coords.end(), _local_coords.begin() + 3 * index, _local_coords.begin() + 3 * index + 3 );
[ # # ][ # # ]
560 : 0 : degree_out = _degrees_out[index];
561 : 0 : interp = _interps[index];
562 : 0 : int ncoeffs = 3 * ( degree_out + 1 );
563 : 0 : size_t istr = _vertID2coeffID[index];
564 [ # # ][ # # ]: 0 : coeffs.insert( coeffs.end(), _local_fit_coeffs.begin() + istr, _local_fit_coeffs.begin() + istr + ncoeffs );
[ # # ][ # # ]
565 : : }
566 : 0 : return true;
567 : : }
568 : : }
569 : :
570 : : /****************************************************************
571 : : * Basic Internal Routines to initialize and set fitting data *
572 : : ****************************************************************/
573 : :
574 : 67790 : int HiReconstruction::estimate_num_rings( int degree, bool interp )
575 : : {
576 [ + + ]: 67790 : return interp ? ( ( degree + 1 ) >> 1 ) + ( ( degree + 1 ) & 1 ) : ( ( degree + 2 ) >> 1 ) + ( ( degree + 2 ) & 1 );
577 : : }
578 : :
579 : 1099548 : ErrorCode HiReconstruction::vertex_get_incident_elements( const EntityHandle& vid, const int elemdim,
580 : : std::vector< EntityHandle >& adjents )
581 : : {
582 : : ErrorCode error;
583 [ - + ]: 1099548 : assert( elemdim == _dim );
584 : : #ifdef HIREC_USE_AHF
585 [ - + ][ # # ]: 1099548 : error = ahf->get_up_adjacencies( vid, elemdim, adjents );MB_CHK_ERR( error );
586 : : #else
587 : : error = mbImpl->get_adjacencies( &vid, 1, elemdim, false, adjents );MB_CHK_ERR( error );
588 : : #endif
589 : 1099548 : return error;
590 : : }
591 : :
592 : 67790 : ErrorCode HiReconstruction::obtain_nring_ngbvs( const EntityHandle vid, int ring, const int minpnts, Range& ngbvs )
593 : : {
594 : : ErrorCode error;
595 [ + - ]: 67790 : std::deque< EntityHandle > todo;
596 [ + - ]: 67790 : todo.push_back( vid );
597 [ + - ]: 67790 : ngbvs.insert( vid );
598 : : EntityHandle pre, nxt;
599 [ + + ]: 260801 : for( int i = 1; i <= ring; ++i )
600 : : {
601 : 193170 : int count = todo.size();
602 [ + + ]: 1284802 : while( count )
603 : : {
604 [ + - ]: 1091632 : EntityHandle center = todo.front();
605 [ + - ]: 1091632 : todo.pop_front();
606 : 1091632 : --count;
607 [ + - ]: 1091632 : std::vector< EntityHandle > adjents;
608 [ + - ][ - + ]: 1091632 : error = vertex_get_incident_elements( center, _dim, adjents );MB_CHK_ERR( error );
[ # # ][ # # ]
609 [ + + ][ + - ]: 6808014 : for( size_t j = 0; j < adjents.size(); ++j )
610 : : {
611 [ + - ]: 5716382 : std::vector< EntityHandle > elemconn;
612 [ + - ][ + - ]: 5716382 : error = mbImpl->get_connectivity( &adjents[j], 1, elemconn );MB_CHK_ERR( error );
[ - + ][ # # ]
[ # # ]
613 : 5716382 : int nvpe = elemconn.size();
614 [ + - ][ + - ]: 17929184 : for( int k = 0; k < nvpe; ++k )
615 : : {
616 [ + - ][ + + ]: 12212802 : if( elemconn[k] == center )
617 : : {
618 [ + + ][ + - ]: 5716382 : pre = k == 0 ? elemconn[nvpe - 1] : elemconn[k - 1];
[ + - ]
619 [ + - ]: 5716382 : nxt = elemconn[( k + 1 ) % nvpe];
620 [ + - ][ + - ]: 5716382 : if( ngbvs.find( pre ) == ngbvs.end() )
[ + - ][ + + ]
621 : : {
622 [ + - ]: 985542 : ngbvs.insert( pre );
623 [ + - ]: 985542 : todo.push_back( pre );
624 : : }
625 [ + - ][ + - ]: 5716382 : if( ngbvs.find( nxt ) == ngbvs.end() )
[ + - ][ + + ]
626 : : {
627 [ + - ]: 983248 : ngbvs.insert( nxt );
628 [ + - ]: 983248 : todo.push_back( nxt );
629 : : }
630 : 5716382 : break;
631 : : }
632 : : }
633 : 5716382 : }
634 : 1091632 : }
635 [ + - ][ - + ]: 193170 : if( _MAXPNTS <= (int)ngbvs.size() )
636 : : {
637 : : // obtain enough points
638 : 0 : return error;
639 : : }
640 [ + + ]: 193170 : if( !todo.size() )
641 : : {
642 : : // current ring cannot introduce any points, return incase deadlock
643 : 159 : return error;
644 : : }
645 [ + + ][ + - ]: 193011 : if( ( i == ring ) && ( minpnts > (int)ngbvs.size() ) )
[ + + ][ + + ]
646 : : {
647 : : // reach maximum ring but not enough points
648 : 209 : ++ring;
649 : : }
650 : : }
651 : 67790 : return error;
652 : : }
653 : :
654 : 159 : void HiReconstruction::initialize_surf_geom( const int degree )
655 : : {
656 [ - + ]: 159 : if( !_hasderiv )
657 : : {
658 : 0 : compute_average_vertex_normals_surf();
659 : 0 : _hasderiv = true;
660 : : }
661 [ + - ]: 159 : if( !_initfittings )
662 : : {
663 : 159 : int ncoeffspv = ( degree + 2 ) * ( degree + 1 ) / 2;
664 [ + - ]: 159 : _degrees_out.assign( _nv2rec, 0 );
665 [ + - ]: 159 : _interps.assign( _nv2rec, false );
666 : 159 : _vertID2coeffID.resize( _nv2rec );
667 [ + - ]: 159 : _local_fit_coeffs.assign( _nv2rec * ncoeffspv, 0 );
668 [ + + ]: 67505 : for( size_t i = 0; i < _nv2rec; ++i )
669 : : {
670 : 67346 : _vertID2coeffID[i] = i * ncoeffspv;
671 : : }
672 : 159 : _initfittings = true;
673 : : }
674 : 159 : }
675 : :
676 : 0 : void HiReconstruction::initialize_surf_geom( const size_t npts, const int* degrees )
677 : : {
678 [ # # ]: 0 : if( !_hasderiv )
679 : : {
680 : 0 : compute_average_vertex_normals_surf();
681 : 0 : _hasderiv = true;
682 : : }
683 [ # # ]: 0 : if( !_initfittings )
684 : : {
685 [ # # ]: 0 : assert( _nv2rec == npts );
686 [ # # ]: 0 : _degrees_out.assign( _nv2rec, 0 );
687 [ # # ]: 0 : _interps.assign( _nv2rec, false );
688 : 0 : _vertID2coeffID.resize( _nv2rec );
689 : 0 : size_t index = 0;
690 [ # # ]: 0 : for( size_t i = 0; i < _nv2rec; ++i )
691 : : {
692 : 0 : _vertID2coeffID[i] = index;
693 : 0 : index += ( degrees[i] + 2 ) * ( degrees[i] + 1 ) / 2;
694 : : }
695 [ # # ]: 0 : _local_fit_coeffs.assign( index, 0 );
696 : 0 : _initfittings = true;
697 : : }
698 : 0 : }
699 : :
700 : 48 : void HiReconstruction::initialize_3Dcurve_geom( const int degree )
701 : : {
702 [ - + ]: 48 : if( !_hasderiv )
703 : : {
704 : 0 : compute_average_vertex_tangents_curve();
705 : 0 : _hasderiv = true;
706 : : }
707 [ + - ]: 48 : if( !_initfittings )
708 : : {
709 : 48 : int ncoeffspvpd = degree + 1;
710 [ + - ]: 48 : _degrees_out.assign( _nv2rec, 0 );
711 [ + - ]: 48 : _interps.assign( _nv2rec, false );
712 : 48 : _vertID2coeffID.resize( _nv2rec );
713 [ + - ]: 48 : _local_fit_coeffs.assign( _nv2rec * ncoeffspvpd * 3, 0 );
714 [ + + ]: 492 : for( size_t i = 0; i < _nv2rec; ++i )
715 : : {
716 : 444 : _vertID2coeffID[i] = i * ncoeffspvpd * 3;
717 : : }
718 : 48 : _initfittings = true;
719 : : }
720 : 48 : }
721 : :
722 : 0 : void HiReconstruction::initialize_3Dcurve_geom( const size_t npts, const int* degrees )
723 : : {
724 [ # # ]: 0 : if( !_hasderiv )
725 : : {
726 : 0 : compute_average_vertex_tangents_curve();
727 : 0 : _hasderiv = true;
728 : : }
729 [ # # ]: 0 : if( !_hasfittings )
730 : : {
731 [ # # ]: 0 : assert( _nv2rec == npts );
732 [ # # ]: 0 : _degrees_out.assign( _nv2rec, 0 );
733 [ # # ]: 0 : _interps.assign( _nv2rec, false );
734 : 0 : _vertID2coeffID.reserve( _nv2rec );
735 : 0 : size_t index = 0;
736 [ # # ]: 0 : for( size_t i = 0; i < _nv2rec; ++i )
737 : : {
738 : 0 : _vertID2coeffID[i] = index;
739 : 0 : index += 3 * ( degrees[i] + 1 );
740 : : }
741 [ # # ]: 0 : _local_fit_coeffs.assign( index, 0 );
742 : 0 : _initfittings = true;
743 : : }
744 : 0 : }
745 : :
746 : : /* ErrorCode HiReconstruction::set_geom_data_surf(const EntityHandle vid, const double* coords,
747 : : const double degree_out, const double* coeffs, bool interp)
748 : : {
749 : : return MB_SUCCESS;
750 : : }
751 : :
752 : : ErrorCode HiReconstruction::set_geom_data_3Dcurve(const EntityHandle vid, const double* coords,
753 : : const double degree_out, const double* coeffs, bool interp)
754 : : {
755 : : return MB_SUCCESS;
756 : : } */
757 : :
758 : : /*********************************************************
759 : : * Routines for vertex normal/tangent vector estimation *
760 : : *********************************************************/
761 : 7879 : ErrorCode HiReconstruction::average_vertex_normal( const EntityHandle vid, double* nrm )
762 : : {
763 : : ErrorCode error;
764 [ + - ]: 7879 : std::vector< EntityHandle > adjfaces;
765 [ + - ][ - + ]: 7879 : error = vertex_get_incident_elements( vid, 2, adjfaces );MB_CHK_ERR( error );
[ # # ][ # # ]
766 : 7879 : int npolys = adjfaces.size();
767 [ - + ][ # # ]: 7879 : if( !npolys ) { MB_SET_ERR( MB_FAILURE, "Vertex has no incident 2D entities" ); }
[ # # ][ # # ]
[ # # ][ # # ]
768 : : else
769 : : {
770 : : double v1[3], v2[3], v3[3], a[3], b[3], c[3];
771 : 7879 : nrm[0] = nrm[1] = nrm[2] = 0;
772 [ + + ]: 49269 : for( int i = 0; i < npolys; ++i )
773 : : {
774 : : // get incident "triangles"
775 [ + - ]: 41390 : std::vector< EntityHandle > elemconn;
776 [ + - ][ + - ]: 41390 : error = mbImpl->get_connectivity( &adjfaces[i], 1, elemconn );MB_CHK_ERR( error );
[ - + ][ # # ]
[ # # ]
777 : : EntityHandle pre, nxt;
778 : 41390 : int nvpe = elemconn.size();
779 [ + - ]: 88160 : for( int j = 0; j < nvpe; ++j )
780 : : {
781 [ + - ][ + + ]: 88160 : if( vid == elemconn[j] )
782 : : {
783 [ + + ][ + - ]: 41390 : pre = j == 0 ? elemconn[nvpe - 1] : elemconn[j - 1];
[ + - ]
784 [ + - ]: 41390 : nxt = elemconn[( j + 1 ) % nvpe];
785 : 41390 : break;
786 : : }
787 : : }
788 : : // compute area weighted normals
789 [ + - ][ - + ]: 41390 : error = mbImpl->get_coords( &pre, 1, a );MB_CHK_ERR( error );
[ # # ][ # # ]
790 [ + - ][ - + ]: 41390 : error = mbImpl->get_coords( &vid, 1, b );MB_CHK_ERR( error );
[ # # ][ # # ]
791 [ + - ][ - + ]: 41390 : error = mbImpl->get_coords( &nxt, 1, c );MB_CHK_ERR( error );
[ # # ][ # # ]
792 [ + - ]: 41390 : DGMSolver::vec_linear_operation( 3, 1, c, -1, b, v1 );
793 [ + - ]: 41390 : DGMSolver::vec_linear_operation( 3, 1, a, -1, b, v2 );
794 [ + - ]: 41390 : DGMSolver::vec_crossprod( v1, v2, v3 );
795 [ + - ][ + - ]: 41390 : DGMSolver::vec_linear_operation( 3, 1, nrm, 1, v3, nrm );
796 : 41390 : }
797 : : #ifndef NDEBUG
798 [ + - ][ - + ]: 7879 : assert( DGMSolver::vec_normalize( 3, nrm, nrm ) );
799 : : #endif
800 : : }
801 : 7879 : return error;
802 : : }
803 : :
804 : 17 : ErrorCode HiReconstruction::compute_average_vertex_normals_surf()
805 : : {
806 [ - + ]: 17 : if( _hasderiv ) { return MB_SUCCESS; }
807 : : ErrorCode error;
808 [ + - ]: 17 : _local_coords.assign( 9 * _nv2rec, 0 );
809 : 17 : size_t index = 0;
810 [ + - ][ + - ]: 7896 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
[ + - ][ + - ]
[ + + ]
811 : : {
812 [ + - ][ + - ]: 7879 : error = average_vertex_normal( *ivert, &( _local_coords[9 * index + 6] ) );MB_CHK_ERR( error );
[ + - ][ - + ]
[ # # ][ # # ]
813 : : }
814 : 17 : return error;
815 : : }
816 : :
817 : 67346 : ErrorCode HiReconstruction::get_normals_surf( const Range& vertsh, double* nrms )
818 : : {
819 : 67346 : ErrorCode error = MB_SUCCESS;
820 [ + - ]: 67346 : if( _hasderiv )
821 : : {
822 : 67346 : size_t id = 0;
823 [ + - ][ + - ]: 2101226 : for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
[ + - ][ + - ]
[ + + ]
824 : : {
825 [ + - ][ + - ]: 2033880 : int index = _verts2rec.index( *ivert );
826 : : #ifdef MOAB_HAVE_MPI
827 [ - + ]: 2033880 : if( -1 == index )
828 : : {
829 : : // ghost vertex
830 [ # # ][ # # ]: 0 : error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
[ # # ][ # # ]
[ # # ]
831 : : }
832 : : else
833 : : {
834 [ + - ]: 2033880 : nrms[3 * id] = _local_coords[9 * index + 6];
835 [ + - ]: 2033880 : nrms[3 * id + 1] = _local_coords[9 * index + 7];
836 [ + - ]: 2033880 : nrms[3 * id + 2] = _local_coords[9 * index + 8];
837 : : }
838 : : #else
839 : : assert( -1 != index );
840 : : nrms[3 * id] = _local_coords[9 * index + 6];
841 : : nrms[3 * id + 1] = _local_coords[9 * index + 7];
842 : : nrms[3 * id + 2] = _local_coords[9 * index + 8];
843 : : #endif
844 : : }
845 : : }
846 : : else
847 : : {
848 : 0 : size_t id = 0;
849 [ # # ][ # # ]: 0 : for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
[ # # ][ # # ]
[ # # ]
850 : : {
851 [ # # ][ # # ]: 0 : error = average_vertex_normal( *ivert, nrms + 3 * id );MB_CHK_ERR( error );
[ # # ][ # # ]
[ # # ]
852 : : }
853 : : }
854 : 67346 : return error;
855 : : }
856 : :
857 : 37 : ErrorCode HiReconstruction::average_vertex_tangent( const EntityHandle vid, double* tang )
858 : : {
859 : : ErrorCode error;
860 [ + - ]: 37 : std::vector< EntityHandle > adjedges;
861 [ + - ][ - + ]: 37 : error = vertex_get_incident_elements( vid, 1, adjedges );MB_CHK_ERR( error );
[ # # ][ # # ]
862 : 37 : int nedges = adjedges.size();
863 [ - + ][ # # ]: 37 : if( !nedges ) { MB_SET_ERR( MB_FAILURE, "Vertex has no incident edges" ); }
[ # # ][ # # ]
[ # # ][ # # ]
864 : : else
865 : : {
866 [ - + ]: 37 : assert( nedges <= 2 );
867 : 37 : tang[0] = tang[1] = tang[2] = 0;
868 [ + + ]: 111 : for( int i = 0; i < nedges; ++i )
869 : : {
870 [ + - ]: 74 : std::vector< EntityHandle > edgeconn;
871 [ + - ][ + - ]: 74 : error = mbImpl->get_connectivity( &adjedges[i], 1, edgeconn );
872 : : double istr[3], iend[3], t[3];
873 [ + - ][ + - ]: 74 : error = mbImpl->get_coords( &( edgeconn[0] ), 1, istr );
874 [ + - ][ + - ]: 74 : error = mbImpl->get_coords( &( edgeconn[1] ), 1, iend );
875 [ + - ]: 74 : DGMSolver::vec_linear_operation( 3, 1, iend, -1, istr, t );
876 [ + - ]: 74 : DGMSolver::vec_linear_operation( 3, 1, tang, 1, t, tang );
877 : 74 : }
878 : : #ifndef NDEBUG
879 [ + - ][ - + ]: 37 : assert( DGMSolver::vec_normalize( 3, tang, tang ) );
880 : : #endif
881 : : }
882 : 37 : return error;
883 : : }
884 : :
885 : 4 : ErrorCode HiReconstruction::compute_average_vertex_tangents_curve()
886 : : {
887 [ - + ]: 4 : if( _hasderiv ) { return MB_SUCCESS; }
888 : : ErrorCode error;
889 [ + - ]: 4 : _local_coords.assign( 3 * _nv2rec, 0 );
890 : 4 : size_t index = 0;
891 [ + - ][ + - ]: 41 : for( Range::iterator ivert = _verts2rec.begin(); ivert != _verts2rec.end(); ++ivert, ++index )
[ + - ][ + - ]
[ + + ]
892 : : {
893 [ + - ][ + - ]: 37 : error = average_vertex_tangent( *ivert, &( _local_coords[3 * index] ) );MB_CHK_ERR( error );
[ + - ][ - + ]
[ # # ][ # # ]
894 : : }
895 : 4 : return error;
896 : : }
897 : :
898 : 444 : ErrorCode HiReconstruction::get_tangents_curve( const Range& vertsh, double* tangs )
899 : : {
900 : 444 : ErrorCode error = MB_SUCCESS;
901 [ + - ]: 444 : if( _hasderiv )
902 : : {
903 : 444 : size_t id = 0;
904 [ + - ][ + - ]: 3144 : for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
[ + - ][ + - ]
[ + + ]
905 : : {
906 [ + - ][ + - ]: 2700 : int index = _verts2rec.index( *ivert );
907 : : #ifdef MOAB_HAVE_MPI
908 [ + - ]: 2700 : if( -1 != index )
909 : : {
910 [ + - ]: 2700 : tangs[3 * id] = _local_coords[3 * index];
911 [ + - ]: 2700 : tangs[3 * id + 1] = _local_coords[3 * index + 1];
912 [ + - ]: 2700 : tangs[3 * id + 2] = _local_coords[3 * index + 2];
913 : : }
914 : : else
915 : : {
916 [ # # ][ # # ]: 0 : error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
[ # # ][ # # ]
[ # # ]
917 : : }
918 : : #else
919 : : assert( -1 != index );
920 : : tangs[3 * id] = _local_coords[3 * index];
921 : : tangs[3 * id + 1] = _local_coords[3 * index + 1];
922 : : tangs[3 * id + 2] = _local_coords[3 * index + 2];
923 : : #endif
924 : : }
925 : : }
926 : : else
927 : : {
928 : 0 : size_t id = 0;
929 [ # # ][ # # ]: 0 : for( Range::iterator ivert = vertsh.begin(); ivert != vertsh.end(); ++ivert, ++id )
[ # # ][ # # ]
[ # # ]
930 : : {
931 [ # # ][ # # ]: 0 : error = average_vertex_tangent( *ivert, tangs + 3 * id );MB_CHK_ERR( error );
[ # # ][ # # ]
[ # # ]
932 : : }
933 : : }
934 : 444 : return error;
935 : : }
936 : :
937 : : /************************************************
938 : : * Internal Routines for local WLS fittings *
939 : : *************************************************/
940 : :
941 : 67346 : void HiReconstruction::polyfit3d_surf_get_coeff( const int nverts, const double* ngbcoords, const double* ngbnrms,
942 : : int degree, const bool interp, const bool safeguard, const int ncoords,
943 : : double* coords, const int ncoeffs, double* coeffs, int* degree_out,
944 : : int* degree_pnt, int* degree_qr )
945 : : {
946 [ - + ]: 134692 : if( nverts <= 0 ) { return; }
947 : :
948 : : // std::cout << "npnts in initial stencil = " << nverts << std::endl;
949 : : // std::cout << "centered at (" << ngbcoords[0] << "," << ngbcoords[1] << "," << ngbcoords[2] <<
950 : : // ")" << std::endl;
951 : :
952 : : // step 1. copmute local coordinate system
953 : 67346 : double nrm[3] = { ngbnrms[0], ngbnrms[1], ngbnrms[2] }, tang1[3] = { 0, 0, 0 }, tang2[3] = { 0, 0, 0 };
954 [ + + ][ + + ]: 67346 : if( fabs( nrm[0] ) > fabs( nrm[1] ) && fabs( nrm[0] ) > fabs( nrm[2] ) ) { tang1[1] = 1.0; }
955 : : else
956 : : {
957 : 47079 : tang1[0] = 1.0;
958 : : }
959 : :
960 [ + - ]: 67346 : DGMSolver::vec_projoff( 3, tang1, nrm, tang1 );
961 : : #ifndef NDEBUG
962 [ + - ][ - + ]: 67346 : assert( DGMSolver::vec_normalize( 3, tang1, tang1 ) );
963 : : #endif
964 [ + - ]: 67346 : DGMSolver::vec_crossprod( nrm, tang1, tang2 );
965 [ + - ][ + - ]: 67346 : if( 9 <= ncoords && coords )
966 : : {
967 : 67346 : coords[0] = tang1[0];
968 : 67346 : coords[1] = tang1[1];
969 : 67346 : coords[2] = tang1[2];
970 : 67346 : coords[3] = tang2[0];
971 : 67346 : coords[4] = tang2[1];
972 : 67346 : coords[5] = tang2[2];
973 : 67346 : coords[6] = nrm[0];
974 : 67346 : coords[7] = nrm[1];
975 : 67346 : coords[8] = nrm[2];
976 : : }
977 [ + - ][ - + ]: 67346 : if( !ncoeffs || !coeffs ) { return; }
978 : : else
979 : : {
980 [ - + ]: 67346 : assert( ncoeffs >= ( degree + 2 ) * ( degree + 1 ) / 2 );
981 : : }
982 : :
983 : : // step 2. project onto local coordinates system
984 : 67346 : int npts2fit = nverts - interp;
985 [ - + ]: 67346 : if( 0 == npts2fit )
986 : : {
987 : 0 : *degree_out = *degree_pnt = *degree_qr = 0;
988 [ # # ]: 0 : for( int i = 0; i < ncoeffs; ++i )
989 : : {
990 : 0 : coeffs[i] = 0;
991 : : }
992 : : // coeffs[0] = 0;
993 : 0 : return;
994 : : }
995 [ + - ]: 67346 : std::vector< double > us( npts2fit * 2 ); // double *us = new double[npts2fit*2];
996 [ + - ][ + - ]: 134692 : std::vector< double > bs( npts2fit ); // double *bs = new double[npts2fit];
997 [ + + ]: 2080864 : for( int i = interp; i < nverts; ++i )
998 : : {
999 : 2013518 : int k = i - interp;
1000 : : double uu[3];
1001 [ + - ]: 2013518 : DGMSolver::vec_linear_operation( 3, 1, ngbcoords + 3 * i, -1, ngbcoords, uu );
1002 [ + - ][ + - ]: 2013518 : us[k * 2] = DGMSolver::vec_innerprod( 3, tang1, uu );
1003 [ + - ][ + - ]: 2013518 : us[k * 2 + 1] = DGMSolver::vec_innerprod( 3, tang2, uu );
1004 [ + - ][ + - ]: 2013518 : bs[k] = DGMSolver::vec_innerprod( 3, nrm, uu );
1005 : : }
1006 : :
1007 : : // step 3. compute weights
1008 [ + - ][ + - ]: 134692 : std::vector< double > ws( npts2fit ); // double *ws = new double[npts2fit];
1009 [ + - ][ + - ]: 67346 : int nzeros = compute_weights( npts2fit, 2, &( us[0] ), nverts, ngbnrms, degree, _MINEPS, &( ws[0] ) );
[ + - ]
1010 : :
1011 : : // step 4. adjust according to zero-weights
1012 [ + + ]: 67346 : if( nzeros )
1013 : : {
1014 [ - + ]: 562 : if( nzeros == npts2fit )
1015 : : {
1016 : 0 : *degree_out = *degree_pnt = *degree_qr = 0;
1017 [ # # ]: 0 : for( int i = 0; i < ncoeffs; ++i )
1018 : : {
1019 : 0 : coeffs[i] = 0;
1020 : : }
1021 : : // coeffs[0] = 0;
1022 : 0 : return;
1023 : : }
1024 : 562 : int index = 0;
1025 [ + + ]: 21188 : for( int i = 0; i < npts2fit; ++i )
1026 : : {
1027 [ + - ][ + + ]: 20626 : if( ws[i] )
1028 : : {
1029 [ + + ]: 15190 : if( i > index )
1030 : : {
1031 [ + - ][ + - ]: 12625 : us[index * 2] = us[i * 2];
1032 [ + - ][ + - ]: 12625 : us[index * 2 + 1] = us[i * 2 + 1];
1033 [ + - ][ + - ]: 12625 : bs[index] = bs[i];
1034 [ + - ][ + - ]: 12625 : ws[index] = ws[i];
1035 : : }
1036 : 15190 : ++index;
1037 : : }
1038 : : }
1039 : 562 : npts2fit -= nzeros;
1040 [ - + ]: 562 : assert( index == npts2fit );
1041 [ + - ]: 562 : us.resize( npts2fit * 2 );
1042 [ + - ]: 562 : bs.resize( npts2fit );
1043 [ + - ]: 562 : ws.resize( npts2fit );
1044 : : /*us = realloc(us,npts2fit*2*sizeof(double));
1045 : : bs = realloc(bs,npts2fit*sizeof(double));
1046 : : ws = realloc(ws,npts2fit*sizeof(double));*/
1047 : : }
1048 : :
1049 : : // std::cout<<"npnts after weighting = "<<npts2fit<<std::endl;
1050 : :
1051 : : // step 5. fitting
1052 [ + - ][ + - ]: 67346 : eval_vander_bivar_cmf( npts2fit, &( us[0] ), 1, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out,
[ + - ]
1053 [ + - ]: 67346 : degree_pnt, degree_qr );
1054 : :
1055 : : // step 6. organize output
1056 : 67346 : int ncoeffs_out = ( *degree_out + 2 ) * ( *degree_out + 1 ) / 2;
1057 [ - + ]: 67346 : assert( ncoeffs_out <= ncoeffs );
1058 : 67346 : coeffs[0] = 0;
1059 [ + + ][ + - ]: 968002 : for( int j = 0; j < ncoeffs_out - interp; ++j )
1060 : : {
1061 [ + - ]: 900656 : coeffs[j + interp] = bs[j];
1062 : 67346 : }
1063 : : // delete [] us; delete [] bs; delete [] ws;
1064 : : }
1065 : :
1066 : 67346 : void HiReconstruction::eval_vander_bivar_cmf( const int npts2fit, const double* us, const int ndim, double* bs,
1067 : : int degree, const double* ws, const bool interp, const bool safeguard,
1068 : : int* degree_out, int* degree_pnt, int* degree_qr )
1069 : : {
1070 : : // step 1. adjust the degree according to number of points to fit
1071 : 67346 : int ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1072 [ + + ][ + + ]: 68662 : while( 1 < degree && npts2fit < ncols )
1073 : : {
1074 : 1316 : --degree;
1075 : 1316 : ncols = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1076 : : }
1077 : 67346 : *degree_pnt = degree;
1078 : :
1079 : : // std::cout << "degree_pnt: " << *degree_pnt << std::endl;
1080 : :
1081 : : // step 2. construct Vandermonde matrix, stored in columnwise
1082 [ + - ]: 67346 : std::vector< double > V; // V(npts2fit*(ncols+interp)); //double *V_init = new double[npts2fit*(ncols+interp)];
1083 [ + - ]: 67346 : DGMSolver::gen_vander_multivar( npts2fit, 2, us, degree, V );
1084 : : // remove the first column of 1s if interpolation
1085 [ + + ][ + - ]: 67346 : if( interp ) { V.erase( V.begin(), V.begin() + npts2fit ); }
[ + - ]
1086 : : /*double* V;
1087 : : if(interp){
1088 : : V = new double[npts2fit*ncols];
1089 : : std::memcpy(V,V_init+npts2fit,ncols*npts2fit*sizeof(double));
1090 : : delete [] V_init; V_init = 0;
1091 : : }else{
1092 : : V = V_init;
1093 : : }*/
1094 : :
1095 : : // step 3. Scale rows to assign different weights to different points
1096 [ + + ]: 2075428 : for( int i = 0; i < npts2fit; ++i )
1097 : : {
1098 [ + + ]: 36978053 : for( int j = 0; j < ncols; ++j )
1099 : : {
1100 [ + - ]: 34969971 : V[j * npts2fit + i] *= ws[i];
1101 : : }
1102 [ + + ]: 4016164 : for( int k = 0; k < ndim; ++k )
1103 : : {
1104 : 2008082 : bs[k * npts2fit + i] *= ws[i];
1105 : : }
1106 : : }
1107 : :
1108 : : // step 4. scale columns to reduce condition number
1109 [ + - ][ + - ]: 134692 : std::vector< double > ts( ncols ); // double *ts = new double[ncols];
1110 [ + - ][ + - ]: 67346 : DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
[ + - ]
1111 : :
1112 : : // step 5. Perform Householder QR factorization
1113 [ + - ][ + - ]: 134692 : std::vector< double > D( ncols ); // double *D = new double[ncols];
1114 : : int rank;
1115 [ + - ][ + - ]: 67346 : DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
[ + - ]
1116 : :
1117 : : // step 6. adjust degree of fitting according to rank of Vandermonde matrix
1118 : 67346 : int ncols_sub = ncols;
1119 [ + + ]: 67810 : while( rank < ncols_sub )
1120 : : {
1121 : 464 : --degree;
1122 [ - + ]: 464 : if( degree == 0 )
1123 : : {
1124 : : // surface is flat, return 0
1125 : 0 : *degree_out = *degree_qr = degree;
1126 [ # # ]: 0 : for( int i = 0; i < npts2fit; ++i )
1127 : : {
1128 [ # # ]: 0 : for( int k = 0; k < ndim; ++k )
1129 : : {
1130 : 0 : bs[k * npts2fit + i] = 0;
1131 : : }
1132 : : }
1133 : 67346 : return;
1134 : : }
1135 : : else
1136 : : {
1137 : 464 : ncols_sub = ( ( ( degree + 2 ) * ( degree + 1 ) ) >> 1 ) - interp;
1138 : : }
1139 : : }
1140 : 67346 : *degree_qr = degree;
1141 : :
1142 : : // std::cout << "degree_qr: " << *degree_qr << std::endl;
1143 : :
1144 : : /* DBG
1145 : : * std::cout<<"before Qtb"<<std::endl;
1146 : : std::cout<<std::endl;
1147 : : std::cout<<"bs = "<<std::endl;
1148 : : std::cout<<std::endl;
1149 : : for (int k=0; k< ndim; k++){
1150 : : for (int j=0; j<npts2fit; ++j){
1151 : : std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1152 : : }
1153 : : }
1154 : : std::cout<<std::endl;*/
1155 : :
1156 : : // step 7. compute Q'b
1157 [ + - ][ + - ]: 67346 : DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1158 : :
1159 : : /* DBG
1160 : : * std::cout<<"after Qtb"<<std::endl;
1161 : : std::cout<<"bs = "<<std::endl;
1162 : : std::cout<<std::endl;
1163 : : for (int k=0; k< ndim; k++){
1164 : : for (int j=0; j<npts2fit; ++j){
1165 : : std::cout<<" "<<bs[npts2fit*k+j]<<std::endl;
1166 : : }
1167 : : }
1168 : : std::cout<<std::endl;*/
1169 : :
1170 : : // step 8. perform backward substitution and scale the solution
1171 : : // assign diagonals of V
1172 [ + + ]: 968002 : for( int i = 0; i < ncols_sub; ++i )
1173 : : {
1174 [ + - ][ + - ]: 900656 : V[i * npts2fit + i] = D[i];
1175 : : }
1176 : :
1177 : : // backsolve
1178 [ + + ]: 67346 : if( safeguard )
1179 : : {
1180 : : // for debug
1181 : : // std::cout << "ts size " << ts.size() << std::endl;
1182 [ + - ]: 4 : DGMSolver::backsolve_polyfit_safeguarded( 2, degree, interp, npts2fit, ncols_sub, &( V[0] ), ndim, bs,
1183 [ + - ][ + - ]: 8 : &( ts[0] ), degree_out );
1184 : : }
1185 : : else
1186 : : {
1187 [ + - ][ + - ]: 67342 : DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), 1, bs, &( ts[0] ) );
[ + - ]
1188 [ + - ]: 67346 : *degree_out = degree;
1189 : 67346 : }
1190 : : /*if(V_init){
1191 : : delete [] V_init;
1192 : : }else{
1193 : : delete [] V;
1194 : : }*/
1195 : : }
1196 : :
1197 : 444 : void HiReconstruction::polyfit3d_curve_get_coeff( const int nverts, const double* ngbcors, const double* ngbtangs,
1198 : : int degree, const bool interp, const bool safeguard,
1199 : : const int ncoords, double* coords, const int ncoeffs, double* coeffs,
1200 : : int* degree_out )
1201 : : {
1202 [ - + ]: 888 : if( !nverts ) { return; }
1203 : : // step 1. compute local coordinates system
1204 : 444 : double tang[3] = { ngbtangs[0], ngbtangs[1], ngbtangs[2] };
1205 [ - + ][ # # ]: 444 : if( coords && ncoords > 2 )
1206 : : {
1207 : 0 : coords[0] = tang[0];
1208 : 0 : coords[1] = tang[1];
1209 : 0 : coords[2] = tang[2];
1210 : : }
1211 [ + - ][ - + ]: 444 : if( !coeffs || !ncoeffs ) { return; }
1212 : : else
1213 : : {
1214 [ - + ]: 444 : assert( ncoeffs >= 3 * ( degree + 1 ) );
1215 : : }
1216 : : // step 2. project onto local coordinate system
1217 : 444 : int npts2fit = nverts - interp;
1218 [ - + ]: 444 : if( !npts2fit )
1219 : : {
1220 : 0 : *degree_out = 0;
1221 [ # # ]: 0 : for( int i = 0; i < ncoeffs; ++i )
1222 : : {
1223 : 0 : coeffs[0] = 0;
1224 : : }
1225 : 0 : return;
1226 : : }
1227 [ + - ]: 444 : std::vector< double > us( npts2fit ); // double *us = new double[npts2fit];
1228 [ + - ][ + + ]: 888 : std::vector< double > bs( npts2fit * 3 ); // double *bs = new double[npts2fit*3];
1229 : : double uu[3];
1230 [ + + ]: 2922 : for( int i = interp; i < nverts; ++i )
1231 : : {
1232 : 2478 : int k = i - interp;
1233 [ + - ]: 2478 : DGMSolver::vec_linear_operation( 3, 1, ngbcors + 3 * i, -1, ngbcors, uu );
1234 [ + - ][ + - ]: 2478 : us[k] = DGMSolver::vec_innerprod( 3, uu, tang );
1235 [ + - ]: 2478 : bs[k] = uu[0];
1236 [ + - ]: 2478 : bs[npts2fit + k] = uu[1];
1237 [ + - ]: 2478 : bs[2 * npts2fit + k] = uu[2];
1238 : : }
1239 : :
1240 : : // step 3. copmute weights
1241 [ + - ][ + + ]: 888 : std::vector< double > ws( npts2fit );
1242 [ + - ][ + - ]: 444 : int nzeros = compute_weights( npts2fit, 1, &( us[0] ), nverts, ngbtangs, degree, _MINEPS, &( ws[0] ) );
[ + - ]
1243 [ - + ]: 444 : assert( nzeros <= npts2fit );
1244 : :
1245 : : // step 4. adjust according to number of zero-weights
1246 [ + + ]: 444 : if( nzeros )
1247 : : {
1248 [ + + ]: 154 : if( nzeros == npts2fit )
1249 : : {
1250 : : // singular case
1251 : 42 : *degree_out = 0;
1252 [ + + ]: 609 : for( int i = 0; i < ncoeffs; ++i )
1253 : : {
1254 : 567 : coeffs[i] = 0;
1255 : : }
1256 : 42 : return;
1257 : : }
1258 : 112 : int npts_new = npts2fit - nzeros;
1259 [ + - ]: 112 : std::vector< double > bs_temp( 3 * npts_new );
1260 : 112 : int index = 0;
1261 [ + + ]: 782 : for( int i = 0; i < npts2fit; ++i )
1262 : : {
1263 [ + - ][ + + ]: 670 : if( ws[i] )
1264 : : {
1265 [ + + ]: 362 : if( i > index )
1266 : : {
1267 [ + - ][ + - ]: 172 : us[index] = us[i];
1268 [ + - ][ + - ]: 172 : ws[index] = ws[i];
1269 : : }
1270 [ + - ][ + - ]: 362 : bs_temp[index] = bs[i];
1271 [ + - ][ + - ]: 362 : bs_temp[index + npts_new] = bs[i + npts2fit];
1272 [ + - ][ + - ]: 362 : bs_temp[index + 2 * npts_new] = bs[i + 2 * npts2fit];
1273 : 362 : ++index;
1274 : : }
1275 : : }
1276 [ - + ]: 112 : assert( index == npts_new );
1277 : 112 : npts2fit = npts_new;
1278 [ + - ]: 112 : us.resize( index );
1279 [ + - ]: 112 : ws.resize( index );
1280 [ + - ]: 112 : bs = bs_temp;
1281 : : // destroy bs_temp;
1282 [ + - ]: 112 : std::vector< double >().swap( bs_temp );
1283 : : }
1284 : :
1285 : : // step 5. fitting
1286 [ + - ][ + - ]: 402 : eval_vander_univar_cmf( npts2fit, &( us[0] ), 3, &( bs[0] ), degree, &( ws[0] ), interp, safeguard, degree_out );
[ + - ][ + - ]
1287 : : // step 6. write results to output pointers
1288 : 402 : int ncoeffs_out_pvpd = *degree_out + 1;
1289 [ - + ]: 402 : assert( ncoeffs >= 3 * ncoeffs_out_pvpd );
1290 : : // write to coeffs consecutively, bs's size is not changed by eval_vander_univar_cmf
1291 : 402 : coeffs[0] = coeffs[ncoeffs_out_pvpd] = coeffs[2 * ncoeffs_out_pvpd] = 0;
1292 [ + + ][ + + ]: 1866 : for( int i = 0; i < ncoeffs_out_pvpd - interp; ++i )
1293 : : {
1294 [ + - ]: 1422 : coeffs[i + interp] = bs[i];
1295 [ + - ]: 1422 : coeffs[i + interp + ncoeffs_out_pvpd] = bs[i + npts2fit];
1296 [ + - ]: 1422 : coeffs[i + interp + 2 * ncoeffs_out_pvpd] = bs[i + 2 * npts2fit];
1297 : 402 : }
1298 : : }
1299 : :
1300 : 402 : void HiReconstruction::eval_vander_univar_cmf( const int npts2fit, const double* us, const int ndim, double* bs,
1301 : : int degree, const double* ws, const bool interp, const bool safeguard,
1302 : : int* degree_out )
1303 : : {
1304 : : // step 1. determine degree of polynomials to fit according to number of points
1305 : 402 : int ncols = degree + 1 - interp;
1306 [ + + ][ + - ]: 609 : while( npts2fit < ncols && degree >= 1 )
1307 : : {
1308 : 207 : --degree;
1309 : 207 : ncols = degree + 1 - interp;
1310 : : }
1311 [ + + ]: 402 : if( !degree )
1312 : : {
1313 [ - + ]: 42 : if( interp )
1314 : : {
1315 [ # # ]: 0 : for( int icol = 0; icol < ndim; ++icol )
1316 : : {
1317 : 0 : bs[icol * npts2fit] = 0;
1318 : : }
1319 : : }
1320 [ - + ]: 42 : for( int irow = 1; irow < npts2fit; ++irow )
1321 : : {
1322 [ # # ]: 0 : for( int icol = 0; icol < ndim; ++icol )
1323 : : {
1324 : 0 : bs[icol * npts2fit + irow] = 0;
1325 : : }
1326 : : }
1327 : 42 : *degree_out = 0;
1328 : 402 : return;
1329 : : }
1330 : : // step 2. construct Vandermonde matrix
1331 [ + - ]: 360 : std::vector< double > V; // V(npts2fit*(ncols+interp));
1332 [ + - ]: 360 : DGMSolver::gen_vander_multivar( npts2fit, 1, us, degree, V );
1333 : :
1334 [ + + ][ + - ]: 360 : if( interp ) { V.erase( V.begin(), V.begin() + npts2fit ); }
[ + - ]
1335 : :
1336 : : // step 3. scale rows with respect to weights
1337 [ + + ]: 2380 : for( int i = 0; i < npts2fit; ++i )
1338 : : {
1339 [ + + ]: 10560 : for( int j = 0; j < ncols; ++j )
1340 : : {
1341 [ + - ]: 8540 : V[j * npts2fit + i] *= ws[i];
1342 : : }
1343 [ + + ]: 8080 : for( int k = 0; k < ndim; ++k )
1344 : : {
1345 : 6060 : bs[k * npts2fit + i] *= ws[i];
1346 : : }
1347 : : }
1348 : :
1349 : : // step 4. scale columns to reduce condition number
1350 [ + - ]: 720 : std::vector< double > ts( ncols );
1351 [ + - ][ + - ]: 360 : DGMSolver::rescale_matrix( npts2fit, ncols, &( V[0] ), &( ts[0] ) );
[ + - ]
1352 : :
1353 : : // step 5. perform Householder QR factorization
1354 [ + - ]: 720 : std::vector< double > D( ncols );
1355 : : int rank;
1356 [ + - ][ + - ]: 360 : DGMSolver::qr_polyfit_safeguarded( npts2fit, ncols, &( V[0] ), &( D[0] ), &rank );
[ + - ]
1357 : :
1358 : : // step 6. adjust degree of fitting
1359 : 360 : int ncols_sub = ncols;
1360 [ - + ]: 360 : while( rank < ncols_sub )
1361 : : {
1362 : 0 : --degree;
1363 [ # # ]: 0 : if( !degree )
1364 : : {
1365 : : // matrix is singular, consider curve on flat plane passing center and orthogonal to
1366 : : // tangent line
1367 : 0 : *degree_out = 0;
1368 [ # # ]: 0 : for( int i = 0; i < npts2fit; ++i )
1369 : : {
1370 [ # # ]: 0 : for( int k = 0; k < ndim; ++k )
1371 : : {
1372 : 0 : bs[k * npts2fit + i] = 0;
1373 : : }
1374 : : }
1375 : : }
1376 : 0 : ncols_sub = degree + 1 - interp;
1377 : : }
1378 : :
1379 : : // step 7. compute Q'*bs
1380 [ + - ][ + - ]: 360 : DGMSolver::compute_qtransposeB( npts2fit, ncols_sub, &( V[0] ), ndim, bs );
1381 : :
1382 : : // step 8. perform backward substitution and scale solutions
1383 : : // assign diagonals of V
1384 [ + + ]: 1740 : for( int i = 0; i < ncols_sub; ++i )
1385 : : {
1386 [ + - ][ + - ]: 1380 : V[i * npts2fit + i] = D[i];
1387 : : }
1388 : : // backsolve
1389 [ - + ]: 360 : if( safeguard )
1390 : : {
1391 [ # # ]: 0 : DGMSolver::backsolve_polyfit_safeguarded( 1, degree, interp, npts2fit, ncols, &( V[0] ), ndim, bs, ws,
1392 [ # # ]: 0 : degree_out );
1393 : : }
1394 : : else
1395 : : {
1396 [ + - ][ + - ]: 360 : DGMSolver::backsolve( npts2fit, ncols_sub, &( V[0] ), ndim, bs, &( ts[0] ) );
[ + - ]
1397 : 360 : *degree_out = degree;
1398 : 360 : }
1399 : : }
1400 : :
1401 : 67790 : int HiReconstruction::compute_weights( const int nrows, const int ncols, const double* us, const int nngbs,
1402 : : const double* ngbnrms, const int degree, const double toler, double* ws )
1403 : : {
1404 [ + - ][ - + ]: 67790 : assert( nrows <= _MAXPNTS && ws );
1405 : 67790 : bool interp = false;
1406 [ + + ]: 67790 : if( nngbs != nrows )
1407 : : {
1408 [ - + ]: 20584 : assert( nngbs == nrows + 1 );
1409 : 20584 : interp = true;
1410 : : }
1411 : 67790 : double epsilon = 1e-2;
1412 : :
1413 : : // First, compute squared distance from each input piont to the center
1414 [ + + ]: 2083786 : for( int i = 0; i < nrows; ++i )
1415 : : {
1416 : 2015996 : ws[i] = DGMSolver::vec_innerprod( ncols, us + i * ncols, us + i * ncols );
1417 : : }
1418 : :
1419 : : // Second, compute a small correction termt o guard against zero
1420 : 67790 : double h = 0;
1421 [ + + ]: 2083786 : for( int i = 0; i < nrows; ++i )
1422 : : {
1423 : 2015996 : h += ws[i];
1424 : : }
1425 : 67790 : h /= (double)nrows;
1426 : :
1427 : : // Finally, compute the weights for each vertex
1428 : 67790 : int nzeros = 0;
1429 [ + + ]: 2083786 : for( int i = 0; i < nrows; ++i )
1430 : : {
1431 : 2015996 : double costheta = DGMSolver::vec_innerprod( 3, ngbnrms, ngbnrms + 3 * ( i + interp ) );
1432 [ + + ]: 2015996 : if( costheta > toler ) { ws[i] = costheta * pow( ws[i] / h + epsilon, -1 * (double)degree / 2.0 ); }
1433 : : else
1434 : : {
1435 : 5852 : ws[i] = 0;
1436 : 5852 : ++nzeros;
1437 : : }
1438 : : }
1439 : 67790 : return nzeros;
1440 : : }
1441 : 581572 : bool HiReconstruction::check_barycentric_coords( const int nws, const double* naturalcoords )
1442 : : {
1443 : 581572 : double sum = 0;
1444 [ + + ]: 2358124 : for( int i = 0; i < nws; ++i )
1445 : : {
1446 [ - + ]: 1776552 : if( naturalcoords[i] < -_MINEPS ) { return false; }
1447 : 1776552 : sum += naturalcoords[i];
1448 : : }
1449 [ - + ]: 581572 : if( fabs( 1 - sum ) > _MINEPS ) { return false; }
1450 : : else
1451 : : {
1452 : 581572 : return true;
1453 : : }
1454 : : }
1455 [ + - ][ + - ]: 8 : } // namespace moab
|
