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 PatchData.cpp
29 : :
30 : : \author Thomas Leurent
31 : : \author Michael Brewer
32 : : \date 2002-01-17
33 : : */
34 : :
35 : : #include "PatchData.hpp"
36 : : #include "MsqMeshEntity.hpp"
37 : : #include "MsqFreeVertexIndexIterator.hpp"
38 : : #include "MeshInterface.hpp"
39 : : #include "MsqTimer.hpp"
40 : : #include "MsqDebug.hpp"
41 : : #include "GlobalPatch.hpp"
42 : : #include "PatchIterator.hpp"
43 : : #include "ExtraData.hpp"
44 : : #include "Settings.hpp"
45 : : #include "MappingFunction.hpp"
46 : :
47 : : #include <list>
48 : : #include <vector>
49 : : #include <map>
50 : : #include <algorithm>
51 : : #include <numeric>
52 : : #include <functional>
53 : : #include <utility>
54 : : #include <iostream>
55 : : #include <iomanip>
56 : : using std::endl;
57 : : using std::internal;
58 : : using std::left;
59 : : using std::list;
60 : : using std::map;
61 : : using std::ostream;
62 : : using std::setfill;
63 : : using std::setw;
64 : : using std::vector;
65 : :
66 : : namespace MBMesquite
67 : : {
68 : :
69 : 30 : const Settings PatchData::defaultSettings;
70 : :
71 : 1263 : PatchData::PatchData()
72 : : : myMesh( 0 ), myDomain( 0 ), numFreeVertices( 0 ), numSlaveVertices( 0 ), haveComputedInfos( 0 ), dataList( 0 ),
73 [ + - ][ + - ]: 1263 : mSettings( &defaultSettings )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
74 : : {
75 : 1263 : }
76 : :
77 : : // Destructor
78 : 2526 : PatchData::~PatchData()
79 : : {
80 : 1263 : notify_patch_destroyed();
81 : 1263 : }
82 : :
83 : 0 : void PatchData::get_minmax_element_unsigned_area( double& min, double& max, MsqError& err )
84 : : {
85 [ # # ][ # # ]: 0 : if( !have_computed_info( MAX_UNSIGNED_AREA ) || !have_computed_info( MIN_UNSIGNED_AREA ) )
[ # # ]
86 : : {
87 : 0 : max = 0;
88 : 0 : min = MSQ_DBL_MAX;
89 : 0 : size_t count = num_elements();
90 [ # # ]: 0 : for( size_t i = 0; i < count; ++i )
91 : : {
92 : : double vol;
93 [ # # ]: 0 : assert( i < elementArray.size() );
94 [ # # ][ # # ]: 0 : vol = elementArray[i].compute_unsigned_area( *this, err );MSQ_ERRRTN( err );
[ # # ]
95 [ # # ]: 0 : if( vol > max ) max = vol;
96 [ # # ]: 0 : if( vol < min ) min = vol;
97 : : }
98 : 0 : note_have_info( MAX_UNSIGNED_AREA );
99 : 0 : note_have_info( MIN_UNSIGNED_AREA );
100 : 0 : computedInfos[MAX_UNSIGNED_AREA] = max;
101 : 0 : computedInfos[MIN_UNSIGNED_AREA] = min;
102 : : }
103 : : else
104 : : {
105 : 0 : max = computedInfos[MAX_UNSIGNED_AREA];
106 : 0 : min = computedInfos[MIN_UNSIGNED_AREA];
107 : : }
108 : :
109 [ # # ][ # # ]: 0 : if( max <= 0 || min < 0 || min == MSQ_DBL_MAX ) MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
[ # # ][ # # ]
110 : : }
111 : :
112 : 58614 : void PatchData::get_minmax_edge_length( double& min, double& max ) const
113 : : {
114 : 58614 : min = HUGE_VAL;
115 : 58614 : max = -HUGE_VAL;
116 : :
117 [ + + ]: 296206 : for( size_t i = 0; i < num_elements(); ++i )
118 : : {
119 : 237592 : const MsqMeshEntity& elem = element_by_index( i );
120 : 237592 : const size_t* conn = elem.get_vertex_index_array();
121 : 237592 : const unsigned num_edges = TopologyInfo::edges( elem.get_element_type() );
122 [ + + ]: 1190877 : for( unsigned e = 0; e < num_edges; ++e )
123 : : {
124 : 953285 : const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), e );
125 : : const double len_sqr =
126 [ + - ]: 953285 : ( vertex_by_index( conn[edge[0]] ) - vertex_by_index( conn[edge[1]] ) ).length_squared();
127 [ + + ]: 953285 : if( len_sqr < min ) min = len_sqr;
128 [ + + ]: 953285 : if( len_sqr > max ) max = len_sqr;
129 : : }
130 : : }
131 : 58614 : min = sqrt( min );
132 : 58614 : max = sqrt( max );
133 : 58614 : }
134 : :
135 : : /*
136 : : double PatchData::get_average_Lambda_3d( MsqError &err)
137 : : {
138 : : double avg;
139 : : if (have_computed_info(AVERAGE_DET3D))
140 : : {
141 : : avg = computedInfos[AVERAGE_DET3D];
142 : : }
143 : : else
144 : : {
145 : : avg =0.;
146 : : int total_num_corners =0;
147 : : Matrix3D A[MSQ_MAX_NUM_VERT_PER_ENT];
148 : : for (size_t i=0; i<elementArray.size(); ++i) {
149 : : int nve = elementArray[i].corner_count();
150 : : elementArray[i].compute_corner_matrices(*this, A, nve, err);
151 : : MSQ_ERRZERO(err);
152 : : total_num_corners += nve;
153 : : for (int c=0; c<nve; ++c) {
154 : : avg += TargetCalculator::compute_Lambda(A[c], err);
155 : : MSQ_ERRZERO(err);
156 : : }
157 : : }
158 : :
159 : : avg = avg / total_num_corners;
160 : : computedInfos[AVERAGE_DET3D] = avg;
161 : : note_have_info(AVERAGE_DET3D);
162 : : }
163 : : return avg;
164 : : }
165 : : */
166 : :
167 : : /*! \fn PatchData::reorder()
168 : : Physically reorder the vertices and elements in the PatchData to improve
169 : : the locality of reference. This method implements a Reverse Breadth First
170 : : Search order starting with the vertex furthest from the origin. Other
171 : : orderings can also be implemented.
172 : : */
173 : 56 : void PatchData::reorder()
174 : : {
175 : : size_t i, j;
176 [ + - ]: 56 : const size_t num_vertex = num_nodes();
177 [ + - ]: 56 : const size_t num_elem = num_elements();
178 : :
179 : : // Step 1: Clear any cached data that will be invalidated by this
180 : 56 : vertexNormalIndices.clear();
181 : 56 : normalData.clear();
182 : : // vertexDomainDOF.clear();
183 : :
184 : : // Step 2: Make sure we have vertex-to-element adjacencies
185 [ + - ][ + - ]: 56 : if( !vertAdjacencyArray.size() ) generate_vertex_to_element_data();
186 : :
187 : : // Step 3: Do breadth-first search
188 [ + - ]: 56 : std::vector< bool > visited( num_vertex, false );
189 [ + - ]: 112 : std::vector< size_t > vertex_order( num_vertex );
190 : 56 : std::vector< size_t >::iterator q1_beg, q1_end, q2_end;
191 : 56 : q1_beg = q1_end = q2_end = vertex_order.begin();
192 : : // Outer loop will be done once for each disconnected chunk of mesh.
193 [ + - ][ + + ]: 112 : while( q1_beg != vertex_order.end() )
194 : : {
195 : : // Find vertex furthest from the origin
196 : 56 : double max = -1.0;
197 : 56 : size_t vtx_idx = num_vertex;
198 [ + + ]: 11669 : for( i = 0; i < num_vertex; ++i )
199 [ + - ][ + - ]: 11613 : if( !visited[i] )
200 : : {
201 [ + - ][ + - ]: 11613 : double dist = vertexArray[i].length_squared();
202 [ + + ]: 11613 : if( dist > max )
203 : : {
204 : 399 : max = dist;
205 : 11613 : vtx_idx = i;
206 : : }
207 : : }
208 [ - + ]: 56 : assert( vtx_idx < num_vertex );
209 : :
210 [ + - ][ + - ]: 56 : *q2_end++ = vtx_idx;
211 : : ;
212 [ + - ]: 56 : visited[vtx_idx] = true;
213 [ + - ][ + + ]: 460 : do
214 : : {
215 : 460 : q1_end = q2_end;
216 [ + - ][ + - ]: 12073 : for( ; q1_beg != q1_end; ++q1_beg )
[ + + ]
217 : : {
218 [ + - ][ + - ]: 11613 : size_t vtx_adj_offset = vertAdjacencyOffsets[*q1_beg];
219 [ + - ][ + - ]: 11613 : size_t vtx_adj_end = vertAdjacencyOffsets[*q1_beg + 1];
220 [ + + ]: 132355 : for( i = vtx_adj_offset; i < vtx_adj_end; ++i )
221 : : {
222 [ + - ]: 120742 : size_t elem = vertAdjacencyArray[i];
223 [ - + ]: 120742 : assert( elem < elementArray.size() );
224 [ + - ][ + - ]: 120742 : size_t num_elem_verts = elementArray[elem].node_count();
225 [ + - ][ + - ]: 120742 : size_t* elem_verts = elementArray[elem].get_vertex_index_array();
226 [ + + ]: 678072 : for( j = 0; j < num_elem_verts; ++j )
227 : : {
228 : 557330 : size_t elem_vert = elem_verts[j];
229 [ + - ][ + + ]: 557330 : if( !visited[elem_vert] )
230 : : {
231 [ + - ][ + - ]: 11557 : *q2_end++ = elem_vert;
232 : : ;
233 [ + - ]: 11557 : visited[elem_vert] = true;
234 : : }
235 : : }
236 : : }
237 : : }
238 : : } while( q2_end != q1_end );
239 : : }
240 : :
241 : : // Step 4: vertex_order contains the list of current vertex indices
242 : : // in the opposite of the order that they will occur in the
243 : : // reorderd patch. The following code will construct veretx_map
244 : : // from vertex_order with the following properties
245 : : // - vertex_map will be indexed by the current vertex index and
246 : : // contain the new index of that vertex (inverse of vertex_order)
247 : : // - the vertices will be grouped by their free/slave/fixed flag.
248 [ + - ]: 112 : std::vector< size_t > vertex_map( num_vertex );
249 : 56 : const size_t fixed_vtx_offset = numFreeVertices + numSlaveVertices;
250 : 56 : size_t free_idx = 0, slave_idx = numFreeVertices, fixed_idx = fixed_vtx_offset;
251 [ + + ]: 11669 : for( i = 1; i <= num_vertex; ++i )
252 : : {
253 [ + - ]: 11613 : size_t vtx_idx = vertex_order[num_vertex - i];
254 [ + + ]: 11613 : if( vtx_idx < numFreeVertices )
255 [ + - ]: 6310 : vertex_map[vtx_idx] = free_idx++;
256 [ - + ]: 5303 : else if( vtx_idx < fixed_vtx_offset )
257 [ # # ]: 0 : vertex_map[vtx_idx] = slave_idx++;
258 : : else
259 [ + - ]: 5303 : vertex_map[vtx_idx] = fixed_idx++;
260 : : }
261 : : // make sure everything adds up
262 [ - + ]: 56 : assert( free_idx == numFreeVertices );
263 [ - + ]: 56 : assert( slave_idx == fixed_vtx_offset );
264 [ - + ]: 56 : assert( fixed_idx == num_vertex );
265 : :
266 : : // Step 5: compute element permutation
267 : : // initialize all to "num_elem" to indicate unvisited
268 [ + - ]: 112 : std::vector< size_t > element_map( num_elem, num_elem );
269 : 56 : size_t elem_idx = 0;
270 [ + + ]: 11669 : for( i = 1; i <= num_vertex; ++i )
271 : : {
272 [ + - ]: 11613 : size_t vtx_idx = vertex_order[num_vertex - i];
273 [ + - ]: 11613 : size_t vtx_adj_offset = vertAdjacencyOffsets[vtx_idx];
274 [ + - ]: 11613 : size_t vtx_adj_end = vertAdjacencyOffsets[vtx_idx + 1];
275 [ + + ]: 132355 : for( j = vtx_adj_offset; j < vtx_adj_end; ++j )
276 : : {
277 [ + - ]: 120742 : size_t elem = vertAdjacencyArray[j];
278 [ + - ][ + + ]: 120742 : if( element_map[elem] == num_elem ) element_map[elem] = elem_idx++;
[ + - ]
279 : : }
280 : : }
281 : : // make sure everything adds up
282 [ - + ]: 56 : assert( elem_idx == num_elem );
283 : :
284 : : // Step 6: Permute the vertices
285 [ + - ]: 112 : std::vector< MsqVertex > new_vertex_array( num_vertex );
286 [ + - ]: 112 : std::vector< Mesh::VertexHandle > new_vtx_handle_array( num_vertex );
287 [ + + ]: 11669 : for( i = 0; i < num_vertex; ++i )
288 : : {
289 [ + - ]: 11613 : size_t new_idx = vertex_map[i];
290 [ + - ][ + - ]: 11613 : new_vertex_array[new_idx] = vertexArray[i];
[ + - ]
291 [ + - ][ + - ]: 11613 : new_vtx_handle_array[new_idx] = vertexHandlesArray[i];
292 : : }
293 : 56 : vertexArray.swap( new_vertex_array );
294 : 56 : vertexHandlesArray.swap( new_vtx_handle_array );
295 : :
296 : : // Step 7: Permute the elements and vertex indices for the elements
297 [ + - ]: 112 : std::vector< MsqMeshEntity > new_elem_array( num_elem );
298 [ + - ]: 112 : std::vector< Mesh::ElementHandle > new_elem_handle_array( num_elem );
299 [ + + ]: 27963 : for( i = 0; i < num_elem; ++i )
300 : : {
301 [ - + ]: 27907 : assert( i < elementArray.size() );
302 [ + - ][ + - ]: 27907 : size_t vert_count = elementArray[i].node_count();
303 [ + - ][ + - ]: 27907 : size_t* conn_array = elementArray[i].get_vertex_index_array();
304 [ + + ]: 148649 : for( j = 0; j < vert_count; ++j )
305 : : {
306 [ + - ]: 120742 : conn_array[j] = vertex_map[conn_array[j]];
307 : : }
308 : :
309 [ + - ]: 27907 : size_t new_idx = element_map[i];
310 [ - + ]: 27907 : assert( new_idx < num_elem );
311 [ + - ][ + - ]: 27907 : new_elem_array[new_idx] = elementArray[i];
312 [ + - ][ + - ]: 27907 : new_elem_handle_array[new_idx] = elementHandlesArray[i];
313 : : }
314 : 56 : elementArray.swap( new_elem_array );
315 : 56 : elementHandlesArray.swap( new_elem_handle_array );
316 : :
317 : : // Step 8: Clear no-longer-valid vertex-to-element adjacency info.
318 [ + - ]: 56 : if( vertAdjacencyOffsets.size() )
319 : : {
320 : 56 : vertAdjacencyOffsets.clear();
321 : 56 : vertAdjacencyArray.clear();
322 [ + - ]: 56 : generate_vertex_to_element_data();
323 : : }
324 : :
325 [ + - ]: 112 : notify_new_patch();
326 : 56 : }
327 : :
328 : : /*!
329 : : PatchData::move_free_vertices_constrained() moves the free vertices
330 : : (see MsqVertex::is_free() ) as specified by the search direction (dk)
331 : : and scale factor (step_size). After being moved, the vertices are
332 : : projected onto the appropriate geometry. Fixed vertices are not moved
333 : : regardless of their corresponding dk direction.
334 : : It is often useful to use the create_coords_momento() function before
335 : : calling this function.
336 : : Compile with -DMSQ_DBG3 to check that fixed vertices
337 : : are not moved with that call.
338 : :
339 : : \param dk: must be a [nb_vtx] array of Vector3D that contains
340 : : the direction in which to move each vertex. Fixed vertices moving
341 : : direction should be zero, although fixed vertices will not be
342 : : moved regardless of their corresponding dk value.
343 : : \param nb_vtx is the number of vertices to move. must corresponds
344 : : to the number of vertices in the PatchData.
345 : : \param step_size will multiply the moving direction given in dk
346 : : for each vertex.
347 : : */
348 : 216870 : void PatchData::move_free_vertices_constrained( Vector3D dk[], size_t nb_vtx, double step_size, MsqError& err )
349 : : {
350 [ - + ]: 216870 : if( nb_vtx != num_free_vertices() )
351 : : {
352 : : MSQ_SETERR( err )
353 [ # # ]: 0 : ( "The directional vector must be of length numVertices.", MsqError::INVALID_ARG );
354 : 0 : return;
355 : : }
356 : :
357 : : size_t i;
358 [ + + ]: 899044 : for( i = 0; i < num_free_vertices(); ++i )
359 : : {
360 [ + - ][ + - ]: 682174 : vertexArray[i] += ( step_size * dk[i] );
361 [ - + ][ # # ]: 682174 : snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
[ - + ]
362 : : }
363 : :
364 [ + + ]: 216870 : if( numSlaveVertices )
365 : : {
366 [ - + ][ # # ]: 216870 : update_slave_node_coordinates( err );MSQ_CHKERR( err );
367 : : }
368 : : }
369 : :
370 : : /*! set_free_vertices_constrained is similar to
371 : : PatchData::move_free_vertices_constrained() except the original vertex positions
372 : : are those stored in the PatchDataVerticesMemento instead of the actual vertex
373 : : position stored in the PatchData Vertex array. The final location is stored
374 : : in the PatchData; the PatchDataVerticesMemento is unchanged.
375 : :
376 : : \param dk: must be a [nb_vtx] array of Vector3D that contains
377 : : the direction in which to move each vertex. Fixed vertices moving
378 : : direction should be zero, although fixed vertices will not be
379 : : moved regardless of their corresponding dk value.
380 : : \param nb_vtx is the number of vertices to move. must corresponds
381 : : to the number of vertices in the PatchData.
382 : : \param step_size will multiply the moving direction given in dk
383 : : for each vertex.
384 : : */
385 : 172104 : void PatchData::set_free_vertices_constrained( PatchDataVerticesMemento* memento, Vector3D dk[], size_t nb_vtx,
386 : : double step_size, MsqError& err )
387 : : {
388 [ + - ][ - + ]: 172104 : if( memento->originator != this || nb_vtx != num_free_vertices() )
[ - + ]
389 : : {
390 [ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_ARG );
391 : 0 : return;
392 : : }
393 : :
394 : : size_t i;
395 [ + + ]: 1429257 : for( i = 0; i < num_free_vertices(); ++i )
396 : : {
397 [ + - ][ + - ]: 1257153 : vertexArray[i] = memento->vertices[i] + ( step_size * dk[i] );
[ + - ][ + - ]
398 [ - + ][ # # ]: 1257153 : snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
[ - + ]
399 : : }
400 : :
401 [ - + ]: 172104 : if( numSlaveVertices )
402 : : {
403 [ # # ][ # # ]: 172104 : update_slave_node_coordinates( err );MSQ_CHKERR( err );
404 : : }
405 : :
406 : : // Checks that moving direction is zero for fixed vertices.
407 : : if( MSQ_DBG( 3 ) )
408 : : {
409 : : for( i = 0; i < num_nodes(); ++i )
410 : : {
411 : : if( dk[i].length_squared() != 0.0 )
412 : : {
413 : : MSQ_DBGOUT( 3 ) << "dk[" << i << "]: " << dk[i] << endl;
414 : : MSQ_DBGOUT( 3 ) << "moving a fixed vertex." << endl;
415 : : }
416 : : }
417 : : }
418 : : }
419 : :
420 : 0 : static void project_to_plane( Vector3D& vect, const Vector3D& norm )
421 : : {
422 : 0 : double len_sqr = norm % norm;
423 [ # # ][ # # ]: 0 : if( len_sqr > DBL_EPSILON ) vect -= norm * ( ( norm % vect ) / len_sqr );
424 : 0 : }
425 : :
426 : 0 : void PatchData::project_gradient( std::vector< Vector3D >& gradient, MsqError& err )
427 : : {
428 [ # # ][ # # ]: 0 : if( !domain_set() ) return;
429 : :
430 [ # # ][ # # ]: 0 : if( gradient.size() != num_free_vertices() )
431 : : {
432 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_ARG );
433 : 0 : return;
434 : : }
435 : :
436 [ # # ]: 0 : if( normalData.empty() )
437 : : {
438 [ # # ][ # # ]: 0 : update_cached_normals( err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
439 : : }
440 : :
441 [ # # ]: 0 : Vector3D norm;
442 [ # # ][ # # ]: 0 : for( size_t i = 0; i < num_free_vertices(); ++i )
443 : : {
444 [ # # ]: 0 : if( vertexNormalIndices.empty() )
445 : : { // whole mesh on single 2D domain
446 [ # # ][ # # ]: 0 : norm = normalData[i];
447 [ # # ][ # # ]: 0 : project_to_plane( gradient[i], norm );
448 : : }
449 [ # # ][ # # ]: 0 : else if( vertexDomainDOF[i] == 2 )
450 : : { // vertex on surface
451 [ # # ][ # # ]: 0 : norm = normalData[vertexNormalIndices[i]];
[ # # ]
452 [ # # ][ # # ]: 0 : project_to_plane( gradient[i], norm );
453 : : }
454 [ # # ][ # # ]: 0 : else if( vertexDomainDOF[i] == 1 )
455 : : {
456 : : size_t j, num_elem;
457 [ # # ][ # # ]: 0 : const size_t* elems = get_vertex_element_adjacencies( i, num_elem, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
458 [ # # ]: 0 : for( j = 0; j < num_elem; ++j )
459 : : {
460 [ # # ][ # # ]: 0 : if( 2 == TopologyInfo::dimension( element_by_index( elems[j] ).get_element_type() ) )
[ # # ][ # # ]
461 : : {
462 [ # # ][ # # ]: 0 : norm = vertexArray[i];
463 [ # # ][ # # ]: 0 : get_domain()->element_normal_at( elementHandlesArray[elems[j]], norm );
[ # # ]
464 [ # # ][ # # ]: 0 : project_to_plane( gradient[i], norm );
465 : : }
466 : : }
467 : : }
468 : : }
469 : : }
470 : :
471 : : /*! Finds the maximum movement (in the distance norm) of the vertices in a
472 : : patch. The previous vertex positions are givena as a
473 : : PatchDataVerticesMemento (memento). The distance squared which each
474 : : vertex has moved is then calculated, and the largest of those distances
475 : : is returned. This function only considers the movement of vertices
476 : : that are currently 'free'.
477 : : \param memento a memento of this patch's vertex position at some
478 : : (prior) time in the optimization.
479 : : */
480 : 827093 : double PatchData::get_max_vertex_movement_squared( PatchDataVerticesMemento* memento, MsqError& )
481 : : {
482 : 827093 : double max_dist = 0.0;
483 [ + + ]: 1887508 : for( size_t i = 0; i < memento->vertices.size(); ++i )
484 : : {
485 [ + - ]: 1060415 : double temp_dist = ( vertexArray[i] - memento->vertices[i] ).length_squared();
486 [ + + ]: 1060415 : if( temp_dist > max_dist ) max_dist = temp_dist;
487 : : }
488 : 827093 : return max_dist;
489 : : }
490 : :
491 : : /*!
492 : : */
493 : 0 : void PatchData::set_all_vertices_soft_fixed( MsqError& /*err*/ )
494 : : {
495 [ # # ]: 0 : for( size_t i = 0; i < num_free_vertices(); ++i )
496 : 0 : vertexArray[i].set_soft_fixed_flag();
497 : 0 : }
498 : :
499 : : /*!
500 : : */
501 : 49295 : void PatchData::set_free_vertices_soft_fixed( MsqError& /*err*/ )
502 : : {
503 [ + + ]: 98590 : for( size_t i = 0; i < num_free_vertices(); ++i )
504 : : {
505 [ + - ]: 49295 : if( vertexArray[i].is_free_vertex() ) vertexArray[i].set_soft_fixed_flag();
506 : : }
507 : 49295 : }
508 : :
509 : : /*!
510 : : */
511 : 776818 : void PatchData::set_all_vertices_soft_free( MsqError& /*err*/ )
512 : : {
513 [ + + ]: 7768598 : for( size_t i = 0; i < num_nodes(); ++i )
514 : 6991780 : vertexArray[i].remove_soft_fixed_flag();
515 : 776818 : }
516 : :
517 : : /*! Get coordinates of element vertices, in canonical order.
518 : :
519 : : \param elem_index The element index in the Patch
520 : : \param coords This vector will have the coordinates appended to it.
521 : : If necessary, make sure to clear the vector before calling the function.
522 : : */
523 : 0 : void PatchData::get_element_vertex_coordinates( size_t elem_index, std::vector< Vector3D >& coords, MsqError& /*err*/ )
524 : : {
525 : : // Check index
526 [ # # ]: 0 : if( elem_index >= num_elements() ) return;
527 : :
528 : : // Ask the element for its vertex indices
529 [ # # ]: 0 : assert( elem_index < elementArray.size() );
530 : 0 : const size_t* vertex_indices = elementArray[elem_index].get_vertex_index_array();
531 : :
532 : : // Get the coords for each indicated vertex
533 : 0 : size_t num_verts = elementArray[elem_index].vertex_count();
534 : 0 : coords.reserve( coords.size() + num_verts );
535 [ # # ]: 0 : for( size_t i = 0; i < num_verts; i++ )
536 [ # # ]: 0 : coords.push_back( Vector3D( vertexArray[vertex_indices[i]] ) );
537 : : }
538 : :
539 : : /*! This is an inefficient way of retrieving vertex_indices.
540 : : Use PatchData::get_element_array followed by
541 : : MsqMeshEntity::get_vertex_index_array() if you don't need
542 : : to fill an STL vector.
543 : : */
544 : 0 : void PatchData::get_element_vertex_indices( size_t elem_index, std::vector< size_t >& vertex_indices,
545 : : MsqError& /*err*/ )
546 : : {
547 : : // Ask the element for its vertex indices
548 [ # # ]: 0 : assert( elem_index < elementArray.size() );
549 : 0 : elementArray[elem_index].get_vertex_indices( vertex_indices );
550 : 0 : }
551 : :
552 : 0 : void PatchData::get_vertex_element_indices( size_t vertex_index, std::vector< size_t >& elem_indices, MsqError& err )
553 : : {
554 : : size_t count;
555 : : const size_t* ptr;
556 [ # # ]: 0 : ptr = get_vertex_element_adjacencies( vertex_index, count, err );
557 [ # # ]: 0 : elem_indices.resize( count );
558 [ # # ]: 0 : std::copy( ptr, ptr + count, elem_indices.begin() );
559 : 0 : }
560 : :
561 : 0 : void PatchData::get_vertex_element_indices( size_t vertex_index, unsigned element_dimension,
562 : : std::vector< size_t >& elem_indices, MsqError& err )
563 : : {
564 : 0 : elem_indices.clear();
565 : : size_t count;
566 : : const size_t* ptr;
567 [ # # ]: 0 : ptr = get_vertex_element_adjacencies( vertex_index, count, err );
568 [ # # ]: 0 : for( const size_t* const end = ptr + count; ptr != end; ++ptr )
569 : : {
570 [ # # ]: 0 : assert( *ptr < elementArray.size() );
571 [ # # ][ # # ]: 0 : const EntityTopology type = elementArray[*ptr].get_element_type();
572 [ # # ]: 0 : const unsigned dim = TopologyInfo::dimension( type );
573 [ # # ][ # # ]: 0 : if( dim == element_dimension ) elem_indices.push_back( *ptr );
574 : : }
575 : 0 : }
576 : :
577 : 61225 : const size_t* PatchData::get_vertex_element_adjacencies( size_t vertex_index, size_t& array_len_out, MsqError& )
578 : : {
579 : : // Make sure we've got the data
580 [ + + ]: 61225 : if( vertAdjacencyArray.empty() ) { generate_vertex_to_element_data(); }
581 : :
582 : 61225 : const size_t begin = vertAdjacencyOffsets[vertex_index];
583 : 61225 : const size_t end = vertAdjacencyOffsets[vertex_index + 1];
584 : 61225 : array_len_out = end - begin;
585 : 61225 : return &vertAdjacencyArray[begin];
586 : : }
587 : :
588 : : /*!
589 : : \brief This function fills a vector<size_t> with the indices
590 : : to vertices connected to the given vertex by an edge. If vert_indices
591 : : is not initially empty, the function will not delete the current
592 : : contents. Instead, it will append the new indices at the end of
593 : : the vector.
594 : :
595 : : */
596 : 1132545 : void PatchData::get_adjacent_vertex_indices( size_t vertex_index, std::vector< size_t >& vert_indices,
597 : : MsqError& err ) const
598 : : {
599 : 1132545 : bitMap.clear();
600 [ + - ][ + - ]: 1132545 : bitMap.resize( num_nodes(), false );
601 : :
602 : : const size_t* conn;
603 : : size_t conn_idx, curr_vtx_idx;
604 : : const unsigned* adj;
605 : : unsigned num_adj, i;
606 : 1132545 : std::vector< MsqMeshEntity >::const_iterator e;
607 [ + - ][ + - ]: 29380576 : for( e = elementArray.begin(); e != elementArray.end(); ++e )
[ + + ]
608 : : {
609 [ + - ][ + - ]: 28248031 : conn = e->get_vertex_index_array();
610 [ + - ][ + - ]: 28248031 : conn_idx = std::find( conn, conn + e->node_count(), vertex_index ) - conn;
[ + - ]
611 [ + - ][ + - ]: 28248031 : if( conn_idx == e->node_count() ) continue;
[ + + ]
612 : :
613 : : // If a higher-order node, return corners of side/face
614 : : // that node is in the center of.
615 [ + - ][ + - ]: 11157403 : EntityTopology type = e->get_element_type();
616 [ + - ][ + + ]: 11157403 : if( conn_idx >= TopologyInfo::corners( type ) && type != POLYGON )
[ - + ][ - + ]
617 : : {
618 : : unsigned dim, id;
619 [ # # ][ # # ]: 1132545 : TopologyInfo::side_from_higher_order( type, e->node_count(), conn_idx, dim, id, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
620 [ # # ]: 0 : adj = TopologyInfo::side_vertices( type, dim, id, num_adj );
621 : : }
622 : : else
623 : : {
624 [ + - ][ + - ]: 11157403 : EntityTopology topo = e->get_element_type();
625 [ + + ]: 11157403 : if( topo == POLYGON )
626 : : {
627 [ + - ][ + - ]: 1500 : unsigned number_of_nodes = e->node_count();
628 : 1500 : num_adj = 2; // always 2 for a polygon
629 : : unsigned vert_adj[2];
630 : 1500 : vert_adj[0] = ( conn_idx + 1 ) % number_of_nodes;
631 : 1500 : vert_adj[1] = ( conn_idx + number_of_nodes - 1 ) % number_of_nodes;
632 [ + + ]: 4500 : for( i = 0; i < num_adj; ++i )
633 : : {
634 : 3000 : curr_vtx_idx = conn[vert_adj[i]]; // get index into patch vertex list
635 [ + - ][ + + ]: 3000 : if( !bitMap[curr_vtx_idx] )
636 : : {
637 [ + - ]: 1500 : vert_indices.push_back( curr_vtx_idx );
638 [ + - ]: 1500 : bitMap[curr_vtx_idx] = true;
639 : : }
640 : : }
641 : : }
642 : : else
643 : : {
644 [ + - ]: 11155903 : adj = TopologyInfo::adjacent_vertices( topo, conn_idx, num_adj );
645 [ + + ]: 41405017 : for( i = 0; i < num_adj; ++i )
646 : : {
647 : 30249114 : curr_vtx_idx = conn[adj[i]]; // get index into patch vertex list
648 [ + - ][ + + ]: 30249114 : if( !bitMap[curr_vtx_idx] )
649 : : {
650 [ + - ]: 7852599 : vert_indices.push_back( curr_vtx_idx );
651 [ + - ]: 7852599 : bitMap[curr_vtx_idx] = true;
652 : : }
653 : : }
654 : : }
655 : : }
656 : : }
657 : : }
658 : :
659 : : /*! Fills a vector of indices into the entities array. The entities
660 : : in the vector are connected the given entity (ent_ind) via an
661 : : n-diminsional entity (where 'n' is a given integer).
662 : : Thus, if n = 0, the entities must be connected via a vertex.
663 : : If n = 1, the entities must be connected via an edge.
664 : : If n = 2, the entities must be connected via a two-dimensional element.
665 : : NOTE: if n is 2 and the elements in the entity array are
666 : : two-dimensional, no entities should meet this criterion.
667 : : The adj_ents vector is cleared at the beginning of the call.
668 : :
669 : : */
670 : 0 : void PatchData::get_adjacent_entities_via_n_dim( int n, size_t ent_ind, std::vector< size_t >& adj_ents, MsqError& err )
671 : : {
672 : : // reset the vector
673 : 0 : adj_ents.clear();
674 : : // vertices of this entity (given by ent_ind)
675 [ # # ]: 0 : vector< size_t > verts;
676 : : // vector to store elements attached to the vertices in verts
677 [ # # ][ # # ]: 0 : vector< size_t > elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
[ # # ][ # #
# # # # ]
678 : : // length of above vectos
679 : : int length_elem_on_vert[MSQ_MAX_NUM_VERT_PER_ENT];
680 : : // get verts on this element
681 [ # # ]: 0 : get_element_vertex_indices( ent_ind, verts, err );
682 : 0 : int num_vert = verts.size();
683 : 0 : int i = 0;
684 : 0 : int j = 0;
685 [ # # ]: 0 : for( i = 0; i < num_vert; ++i )
686 : : {
687 : : // get elements on the vertices in verts and the number of vertices
688 [ # # ][ # # ]: 0 : get_vertex_element_indices( verts[i], elem_on_vert[i], err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
689 : 0 : length_elem_on_vert[i] = elem_on_vert[i].size();
690 : : }
691 : : // this_ent is the index for an entity which is a candidate to be placed
692 : : // into adj_ents
693 : : size_t this_ent;
694 : : // num of times this_ent has been found in the vectors of entity indices
695 : 0 : int counter = 0;
696 : 0 : i = 0;
697 : : // loop of each vert on ent_ind
698 [ # # ][ # # ]: 0 : while( i < num_vert )
699 : : {
700 : : // loop over each ent connected to vert
701 : 0 : j = 0;
702 [ # # ]: 0 : while( j < length_elem_on_vert[i] )
703 : : {
704 : : // get candidate element
705 [ # # ]: 0 : this_ent = elem_on_vert[i][j];
706 : : // if we haven't already consider this ent
707 [ # # ]: 0 : if( this_ent != ent_ind )
708 : : {
709 : : // if this_ent occurred earlier we would have already considered it
710 : : // so start at i and j+1
711 : 0 : int k1 = i;
712 : 0 : int k2 = j + 1;
713 : : // this_ent has occured once so far
714 : 0 : counter = 1;
715 : : // while k1 < num_vert
716 [ # # ]: 0 : while( k1 < num_vert )
717 : : {
718 : : // loop over entries in the elem on vert vector
719 [ # # ]: 0 : while( k2 < length_elem_on_vert[k1] )
720 : : {
721 : : // if it matches this_ent
722 [ # # ][ # # ]: 0 : if( elem_on_vert[k1][k2] == this_ent )
723 : : {
724 : : // mark it as 'seen', by making it the same as ent_ind
725 : : // i.e., the entity passed to us.
726 [ # # ]: 0 : elem_on_vert[k1][k2] = ent_ind;
727 : 0 : ++counter;
728 : : // do not look at remaining elems in this vector
729 : 0 : k2 += length_elem_on_vert[k1];
730 : : }
731 : : else
732 : 0 : ++k2;
733 : : }
734 : 0 : ++k1;
735 : 0 : k2 = 0;
736 : : }
737 : : // if this_ent occured enough times and isn't ent_ind
738 [ # # ][ # # ]: 0 : if( counter > n && this_ent != ent_ind ) { adj_ents.push_back( this_ent ); }
[ # # ]
739 : : }
740 : 0 : ++j;
741 : : }
742 : 0 : ++i;
743 : 0 : }
744 : : }
745 : :
746 : : /*! \fn PatchData::update_mesh(MsqError &err)
747 : :
748 : : \brief This function copies to the TSTT mesh the changes made to the
749 : : free vertices / elements of the PatchData object.
750 : :
751 : : */
752 : 829946 : void PatchData::update_mesh( MsqError& err, const TagHandle* tag )
753 : : {
754 [ - + ]: 829946 : if( !myMesh )
755 : : {
756 [ # # ]: 0 : MSQ_SETERR( err )( "Cannot update mesh on temporary patch.\n", MsqError::INTERNAL_ERROR );
757 : 0 : return;
758 : : }
759 : :
760 : 829946 : const size_t not_fixed = numFreeVertices + numSlaveVertices;
761 [ + + ]: 829946 : if( tag )
762 : : {
763 [ + + ]: 8 : for( size_t i = 0; i < not_fixed; ++i )
764 : : {
765 [ - + ][ # # ]: 4 : myMesh->tag_set_vertex_data( *tag, 1, &vertexHandlesArray[i], vertexArray[i].to_array(), err );MSQ_ERRRTN( err );
[ - + ]
766 : : }
767 : : }
768 : : else
769 : : {
770 [ + + ]: 1691143 : for( size_t i = 0; i < not_fixed; ++i )
771 : : {
772 [ - + ][ # # ]: 861201 : myMesh->vertex_set_coordinates( vertexHandlesArray[i], vertexArray[i], err );MSQ_ERRRTN( err );
[ - + ]
773 : : }
774 : : }
775 : :
776 [ + + ]: 8363532 : for( size_t i = 0; i < vertexArray.size(); ++i )
777 : : {
778 [ - + ][ # # ]: 7533586 : myMesh->vertex_set_byte( vertexHandlesArray[i], vertexArray[i].get_flags(), err );MSQ_ERRRTN( err );
[ - + ]
779 : : }
780 : : }
781 : :
782 : 182 : void PatchData::update_slave_node_coordinates( MsqError& err )
783 : : {
784 : : // update slave vertices
785 [ + - ][ - + ]: 364 : if( 0 == num_slave_vertices() ) return;
786 : :
787 : : // Set a mark on every slave vertex. We'll clear the marks as we
788 : : // set the vertex coordinates. This way we can check that all
789 : : // vertices got updated.
790 [ + - ][ + - ]: 182 : const size_t vert_end = num_free_vertices() + num_slave_vertices();
791 [ + - ][ + + ]: 4837 : for( size_t i = num_free_vertices(); i < vert_end; ++i )
792 [ + - ][ + - ]: 4655 : vertexArray[i].flags() |= MsqVertex::MSQ_MARK;
793 : :
794 : : // For each element, calculate slave vertex coordinates from
795 : : // mapping function.
796 : 182 : const int ARR_SIZE = 27;
797 : : double coeff[ARR_SIZE];
798 : : size_t index[ARR_SIZE];
799 [ + - ][ + + ]: 6863 : for( size_t i = 0; i < num_elements(); ++i )
800 : : {
801 [ + - ]: 6681 : MsqMeshEntity& elem = element_by_index( i );
802 [ + - ]: 6681 : const int num_corner = elem.vertex_count();
803 [ + - ]: 6681 : const int num_node = elem.node_count();
804 [ - + ]: 6681 : assert( num_node < ARR_SIZE );
805 : :
806 [ + - ]: 6681 : const EntityTopology type = elem.get_element_type();
807 [ + - ]: 6681 : const MappingFunction* const mf = get_mapping_function( type );
808 [ + - ][ - + ]: 6681 : if( 0 == mf || num_node == num_corner ) continue;
809 : :
810 [ + - ]: 6681 : const size_t* conn = elem.get_vertex_index_array();
811 : :
812 : : // for each higher-order non-slave node, set bit indicating
813 : : // that mapping function is a function of the non-slave node
814 : : // coordinates
815 [ + - ]: 6681 : NodeSet ho_bits = non_slave_node_set( i );
816 : :
817 : : // for each higher-order slave node
818 [ + + ]: 27255 : for( int k = num_corner; k < num_node; ++k )
819 : : {
820 [ + - ][ + + ]: 25229 : if( !is_vertex_slave( conn[k] ) ) continue;
821 : :
822 : : // check if we already did this one for an adjacent element
823 [ + - ]: 9310 : MsqVertex& vert = vertexArray[conn[k]];
824 [ + - ][ + + ]: 9310 : if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;
825 : :
826 : : // what is this node a mid node of (i.e. face 1, edge 2, etc.)
827 [ + - ][ + - ]: 4655 : Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
828 : :
829 : : // evaluate mapping function at logical loction of HO node.
830 : : size_t num_coeff;
831 [ + - ][ + - ]: 4655 : mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
832 [ + - ][ + - ]: 4655 : mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
833 : :
834 : : // calulate new coordinates for slave node
835 [ - + ]: 4655 : assert( num_coeff > 0 );
836 [ + - ][ + - ]: 4655 : vert = coeff[0] * vertex_by_index( conn[index[0]] );
[ + - ]
837 [ + + ]: 9310 : for( size_t j = 1; j < num_coeff; ++j )
838 [ + - ][ + - ]: 4655 : vert += coeff[j] * vertex_by_index( conn[index[j]] );
[ + - ]
839 : :
840 : : // clear mark
841 [ + - ]: 4655 : vert.flags() &= ~MsqVertex::MSQ_MARK;
842 : : }
843 : : }
844 : :
845 : : // make sure we set the coordinates for every slave node
846 [ + - ][ + + ]: 4837 : for( size_t i = num_free_vertices(); i < vert_end; ++i )
847 : : {
848 [ + - ][ + - ]: 4655 : if( vertex_by_index( i ).is_flag_set( MsqVertex::MSQ_MARK ) )
[ - + ]
849 : : {
850 : : MSQ_SETERR( err )
851 : : ( MsqError::INVALID_MESH, "No element with mapping function adjacent to slave vertex %lu (%lu)\n",
852 [ # # ][ # # ]: 0 : (unsigned long)i, (unsigned long)get_vertex_handles_array()[i] );
[ # # ]
853 : : // make sure we finish with all marks cleared
854 [ # # ][ # # ]: 0 : vertexArray[i].flags() &= ~MsqVertex::MSQ_MARK;
855 : : }
856 : : }
857 : :
858 : : // snap vertices to domain
859 [ + - ][ + - ]: 182 : if( domain_set() )
860 : : {
861 [ + - ][ + + ]: 4837 : for( size_t i = num_free_vertices(); i < vert_end; ++i )
862 : : {
863 [ + - ][ + - ]: 4655 : snap_vertex_to_domain( i, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
864 : : }
865 : : }
866 : : }
867 : :
868 : 0 : void PatchData::update_slave_node_coordinates( const size_t* elements, size_t num_elems, MsqError& err )
869 : : {
870 : : // update slave vertices
871 [ # # ][ # # ]: 0 : if( 0 == num_slave_vertices() ) return;
872 : :
873 : : // set a mark on each vertex so we don't update shared
874 : : // vertices more than once.
875 [ # # ]: 0 : for( size_t i = 0; i < num_elems; ++i )
876 : : {
877 [ # # ]: 0 : MsqMeshEntity& elem = element_by_index( elements[i] );
878 [ # # ]: 0 : const int num_corner = elem.vertex_count();
879 [ # # ]: 0 : const int num_node = elem.node_count();
880 [ # # ]: 0 : const size_t* conn = elem.get_vertex_index_array();
881 [ # # ]: 0 : for( int j = num_corner; j < num_node; ++j )
882 [ # # ][ # # ]: 0 : vertexArray[conn[j]].flags() |= MsqVertex::MSQ_MARK;
883 : : }
884 : :
885 : : // For each element, calculate slave vertex coordinates from
886 : : // mapping function.
887 : 0 : const int ARR_SIZE = 27;
888 : : double coeff[ARR_SIZE];
889 : : size_t index[ARR_SIZE];
890 [ # # ]: 0 : for( size_t i = 0; i < num_elems; ++i )
891 : : {
892 [ # # ]: 0 : MsqMeshEntity& elem = element_by_index( elements[i] );
893 [ # # ]: 0 : const int num_corner = elem.vertex_count();
894 [ # # ]: 0 : const int num_node = elem.node_count();
895 [ # # ]: 0 : assert( num_node < ARR_SIZE );
896 : :
897 [ # # ]: 0 : const EntityTopology type = elem.get_element_type();
898 [ # # ]: 0 : const MappingFunction* const mf = get_mapping_function( type );
899 [ # # ][ # # ]: 0 : if( 0 == mf || num_node == num_corner ) continue;
900 : :
901 [ # # ]: 0 : const size_t* conn = elem.get_vertex_index_array();
902 : :
903 : : // for each higher-order non-slave node, set bit indicating
904 : : // that mapping function is a function of the non-slave node
905 : : // coordinates
906 [ # # ]: 0 : NodeSet ho_bits = non_slave_node_set( i );
907 : :
908 : : // for each higher-order slave node
909 [ # # ]: 0 : for( int k = num_corner; k < num_node; ++k )
910 : : {
911 [ # # ][ # # ]: 0 : if( !is_vertex_slave( conn[k] ) ) continue;
912 : :
913 : : // check if we already did this one for an adjacent element
914 [ # # ]: 0 : MsqVertex& vert = vertexArray[conn[k]];
915 [ # # ][ # # ]: 0 : if( !vert.is_flag_set( MsqVertex::MSQ_MARK ) ) continue;
916 : :
917 : : // what is this node a mid node of (i.e. face 1, edge 2, etc.)
918 [ # # ][ # # ]: 0 : Sample loc = TopologyInfo::sample_from_node( type, elem.node_count(), k, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
919 : :
920 : : // evaluate mapping function at logical loction of HO node.
921 : : size_t num_coeff;
922 [ # # ][ # # ]: 0 : mf->coefficients( loc, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
923 [ # # ][ # # ]: 0 : mf->convert_connectivity_indices( num_node, index, num_coeff, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
924 : :
925 : : // calulate new coordinates for slave node
926 [ # # ]: 0 : assert( num_coeff > 0 );
927 [ # # ][ # # ]: 0 : vert = coeff[0] * vertex_by_index( conn[index[0]] );
[ # # ]
928 [ # # ]: 0 : for( size_t j = 1; j < num_coeff; ++j )
929 [ # # ][ # # ]: 0 : vert += coeff[j] * vertex_by_index( conn[index[j]] );
[ # # ]
930 : :
931 : : // snap vertices to domain
932 [ # # ][ # # ]: 0 : if( domain_set() )
933 : : {
934 [ # # ][ # # ]: 0 : snap_vertex_to_domain( conn[k], err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
935 : : }
936 : :
937 : : // clear mark
938 [ # # ]: 0 : vert.flags() &= ~MsqVertex::MSQ_MARK;
939 : : }
940 : : }
941 : : }
942 : :
943 : 36783 : void PatchData::generate_vertex_to_element_data()
944 : : {
945 : : MSQ_FUNCTION_TIMER( "PatchData::generate_vertex_to_element_data" );
946 : :
947 : : // Skip if data already exists
948 [ - + ]: 73566 : if( !vertAdjacencyArray.empty() ) return;
949 : :
950 : : // Skip if patch is empty
951 [ + - ][ - + ]: 36783 : if( 0 == num_nodes() ) return;
952 : :
953 : : // Allocate offset array
954 : 36783 : vertAdjacencyOffsets.clear();
955 [ + - ][ + - ]: 36783 : vertAdjacencyOffsets.resize( num_nodes() + 1, 0 );
956 : :
957 : : // Temporarily use offsets array to hold per-vertex element count
958 : 36783 : std::vector< MsqMeshEntity >::iterator elem_iter;
959 : 36783 : const std::vector< MsqMeshEntity >::iterator elem_end = elementArray.end();
960 [ + - ][ + - ]: 270463 : for( elem_iter = elementArray.begin(); elem_iter != elem_end; ++elem_iter )
[ + + ]
961 : : {
962 [ + - ][ + - ]: 233680 : size_t* conn_iter = elem_iter->get_vertex_index_array();
963 [ + - ][ + - ]: 233680 : const size_t* conn_end = conn_iter + elem_iter->node_count();
964 [ + + ]: 1190628 : for( ; conn_iter != conn_end; ++conn_iter )
965 [ + - ]: 956948 : ++vertAdjacencyOffsets[*conn_iter];
966 : : }
967 : :
968 : : // Convert counts to end indices.
969 : : // When done, vertAdjacencyOffsets will contain, for each vertex,
970 : : // one more than the *last* index for that vertex's data in the
971 : : // adjacency array. This is *not* the final state for this data.
972 : : // See comments for next loop.
973 : 36783 : std::vector< size_t >::iterator off_iter = vertAdjacencyOffsets.begin();
974 : 36783 : const std::vector< size_t >::iterator off_end = vertAdjacencyOffsets.end();
975 [ + - ]: 36783 : size_t prev = *off_iter;
976 [ + - ]: 36783 : ++off_iter;
977 [ + - ][ + - ]: 414178 : for( ; off_iter != off_end; ++off_iter )
[ + + ]
978 : : {
979 [ + - ]: 377395 : prev += *off_iter;
980 [ + - ]: 377395 : *off_iter = prev;
981 : : }
982 : :
983 : : // Allocate space for element numbers
984 [ + - ][ + - ]: 36783 : const size_t num_vert_uses = vertAdjacencyOffsets[num_nodes() - 1];
985 [ - + ]: 36783 : assert( num_vert_uses == elemConnectivityArray.size() );
986 [ + - ]: 36783 : vertAdjacencyArray.resize( num_vert_uses );
987 : :
988 : : // Fill vertAdjacencyArray, using the indices in vertAdjacencyOffsets
989 : : // as the location to insert the next element number in
990 : : // vertAdjacencyArray. When done, vertAdjacenyOffsets will contain
991 : : // the start index for each vertex, rather than one past the last
992 : : // index.
993 [ + + ]: 270463 : for( size_t i = 0; i < elementArray.size(); ++i )
994 : : {
995 [ + - ][ + - ]: 233680 : size_t* conn_iter = elementArray[i].get_vertex_index_array();
996 [ + - ][ + - ]: 233680 : const size_t* conn_end = conn_iter + elementArray[i].node_count();
997 [ + + ]: 1190628 : for( ; conn_iter != conn_end; ++conn_iter )
998 : : {
999 [ + - ]: 956948 : const size_t array_index = --vertAdjacencyOffsets[*conn_iter];
1000 [ + - ]: 956948 : vertAdjacencyArray[array_index] = i;
1001 : : }
1002 : : }
1003 : :
1004 : : // Last entry should be number of vertex uses (one past the
1005 : : // last index of the last vertex.)
1006 [ + - ][ + - ]: 36783 : vertAdjacencyOffsets[num_nodes()] = num_vert_uses;
1007 : : }
1008 : :
1009 : 0 : void PatchData::get_subpatch( size_t center_vertex_index, unsigned num_adj_elem_layers, PatchData& subpatch,
1010 : : MsqError& err )
1011 : : {
1012 : : // Make sure we're in range
1013 [ # # ][ # # ]: 0 : if( center_vertex_index >= num_free_vertices() )
1014 : : {
1015 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Invalid index for center vertex", MsqError::INVALID_ARG );
1016 : 0 : return;
1017 : : }
1018 : :
1019 : : // Notify any observers of the existing subpatch that the mesh
1020 : : // in the patch is to be changed.
1021 [ # # ]: 0 : subpatch.notify_new_patch();
1022 : :
1023 : : // Get list of vertices and elements in subpatch.
1024 : : // Ultimately, end up with arrays of unique, sorted indices.
1025 : : // It is important that the vertex indices be sorted so later
1026 : : // a reverse lookup can be done using a binary search (std::lower_bound).
1027 [ # # ][ # # ]: 0 : std::vector< size_t > elements, vertices, offsets;
[ # # ][ # # ]
[ # # ]
1028 [ # # ]: 0 : vertices.push_back( center_vertex_index );
1029 [ # # ]: 0 : for( unsigned i = 0; i < num_adj_elem_layers; ++i )
1030 : : {
1031 : 0 : elements.clear();
1032 [ # # ]: 0 : for( unsigned v = 0; v < vertices.size(); ++v )
1033 : : {
1034 : : size_t num_elem;
1035 [ # # ][ # # ]: 0 : const size_t* vert_elems = get_vertex_element_adjacencies( vertices[v], num_elem, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1036 [ # # ]: 0 : elements.insert( elements.end(), vert_elems, vert_elems + num_elem );
1037 : : }
1038 [ # # ]: 0 : std::sort( elements.begin(), elements.end() );
1039 [ # # ][ # # ]: 0 : elements.erase( std::unique( elements.begin(), elements.end() ), elements.end() );
1040 : :
1041 : 0 : vertices.clear();
1042 [ # # ]: 0 : for( unsigned e = 0; e < elements.size(); ++e )
1043 : : {
1044 [ # # ][ # # ]: 0 : MsqMeshEntity& elem = element_by_index( elements[e] );
1045 [ # # ]: 0 : size_t num_vert = elem.node_count();
1046 [ # # ]: 0 : const size_t* elem_verts = elem.get_vertex_index_array();
1047 [ # # ]: 0 : vertices.insert( vertices.end(), elem_verts, elem_verts + num_vert );
1048 : : }
1049 [ # # ]: 0 : std::sort( vertices.begin(), vertices.end() );
1050 [ # # ][ # # ]: 0 : vertices.erase( std::unique( vertices.begin(), vertices.end() ), vertices.end() );
1051 : : }
1052 : :
1053 : : // Allocate space for element connectivity info.
1054 : 0 : size_t num_vert_uses = 0;
1055 [ # # ]: 0 : for( unsigned i = 0; i < elements.size(); ++i )
1056 [ # # ][ # # ]: 0 : num_vert_uses += element_by_index( elements[i] ).node_count();
[ # # ]
1057 [ # # ]: 0 : subpatch.elementArray.resize( elements.size() );
1058 [ # # ]: 0 : subpatch.elementHandlesArray.resize( elements.size() );
1059 [ # # ]: 0 : subpatch.elemConnectivityArray.resize( num_vert_uses );
1060 [ # # ]: 0 : offsets.resize( elements.size() + 1 );
1061 : :
1062 : : // Construct element connectivity data in new patch,
1063 : : // and copy element type into new patch
1064 : 0 : size_t curr_offset = 0;
1065 [ # # ]: 0 : for( unsigned i = 0; i < elements.size(); ++i )
1066 : : {
1067 [ # # ][ # # ]: 0 : MsqMeshEntity& elem = element_by_index( elements[i] );
1068 [ # # ]: 0 : assert( i < elementArray.size() );
1069 [ # # ][ # # ]: 0 : subpatch.elementArray[i].set_element_type( elem.get_element_type() );
[ # # ]
1070 [ # # ][ # # ]: 0 : subpatch.elementHandlesArray[i] = elementHandlesArray[elements[i]];
[ # # ]
1071 [ # # ]: 0 : const size_t* verts = elem.get_vertex_index_array();
1072 [ # # ]: 0 : offsets[i] = curr_offset;
1073 [ # # ][ # # ]: 0 : for( unsigned j = 0; j < elem.node_count(); ++j )
1074 : : {
1075 [ # # ]: 0 : subpatch.elemConnectivityArray[curr_offset++] =
1076 [ # # ][ # # ]: 0 : std::lower_bound( vertices.begin(), vertices.end(), verts[j] ) - vertices.begin();
1077 : : }
1078 : : }
1079 [ # # ]: 0 : offsets[elements.size()] = curr_offset;
1080 : :
1081 : : // Store index in this patch in vertex handle array of subpatch
1082 : : // so we can determine how vertices were reordered when setting
1083 : : // vertex coordinates.
1084 : : assert( sizeof( size_t ) == sizeof( void* ) );
1085 [ # # ]: 0 : subpatch.vertexHandlesArray.resize( vertices.size() );
1086 [ # # ]: 0 : size_t* vert_handles = reinterpret_cast< size_t* >( &subpatch.vertexHandlesArray[0] );
1087 [ # # ]: 0 : std::copy( vertices.begin(), vertices.end(), vert_handles );
1088 : :
1089 : : // All vertices except vertex at center_vertex_index are fixed.
1090 [ # # ]: 0 : subpatch.byteArray.resize( vertices.size() );
1091 [ # # ]: 0 : for( size_t pi = 0; pi < vertices.size(); ++pi )
1092 : : {
1093 [ # # ][ # # ]: 0 : if( vertices[pi] == center_vertex_index )
1094 [ # # ][ # # ]: 0 : subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() & ~MsqVertex::MSQ_PATCH_FIXED;
[ # # ][ # # ]
1095 : : else
1096 [ # # ][ # # ]: 0 : subpatch.byteArray[pi] = vertexArray[vertices[pi]].get_flags() | MsqVertex::MSQ_PATCH_FIXED;
[ # # ][ # # ]
1097 : : }
1098 : :
1099 : : // Re-order vertices and initialize other data in subpatch
1100 [ # # ][ # # ]: 0 : subpatch.initialize_data( arrptr( offsets ), &subpatch.byteArray[0], err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1101 : :
1102 : : // Copy vertex data into subpatch. subpatch.vertexHandlesArray contains
1103 : : // the indices into this PatchData for each vertex, as reordered by the
1104 : : // call to initialize_data.
1105 [ # # ]: 0 : subpatch.vertexArray.resize( vertices.size() );
1106 [ # # ]: 0 : for( unsigned i = 0; i < vertices.size(); ++i )
1107 : : {
1108 [ # # ]: 0 : size_t vert_index = ( size_t )( subpatch.vertexHandlesArray[i] );
1109 [ # # ]: 0 : vertices[i] = vert_index;
1110 [ # # ][ # # ]: 0 : subpatch.vertexHandlesArray[i] = vertexHandlesArray[vert_index];
1111 [ # # ][ # # ]: 0 : subpatch.vertexArray[i] = vertexArray[vert_index];
[ # # ]
1112 [ # # ][ # # ]: 0 : subpatch.vertexArray[i].flags() = subpatch.byteArray[i];
[ # # ]
1113 : : }
1114 : :
1115 : 0 : subpatch.myMesh = myMesh;
1116 : 0 : subpatch.myDomain = myDomain;
1117 : 0 : subpatch.mSettings = mSettings;
1118 : :
1119 [ # # ][ # # ]: 0 : notify_sub_patch( subpatch, arrptr( vertices ), elements.empty() ? 0 : arrptr( elements ), err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1120 : : }
1121 : :
1122 : : //! Adjust the position of the specified vertex so that it
1123 : : //! lies on its constraining domain. The actual domain constraint
1124 : : //! is managed by the TSTT mesh implementation
1125 : 2734298 : void PatchData::snap_vertex_to_domain( size_t vertex_index, MsqError& err )
1126 : : {
1127 [ + + ]: 2734298 : if( domain_set() )
1128 : : {
1129 : : // if not doing normal caching or vertex is not on a single surface
1130 [ + + ]: 1908516 : if( normalData.empty() )
1131 : 1420726 : { get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] ); }
1132 : : // entire domain is 2-D (all vertices have a single normal)
1133 [ + + ]: 487790 : else if( vertexNormalIndices.empty() )
1134 : : {
1135 [ + - ]: 1462746 : get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
1136 [ + - ][ - + ]: 1462746 : vertexArray[vertex_index], normalData[vertex_index], err );MSQ_ERRRTN( err );
[ # # ][ - + ]
1137 : : }
1138 [ + + ]: 208 : else if( vertexNormalIndices[vertex_index] < normalData.size() )
1139 : : { // vertex has a unique normal
1140 [ + - ]: 129 : get_domain()->closest_point( vertexHandlesArray[vertex_index], Vector3D( vertexArray[vertex_index] ),
1141 : 43 : vertexArray[vertex_index], normalData[vertexNormalIndices[vertex_index]],
1142 [ + - ][ - + ]: 86 : err );MSQ_ERRRTN( err );
[ # # ][ - + ]
1143 : : }
1144 : : else
1145 : : { // vertex has multiple normals
1146 : 2734298 : get_domain()->snap_to( vertexHandlesArray[vertex_index], vertexArray[vertex_index] );
1147 : : }
1148 : : }
1149 : : }
1150 : :
1151 : 118608 : void PatchData::update_cached_normals( MsqError& err )
1152 : : {
1153 : : size_t i;
1154 : :
1155 [ + - ]: 118608 : MeshDomain* domain = get_domain();
1156 [ - + ]: 118608 : if( !domain )
1157 : : {
1158 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
1159 : 118608 : return;
1160 : : }
1161 : :
1162 : : // Determine which vertices lie on surfaces
1163 [ + - ][ + - ]: 118608 : vertexDomainDOF.resize( num_nodes() );
1164 [ + - ][ + - ]: 118608 : domain->domain_DoF( arrptr( vertexHandlesArray ), arrptr( vertexDomainDOF ), num_nodes(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
1165 : :
1166 : : // Count how many vertices have a single normal
1167 : : // Sun doesn't support partial template specialization, so can't use std::count
1168 : : // size_t n = std::count( vertexDomainDOF.begin(), vertexDomainDOF.end(), 2 );
1169 : 118608 : size_t n = 0;
1170 : 118608 : std::vector< unsigned short >::iterator k;
1171 [ + - ][ + - ]: 918543 : for( k = vertexDomainDOF.begin(); k != vertexDomainDOF.end(); ++k )
[ + + ]
1172 [ + - ][ + + ]: 799935 : if( *k == 2 ) ++n;
1173 : :
1174 [ + - ]: 118608 : normalData.resize( n );
1175 : :
1176 : : // If all vertices are on a surface, pass in the existing handles array
1177 : : // and store a single normal per vertex.
1178 [ + - ][ + + ]: 118608 : if( n == num_nodes() )
1179 : : {
1180 [ + - ]: 118537 : std::copy( vertexArray.begin(), vertexArray.end(), normalData.begin() );
1181 [ + - ][ + - ]: 118537 : domain->vertex_normal_at( arrptr( vertexHandlesArray ), arrptr( normalData ), num_nodes(), err );
[ + - ][ + - ]
1182 : 118537 : vertexNormalIndices.clear();
1183 [ + - ][ - + ]: 118537 : vertexDomainDOF.clear();MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
1184 : : }
1185 : : else
1186 : : {
1187 [ + - ][ + - ]: 71 : vertexNormalIndices.resize( num_nodes() );
1188 : 71 : size_t nn = 0;
1189 [ + - ][ + + ]: 498 : for( i = 0; i < num_nodes(); ++i )
1190 : : {
1191 [ + - ][ + + ]: 427 : if( vertexDomainDOF[i] == 2 )
1192 : : {
1193 [ + - ][ + - ]: 71 : normalData[nn] = vertexArray[i];
[ + - ]
1194 [ + - ][ + - ]: 71 : domain->vertex_normal_at( vertexHandlesArray[i], normalData[nn] );
[ + - ]
1195 [ + - ]: 71 : vertexNormalIndices[i] = nn;
1196 : 71 : ++nn;
1197 : : }
1198 : : else
1199 : : {
1200 [ + - ]: 356 : vertexNormalIndices[i] = n + 1;
1201 : : }
1202 : : }
1203 [ - + ]: 118608 : assert( nn == n );
1204 : : }
1205 : : }
1206 : :
1207 : 335285 : void PatchData::get_domain_normal_at_element( size_t elem_index, Vector3D& surf_norm, MsqError& err )
1208 : : {
1209 : : // check if element as a mid-face node
1210 : 335285 : const MsqMeshEntity& elem = element_by_index( elem_index );
1211 [ - + ][ # # ]: 335285 : const int mid_node = TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 2, 0, err );MSQ_ERRRTN( err );
[ - + ]
1212 : : // if we have the mid-element vertex, get cached normal for it
1213 [ - + ]: 335285 : if( mid_node > 0 )
1214 : : {
1215 : 0 : get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index],
1216 [ # # ][ # # ]: 0 : surf_norm, err );MSQ_ERRRTN( err );
[ # # ]
1217 : : }
1218 : : // otherwise query domain for normal at element centroid
1219 [ + - ]: 335285 : else if( domain_set() )
1220 : : {
1221 [ - + ]: 335285 : assert( elem_index < elementArray.size() );
1222 [ - + ][ # # ]: 335285 : elementArray[elem_index].get_centroid( surf_norm, *this, err );MSQ_ERRRTN( err );
[ - + ]
1223 : 335285 : get_domain()->element_normal_at( elementHandlesArray[elem_index], surf_norm );
1224 : : }
1225 : : else
1226 [ # # ]: 335285 : MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
1227 : : }
1228 : :
1229 : 33686 : void PatchData::get_domain_normal_at_mid_edge( size_t elem_index, unsigned edge_num, Vector3D& normal, MsqError& err )
1230 : : {
1231 : : // check if element as a mid-edge node
1232 : 33686 : const MsqMeshEntity& elem = element_by_index( elem_index );
1233 : : const int mid_node =
1234 [ - + ][ # # ]: 33686 : TopologyInfo::higher_order_from_side( elem.get_element_type(), elem.node_count(), 1, edge_num, err );MSQ_ERRRTN( err );
[ - + ]
1235 : : // if we have the mid-edge vertex, get cached normal for it
1236 [ + - ]: 33686 : if( mid_node > 0 )
1237 : : {
1238 : 33686 : get_domain_normal_at_vertex( elem.get_vertex_index_array()[mid_node], elementHandlesArray[elem_index], normal,
1239 [ - + ][ # # ]: 33686 : err );MSQ_ERRRTN( err );
[ - + ]
1240 : : }
1241 : : // otherwise query domain for normal at center of edge
1242 [ # # ]: 0 : else if( domain_set() )
1243 : : {
1244 [ # # ][ # # ]: 0 : const unsigned* edge = TopologyInfo::edge_vertices( elem.get_element_type(), edge_num, err );MSQ_ERRRTN( err );
[ # # ]
1245 : 0 : const MsqVertex& v1 = vertex_by_index( elem.get_vertex_index_array()[edge[0]] );
1246 : 0 : const MsqVertex& v2 = vertex_by_index( elem.get_vertex_index_array()[edge[1]] );
1247 [ # # ][ # # ]: 0 : normal = 0.5 * ( v1 + v2 );
1248 : 0 : get_domain()->element_normal_at( elementHandlesArray[elem_index], normal );
1249 : : }
1250 : : else
1251 : : {
1252 [ # # ]: 0 : MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
1253 : 33686 : return;
1254 : : }
1255 : : }
1256 : :
1257 : 0 : void PatchData::get_domain_normals_at_corners( size_t elem_index, Vector3D normals_out[], MsqError& err )
1258 : : {
1259 [ # # ]: 0 : if( !domain_set() )
1260 : : {
1261 [ # # ]: 0 : MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
1262 : 0 : return;
1263 : : }
1264 : :
1265 [ # # ]: 0 : assert( elem_index < elementArray.size() );
1266 [ # # ]: 0 : if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
1267 : : {
1268 [ # # ]: 0 : MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
1269 : 0 : return;
1270 : : }
1271 : :
1272 [ # # ]: 0 : if( normalData.empty() )
1273 : : {
1274 [ # # ][ # # ]: 0 : update_cached_normals( err );MSQ_ERRRTN( err );
[ # # ]
1275 : : }
1276 : :
1277 : 0 : MsqMeshEntity& elem = elementArray[elem_index];
1278 : 0 : const unsigned count = elem.vertex_count();
1279 : 0 : const size_t* const vertex_indices = elem.get_vertex_index_array();
1280 [ # # ]: 0 : for( size_t i = 0; i < count; ++i )
1281 : : {
1282 : 0 : const size_t v = vertex_indices[i];
1283 [ # # ]: 0 : if( vertexNormalIndices.empty() ) { normals_out[i] = normalData[v]; }
1284 [ # # ]: 0 : else if( vertexNormalIndices[v] < normalData.size() )
1285 : : {
1286 : 0 : normals_out[i] = normalData[vertexNormalIndices[v]];
1287 : : }
1288 : : else
1289 : : {
1290 : 0 : normals_out[i] = vertexArray[v];
1291 : 0 : get_domain()->element_normal_at( elementHandlesArray[elem_index], normals_out[i] );
1292 : : }
1293 : : }
1294 : : }
1295 : :
1296 : 12713940 : void PatchData::get_domain_normal_at_vertex( size_t vert_index, Mesh::EntityHandle handle, Vector3D& normal,
1297 : : MsqError& err )
1298 : : {
1299 [ - + ]: 12713940 : if( !domain_set() )
1300 : : {
1301 [ # # ]: 0 : MSQ_SETERR( err )( "No domain constraint set.", MsqError::INVALID_STATE );
1302 : 0 : return;
1303 : : }
1304 : :
1305 [ + + ]: 12713940 : if( normalData.empty() )
1306 : : {
1307 [ - + ][ # # ]: 118608 : update_cached_normals( err );MSQ_ERRRTN( err );
[ - + ]
1308 : : }
1309 : :
1310 [ + + ]: 12713940 : if( vertexNormalIndices.empty() ) { normal = normalData[vert_index]; }
1311 [ + + ]: 2621 : else if( vertexNormalIndices[vert_index] < normalData.size() )
1312 : : {
1313 : 657 : normal = normalData[vertexNormalIndices[vert_index]];
1314 : : }
1315 : : else
1316 : : {
1317 : 1964 : normal = vertexArray[vert_index];
1318 : 12713940 : get_domain()->element_normal_at( handle, normal );
1319 : : }
1320 : : }
1321 : :
1322 : 12680254 : void PatchData::get_domain_normal_at_corner( size_t elem_index, unsigned corner, Vector3D& normal, MsqError& err )
1323 : : {
1324 [ - + ]: 12680254 : assert( elem_index < elementArray.size() );
1325 [ - + ]: 12680254 : if( 2 != TopologyInfo::dimension( elementArray[elem_index].get_element_type() ) )
1326 : : {
1327 [ # # ]: 0 : MSQ_SETERR( err )( "Attempt to get corners of non-surface element", MsqError::INVALID_ARG );
1328 : 12680254 : return;
1329 : : }
1330 : :
1331 : 12680254 : MsqMeshEntity& elem = elementArray[elem_index];
1332 : 12680254 : const size_t* const vertex_indices = elem.get_vertex_index_array();
1333 [ - + ][ # # ]: 12680254 : get_domain_normal_at_vertex( vertex_indices[corner], elementHandlesArray[elem_index], normal, err );MSQ_CHKERR( err );
1334 : : }
1335 : :
1336 : 328 : void PatchData::set_mesh( Mesh* ms )
1337 : : {
1338 : 328 : myMesh = ms;
1339 : : // observers should treat this the same as if the
1340 : : // instance of this object wzs being deleted.
1341 : 328 : notify_patch_destroyed();
1342 : 328 : }
1343 : :
1344 : 280 : void PatchData::set_domain( MeshDomain* d )
1345 : : {
1346 : 280 : myDomain = d;
1347 : :
1348 : : // clear all cached data from the previous domain
1349 : 280 : vertexNormalIndices.clear();
1350 : 280 : normalData.clear();
1351 : : // vertexDomainDOF.clear();
1352 : :
1353 : : // observers should treat this the same as if the
1354 : : // instance of this object wzs being deleted.
1355 : 280 : notify_patch_destroyed();
1356 : 280 : }
1357 : :
1358 : 0 : static int width( double d )
1359 : : {
1360 [ # # ]: 0 : if( d == 0.0 ) return 1;
1361 : 0 : const int max_precision = 6;
1362 : 0 : int w = (int)ceil( log10( 0.001 + fabs( d ) ) );
1363 [ # # ][ # # ]: 0 : if( w < 0 ) w = 2 + std::min( max_precision, -w );
1364 [ # # ]: 0 : if( d < 0.0 ) ++w;
1365 : 0 : return w;
1366 : : }
1367 : 0 : static int width( size_t t )
1368 : : {
1369 [ # # ]: 0 : return t ? (int)ceil( log10( (double)( 1 + t ) ) ) : 1;
1370 : : }
1371 : 0 : static int width( const void* ptr )
1372 : : {
1373 : 0 : return width( (size_t)ptr );
1374 : : }
1375 : :
1376 : 0 : ostream& operator<<( ostream& stream, const PatchData& pd )
1377 : : {
1378 : : size_t i;
1379 : 0 : int fw = 5; // width of bit flags
1380 : 0 : int hw = 6; // width of a handle
1381 : 0 : int cw = 4; // with of a coordinate value
1382 : 0 : int iw = 3; // with of an index
1383 : 0 : int tw = 3; // with of the string name of an element type
1384 : 0 : int xw = cw, yw = cw, zw = cw;
1385 : :
1386 [ # # ]: 0 : for( i = 0; i < pd.num_nodes(); ++i )
1387 : : {
1388 : 0 : int w = 2 + width( pd.vertexHandlesArray[i] );
1389 [ # # ]: 0 : if( w > hw ) hw = w;
1390 : 0 : w = width( pd.vertexArray[i].x() );
1391 [ # # ]: 0 : if( w > xw ) xw = w;
1392 : 0 : w = width( pd.vertexArray[i].y() );
1393 [ # # ]: 0 : if( w > yw ) yw = w;
1394 : 0 : w = width( pd.vertexArray[i].z() );
1395 [ # # ]: 0 : if( w > zw ) zw = w;
1396 : : }
1397 [ # # ]: 0 : for( i = 0; i < pd.num_elements(); ++i )
1398 : : {
1399 : 0 : int w = 2 + width( pd.elementHandlesArray[i] );
1400 [ # # ]: 0 : if( w > hw ) hw = w;
1401 : 0 : const char* name = TopologyInfo::short_name( pd.elementArray[i].get_element_type() );
1402 [ # # ][ # # ]: 0 : if( name && (int)strlen( name ) > tw ) tw = strlen( name );
1403 : : }
1404 [ # # ]: 0 : if( iw < (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) ) )
1405 : 0 : iw = (int)ceil( log10( (double)( 1 + pd.num_nodes() ) ) );
1406 [ # # ]: 0 : if( iw < (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) ) )
1407 : 0 : iw = (int)ceil( log10( (double)( 1 + pd.num_elements() ) ) );
1408 : :
1409 : 0 : stream << "Vertices: " << endl;
1410 : 0 : stream << "Flags: C: culled, F: fixed, S: slave, P: patch vertex, M: marked" << endl;
1411 : 0 : stream << setw( iw ) << "Idx"
1412 : 0 : << " " << setw( hw ) << "Handle"
1413 : 0 : << " " << setw( cw ) << "X"
1414 : 0 : << "," << setw( cw ) << "Y"
1415 : 0 : << "," << setw( cw ) << "Z"
1416 : 0 : << " " << setw( fw ) << "Flags"
1417 : 0 : << " "
1418 : 0 : << "Adj.Elems" << endl
1419 : 0 : << setw( iw ) << setfill( '-' ) << ""
1420 : 0 : << " " << setw( hw ) << setfill( '-' ) << ""
1421 : 0 : << " " << setw( cw ) << setfill( '-' ) << ""
1422 : 0 : << "," << setw( cw ) << setfill( '-' ) << ""
1423 : 0 : << "," << setw( cw ) << setfill( '-' ) << ""
1424 : 0 : << " " << setw( fw ) << setfill( '-' ) << ""
1425 : 0 : << " " << setfill( ' ' ) << "-------------" << std::endl;
1426 [ # # ]: 0 : for( i = 0; i < pd.num_nodes(); ++i )
1427 : : {
1428 : 0 : stream << setw( iw ) << i << " " << setw( hw ) << pd.vertexHandlesArray[i] << " " << setw( cw )
1429 : 0 : << pd.vertexArray[i].x() << "," << setw( cw ) << pd.vertexArray[i].y() << "," << setw( cw )
1430 : 0 : << pd.vertexArray[i].z() << " ";
1431 [ # # ]: 0 : if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_CULLED ) )
1432 : 0 : stream << "C";
1433 : : else
1434 : 0 : stream << " ";
1435 [ # # ]: 0 : if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_HARD_FIXED ) )
1436 : 0 : stream << "F";
1437 : : else
1438 : 0 : stream << " ";
1439 [ # # ]: 0 : if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_DEPENDENT ) )
1440 : 0 : stream << "S";
1441 : : else
1442 : 0 : stream << " ";
1443 [ # # ]: 0 : if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_PATCH_FIXED ) )
1444 : 0 : stream << "f";
1445 : : else
1446 : 0 : stream << " ";
1447 [ # # ]: 0 : if( pd.vertexArray[i].is_flag_set( MsqVertex::MSQ_MARK ) )
1448 : 0 : stream << "M";
1449 : : else
1450 : 0 : stream << " ";
1451 : :
1452 [ # # ]: 0 : if( pd.vertAdjacencyArray.size() )
1453 : : {
1454 : 0 : size_t j = pd.vertAdjacencyOffsets[i];
1455 : 0 : size_t end = pd.vertAdjacencyOffsets[i + 1];
1456 [ # # ]: 0 : if( j != end ) stream << " " << pd.vertAdjacencyArray[j++];
1457 [ # # ]: 0 : for( ; j < end; ++j )
1458 : 0 : stream << "," << pd.vertAdjacencyArray[j];
1459 : : }
1460 : :
1461 : 0 : stream << endl;
1462 : : }
1463 : :
1464 : 0 : stream << "Elements: " << endl;
1465 : 0 : stream << setw( iw ) << "Idx"
1466 : 0 : << " " << setw( hw ) << "Handle"
1467 : 0 : << " " << setw( tw + 2 ) << "Type"
1468 : 0 : << " "
1469 : 0 : << "Connectivity" << std::endl
1470 : 0 : << setw( iw ) << setfill( '-' ) << ""
1471 : 0 : << " " << setw( hw ) << setfill( '-' ) << ""
1472 : 0 : << " " << setw( tw + 2 ) << setfill( '-' ) << ""
1473 : 0 : << " " << setfill( ' ' ) << "--------------------------" << std::endl;
1474 [ # # ]: 0 : for( i = 0; i < pd.num_elements(); ++i )
1475 : : {
1476 : 0 : EntityTopology type = pd.elementArray[i].get_element_type();
1477 : 0 : stream << setw( iw ) << i << " " << setw( hw ) << pd.elementHandlesArray[i] << " " << setw( tw )
1478 : 0 : << TopologyInfo::short_name( type ) << left << setw( 2 ) << pd.elementArray[i].node_count() << internal
1479 : 0 : << " " << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[0];
1480 [ # # ]: 0 : for( size_t j = 1; j < pd.elementArray[i].node_count(); ++j )
1481 : 0 : stream << "," << setw( iw ) << pd.elementArray[i].get_vertex_index_array()[j];
1482 : 0 : stream << endl;
1483 : : }
1484 : 0 : stream << endl;
1485 : :
1486 [ # # ]: 0 : stream << "Mesh: " << ( pd.myMesh ? "yes" : "no" ) << endl;
1487 [ # # ]: 0 : stream << "Domain: " << ( pd.myDomain ? "yes" : "no" ) << endl;
1488 : : // stream << "mType: " << (pd.mType==PatchData::VERTICES_ON_VERTEX_PATCH?"vert-on-vert":
1489 : : // pd.mType==PatchData::ELEMENTS_ON_VERTEX_PATCH?"elem-on-vert":
1490 : : // pd.mType==PatchData::GLOBAL_PATCH?"global":"unknown") << endl;
1491 : :
1492 [ # # ]: 0 : if( pd.haveComputedInfos )
1493 : : {
1494 : 0 : stream << "ComputedInfos:" << endl;
1495 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MIN_UNSIGNED_AREA ) )
1496 : 0 : stream << "\t MIN_UNSINGED_AREA = " << pd.computedInfos[PatchData::MIN_UNSIGNED_AREA] << endl;
1497 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MAX_UNSIGNED_AREA ) )
1498 : 0 : stream << "\t MAX_UNSIGNED_AREA = " << pd.computedInfos[PatchData::MAX_UNSIGNED_AREA] << endl;
1499 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MIN_EDGE_LENGTH ) )
1500 : 0 : stream << "\t MIN_EDGE_LENGTH = " << pd.computedInfos[PatchData::MIN_EDGE_LENGTH] << endl;
1501 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MAX_EDGE_LENGTH ) )
1502 : 0 : stream << "\t MAX_EDGE_LENGTH = " << pd.computedInfos[PatchData::MAX_EDGE_LENGTH] << endl;
1503 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET2D ) )
1504 : 0 : stream << "\t MINMAX_SIGNED_DET2D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET2D] << endl;
1505 [ # # ]: 0 : if( pd.have_computed_info( PatchData::MINMAX_SIGNED_DET3D ) )
1506 : 0 : stream << "\t MINMAX_SIGNED_DET3D = " << pd.computedInfos[PatchData::MINMAX_SIGNED_DET3D] << endl;
1507 [ # # ]: 0 : if( pd.have_computed_info( PatchData::AVERAGE_DET3D ) )
1508 : 0 : stream << "\t AVERAGE_DET3D = " << pd.computedInfos[PatchData::AVERAGE_DET3D] << endl;
1509 : : }
1510 : :
1511 : 0 : return stream << endl;
1512 : : }
1513 : :
1514 : 0 : void print_patch_data( const PatchData& pd )
1515 : : {
1516 : 0 : std::cout << pd << std::endl;
1517 : 0 : }
1518 : :
1519 : 0 : void PatchData::enslave_higher_order_nodes( const size_t* elem_offset_array, unsigned char* vertex_flags,
1520 : : MsqError& ) const
1521 : : {
1522 [ # # ]: 0 : for( size_t i = 0; i < elementArray.size(); ++i )
1523 : : {
1524 : 0 : size_t start = elem_offset_array[i];
1525 : 0 : size_t conn_len = elem_offset_array[i + 1] - start;
1526 [ # # ]: 0 : for( size_t j = elementArray[i].vertex_count(); j < conn_len; ++j )
1527 : : {
1528 : 0 : const size_t vert_idx = elemConnectivityArray[start + j];
1529 [ # # ]: 0 : assert( vert_idx < vertexHandlesArray.size() );
1530 [ # # ]: 0 : if( !( vertex_flags[vert_idx] & MsqVertex::MSQ_HARD_FIXED ) )
1531 : 0 : vertex_flags[vert_idx] |= MsqVertex::MSQ_DEPENDENT;
1532 : : }
1533 : : }
1534 : 0 : }
1535 : :
1536 : 1128315 : void PatchData::initialize_data( size_t* elem_offset_array, unsigned char* vertex_flags, MsqError& )
1537 : : {
1538 : : // Clear out data specific to patch
1539 : 1128315 : vertexNormalIndices.clear();
1540 : 1128315 : normalData.clear();
1541 : : // vertexDomainDOF.clear();
1542 : :
1543 : : // Clear any vertex->element adjacency data. It
1544 : : // is probably invalid, and certainly will be by the time
1545 : : // this function completes if the mesh contains higher-order
1546 : : // elements.
1547 : 1128315 : vertAdjacencyArray.clear();
1548 : 1128315 : vertAdjacencyOffsets.clear();
1549 : : size_t i, j;
1550 [ + + ]: 5069195 : for( i = 0; i < elementArray.size(); ++i )
1551 : : {
1552 : 3940880 : size_t start = elem_offset_array[i];
1553 : 3940880 : size_t conn_len = elem_offset_array[i + 1] - start;
1554 [ - + ]: 3940880 : assert( conn_len > 0 );
1555 : 3940880 : elementArray[i].set_connectivity( &elemConnectivityArray[start], conn_len );
1556 : : }
1557 : :
1558 : : // Use vertAdjacencyOffsets array as temporary storage.
1559 : 1128315 : vertAdjacencyOffsets.resize( vertexHandlesArray.size() + 1 );
1560 : 1128315 : size_t* vertex_index_map = arrptr( vertAdjacencyOffsets );
1561 : :
1562 : : // Count number of free vertices and initialize vertex_index_map
1563 : 1128315 : numFreeVertices = 0;
1564 [ + + ]: 10156385 : for( i = 0; i < vertexHandlesArray.size(); ++i )
1565 : : {
1566 [ + + ]: 9028070 : if( !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) ) ++numFreeVertices;
1567 : 9028070 : vertex_index_map[i] = i;
1568 : : }
1569 : :
1570 : : // Re-order vertices such that all free vertices are
1571 : : // first in the list. Construct map from old to new
1572 : : // position in list for updating element connectivity.
1573 : 1128315 : i = 0;
1574 : 1128315 : j = numFreeVertices;
1575 : 989291 : for( ;; ++i, ++j )
1576 : : {
1577 : : // find next fixed vertex in the range [i,numFreeVertices)
1578 [ + + ][ + + ]: 2944869 : for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ); ++i )
1579 : : ;
1580 : : // if no more fixed vertices in the free vertex range [0,numFreeVertices)
1581 [ + + ]: 2117606 : if( i == numFreeVertices ) break;
1582 : : // find the next free vertex in the range [j,num_nodes)
1583 [ + + ]: 3664171 : for( ; vertex_flags[j] & MsqVertex::MSQ_FIXED; ++j )
1584 : : ;
1585 : : // swap fixed (i) and free (j) vertices
1586 : 989291 : vertex_index_map[i] = j;
1587 : 989291 : vertex_index_map[j] = i;
1588 : 989291 : std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
1589 : 989291 : std::swap( vertex_flags[i], vertex_flags[j] );
1590 : : }
1591 [ - + ]: 1128315 : assert( i == numFreeVertices );
1592 [ - + ]: 1128315 : assert( j <= vertexHandlesArray.size() );
1593 : :
1594 : : // Update element connectivity data for new vertex indices
1595 [ + + ]: 16993745 : for( i = 0; i < elemConnectivityArray.size(); ++i )
1596 : 15865430 : elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
1597 : :
1598 : : // Reorder vertices such that free, slave vertices
1599 : : // occur after free, non-slave vertices in list.
1600 : 1128315 : numSlaveVertices = 0;
1601 [ + + ]: 10156385 : for( i = 0; i < vertexHandlesArray.size(); ++i )
1602 : : {
1603 [ + + ][ + + ]: 9028070 : if( ( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ) && !( vertex_flags[i] & MsqVertex::MSQ_FIXED ) )
1604 : 366 : ++numSlaveVertices;
1605 : : }
1606 : 1128315 : numFreeVertices -= numSlaveVertices;
1607 : :
1608 [ + + ]: 1128418 : if( numSlaveVertices )
1609 : : {
1610 : : // re-initialize vertex index map
1611 [ + + ]: 1172 : for( i = 0; i < vertexHandlesArray.size(); ++i )
1612 : 1080 : vertex_index_map[i] = i;
1613 : :
1614 : : // Re-order free vertices such that all slave vertices are
1615 : : // last in the list. Construct map from old to new
1616 : : // position in list for updating element connectivity.
1617 : 92 : i = 0;
1618 : 92 : j = numFreeVertices;
1619 : 103 : for( ;; ++i, ++j )
1620 : : {
1621 : : // find next slave vertex in the range [i,numFreeVertices)
1622 [ + + ][ + + ]: 522 : for( ; i < numFreeVertices && !( vertex_flags[i] & MsqVertex::MSQ_DEPENDENT ); ++i )
1623 : : ;
1624 : : // if no more slave vertices in [0,numFreeVertices), then done.
1625 [ + + ]: 195 : if( i == numFreeVertices ) break;
1626 : : // find the next free (non-slave) vertex in the range
1627 : : // [numFreeVertices,numFreeVertices+numSlaveVertices)
1628 [ + + ]: 137 : for( ; vertex_flags[j] & MsqVertex::MSQ_DEPENDENT; ++j )
1629 : : ;
1630 : : // swap free (j) and slave (i) vertices
1631 : 103 : vertex_index_map[i] = j;
1632 : 103 : vertex_index_map[j] = i;
1633 : 103 : std::swap( vertexHandlesArray[i], vertexHandlesArray[j] );
1634 : 103 : std::swap( vertex_flags[i], vertex_flags[j] );
1635 : : }
1636 [ - + ]: 92 : assert( i == numFreeVertices );
1637 [ - + ]: 92 : assert( j <= numFreeVertices + numSlaveVertices );
1638 : :
1639 : : // Update element connectivity data for new vertex indices
1640 [ + + ]: 1952 : for( i = 0; i < elemConnectivityArray.size(); ++i )
1641 : 1860 : elemConnectivityArray[i] = vertex_index_map[elemConnectivityArray[i]];
1642 : : }
1643 : :
1644 : : // Clear temporary data
1645 : 1128315 : vertAdjacencyOffsets.clear();
1646 : :
1647 : 2117606 : notify_new_patch();
1648 : 1128315 : }
1649 : :
1650 : 0 : size_t PatchData::num_corners() const
1651 : : {
1652 : 0 : size_t result = 0;
1653 [ # # ]: 0 : for( unsigned i = 0; i < elementArray.size(); ++i )
1654 : 0 : result += elementArray[i].vertex_count();
1655 : 0 : return result;
1656 : : }
1657 : :
1658 : 0 : void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, EntityTopology type,
1659 : : const size_t* connectivity, const bool* fixed, MsqError& err )
1660 : : {
1661 [ # # ]: 0 : std::vector< EntityTopology > types( num_elem );
1662 [ # # ]: 0 : std::fill( types.begin(), types.end(), type );
1663 [ # # ][ # # ]: 0 : const EntityTopology* type_ptr = num_elem ? arrptr( types ) : 0;
1664 [ # # ][ # # ]: 0 : this->fill( num_vertex, coords, num_elem, type_ptr, connectivity, fixed, err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ]
1665 : 0 : }
1666 : :
1667 : 0 : void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
1668 : : const size_t* conn, const bool* fixed, MsqError& err )
1669 : : {
1670 [ # # ]: 0 : std::vector< size_t > lengths( num_elem );
1671 [ # # ][ # # ]: 0 : std::transform( types, types + num_elem, lengths.begin(), std::ptr_fun( TopologyInfo::corners ) );
1672 [ # # ][ # # ]: 0 : const size_t* len_ptr = num_elem ? arrptr( lengths ) : 0;
1673 [ # # ][ # # ]: 0 : this->fill( num_vertex, coords, num_elem, types, len_ptr, conn, fixed, err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ]
1674 : 0 : }
1675 : :
1676 : 0 : void PatchData::fill( size_t num_vertex, const double* coords, size_t num_elem, const EntityTopology* types,
1677 : : const size_t* lengths, const size_t* conn, const bool* fixed, MsqError& err )
1678 : : {
1679 : : size_t i;
1680 : :
1681 : : // count vertex uses
1682 [ # # ]: 0 : size_t num_uses = std::accumulate( lengths, lengths + num_elem, 0 );
1683 : :
1684 : : // Allocate storage for data
1685 [ # # ]: 0 : vertexArray.resize( num_vertex );
1686 [ # # ]: 0 : vertexHandlesArray.resize( num_vertex );
1687 [ # # ]: 0 : elementArray.resize( num_elem );
1688 [ # # ]: 0 : elementHandlesArray.resize( num_elem );
1689 [ # # ]: 0 : elemConnectivityArray.resize( num_uses );
1690 : 0 : numFreeVertices = 0;
1691 : 0 : numSlaveVertices = 0;
1692 : :
1693 : : // Must call clear() first so that any stale values get
1694 : : // zero'd when we call resize.
1695 : 0 : byteArray.clear();
1696 [ # # ]: 0 : if( fixed )
1697 : : {
1698 [ # # ]: 0 : byteArray.resize( num_vertex, 0 );
1699 [ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
1700 [ # # ][ # # ]: 0 : if( fixed[i] ) byteArray[i] |= ( MsqVertex::MSQ_HARD_FIXED | MsqVertex::MSQ_PATCH_FIXED );
1701 : : }
1702 : :
1703 [ # # ]: 0 : for( i = 0; i < num_elem; ++i )
1704 : : {
1705 [ # # ][ # # ]: 0 : element_by_index( i ).set_element_type( types[i] );
1706 [ # # ]: 0 : elementHandlesArray[i] = (Mesh::ElementHandle)i;
1707 : : }
1708 [ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
1709 [ # # ]: 0 : vertexHandlesArray[i] = (Mesh::VertexHandle)i;
1710 : :
1711 [ # # ]: 0 : memcpy( get_connectivity_array(), conn, num_uses * sizeof( size_t ) );
1712 : :
1713 [ # # ]: 0 : std::vector< size_t > offsets( num_elem + 1 );
1714 [ # # ]: 0 : size_t sum = offsets[0] = 0;
1715 [ # # ]: 0 : for( i = 1; i <= num_elem; ++i )
1716 [ # # ]: 0 : offsets[i] = sum += lengths[i - 1];
1717 : :
1718 : : const Settings::HigherOrderSlaveMode ho_mode =
1719 [ # # ][ # # ]: 0 : mSettings ? mSettings->get_slaved_ho_node_mode() : Settings::SLAVE_ALL;
1720 [ # # # ]: 0 : switch( ho_mode )
1721 : : {
1722 : : case Settings::SLAVE_ALL:
1723 [ # # ]: 0 : byteArray.resize( num_vertex, 0 );
1724 [ # # ][ # # ]: 0 : enslave_higher_order_nodes( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1725 : 0 : break;
1726 : : case Settings::SLAVE_NONE:
1727 : : // Do nothing. We clear other bits when processing the 'fixed' array above.
1728 : 0 : break;
1729 : : default:
1730 : : MSQ_SETERR( err )
1731 : : ( "Specified higher-order noded slaving scheme not supported "
1732 : : "when initializind PatchData using PatchData::fill",
1733 [ # # ][ # # ]: 0 : MsqError::NOT_IMPLEMENTED );
1734 : 0 : return;
1735 : : }
1736 : :
1737 [ # # ][ # # ]: 0 : this->initialize_data( arrptr( offsets ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1738 : :
1739 : : // initialize_data will re-order vertex handles and
1740 : : // update element connectivity accordingly. Use
1741 : : // the values we stored in vertexHandlesArray to
1742 : : // figure out the new index of each vertex, and initialize
1743 : : // the vertex.
1744 [ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
1745 [ # # ][ # # ]: 0 : vertexArray[i] = coords + 3 * (size_t)vertexHandlesArray[i];
[ # # ][ # # ]
1746 : :
1747 [ # # ][ # # ]: 0 : for( i = 0; i < num_vertex; ++i )
1748 [ # # ][ # # ]: 0 : vertexArray[i].flags() = byteArray[i];
[ # # ]
1749 : : }
1750 : :
1751 : 1128315 : void PatchData::make_handles_unique( Mesh::EntityHandle* handles, size_t& count, size_t* index_map )
1752 : : {
1753 [ - + ]: 2256630 : if( count < 2 ) { return; }
1754 : : // save this now, as we'll be changing count later
1755 : 1128315 : const size_t* index_end = index_map + count;
1756 : :
1757 [ + - ]: 1128315 : if( index_map )
1758 : : {
1759 : : // Copy input handle list into index map as a temporary
1760 : : // copy of the input list.
1761 : : assert( sizeof( Mesh::EntityHandle ) == sizeof( size_t ) );
1762 : 1128315 : memcpy( index_map, handles, count * sizeof( size_t ) );
1763 : : }
1764 : :
1765 : : // Make handles a unique list
1766 : 1128315 : std::sort( handles, handles + count );
1767 : 1128315 : Mesh::EntityHandle* end = std::unique( handles, handles + count );
1768 : 1128315 : count = end - handles;
1769 : :
1770 [ + - ]: 1128315 : if( index_map )
1771 : : {
1772 : : // Replace each handle in index_map with the index of
1773 : : // its position in the handles array
1774 : : Mesh::EntityHandle* pos;
1775 [ + + ]: 16993745 : for( size_t* iter = index_map; iter != index_end; ++iter )
1776 : : {
1777 [ + - ]: 15865430 : pos = std::lower_bound( handles, handles + count, (Mesh::EntityHandle)*iter );
1778 : 15865430 : *iter = pos - handles;
1779 : : }
1780 : : }
1781 : : }
1782 : :
1783 : 73 : void PatchData::fill_global_patch( MsqError& err )
1784 : : {
1785 [ + - ]: 73 : GlobalPatch gp;
1786 [ + - ][ + - ]: 73 : gp.set_mesh( get_mesh() );
1787 [ + - ][ + - ]: 146 : PatchIterator iter( &gp );
1788 [ + - ][ + - ]: 146 : bool b = iter.get_next_patch( *this, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
1789 [ - + ][ # # ]: 73 : if( !b ) MSQ_SETERR( err )( "Empty patch in PatchData::fill_global_patch", MsqError::INVALID_MESH );
[ # # ]
1790 [ - + ][ + - ]: 146 : assert( b );
1791 : : }
1792 : :
1793 : 1128315 : void PatchData::set_mesh_entities( std::vector< Mesh::ElementHandle >& elements,
1794 : : std::vector< Mesh::VertexHandle >& free_vertices, MsqError& err )
1795 : : {
1796 [ + - ]: 1128315 : Mesh* current_mesh = get_mesh();
1797 [ - + ]: 1128315 : if( !current_mesh )
1798 : : {
1799 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No Mesh associated with PatchData.", MsqError::INVALID_STATE );
1800 : 1128315 : return;
1801 : : }
1802 : :
1803 [ - + ]: 1128315 : if( elements.empty() )
1804 : : {
1805 [ # # ]: 0 : clear();
1806 : 0 : return;
1807 : : }
1808 : :
1809 [ + - ]: 1128315 : elementHandlesArray = elements;
1810 [ + - + - ]: 2256630 : get_mesh()->elements_get_attached_vertices( arrptr( elementHandlesArray ), elementHandlesArray.size(),
1811 [ + - ][ + - ]: 2256630 : vertexHandlesArray, offsetArray, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
1812 : :
1813 : : // Construct CSR-rep element connectivity data
1814 : 1128315 : size_t num_vert = vertexHandlesArray.size();
1815 [ + - ]: 1128315 : elemConnectivityArray.resize( num_vert );
1816 [ + - ][ + - ]: 1128315 : make_handles_unique( arrptr( vertexHandlesArray ), num_vert, arrptr( elemConnectivityArray ) );
[ + - ]
1817 [ + - ]: 1128315 : vertexHandlesArray.resize( num_vert );
1818 : :
1819 : : // Get element topologies
1820 [ + - ]: 1128315 : std::vector< EntityTopology > elem_topologies( elementHandlesArray.size() );
1821 [ + - ][ + - ]: 2256630 : get_mesh()->elements_get_topologies( arrptr( elementHandlesArray ), arrptr( elem_topologies ),
1822 [ + - ][ + - ]: 2256630 : elementHandlesArray.size(), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
1823 : :
1824 : : // get vertex bits from mesh
1825 [ + - ]: 1128315 : byteArray.resize( vertexHandlesArray.size() );
1826 [ + - ][ + - ]: 1128315 : get_mesh()->vertices_get_byte( arrptr( vertexHandlesArray ), arrptr( byteArray ), vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
1827 : :
1828 : : // if free_vertices is not empty, then we need to mark as
1829 : : // fixed any vertices *not* in that list.
1830 [ + + ]: 1128315 : if( free_vertices.size() == 1 )
1831 : : {
1832 [ + + ]: 8694834 : for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
1833 [ + - ][ + - ]: 7827359 : if( vertexHandlesArray[i] == free_vertices.front() )
[ + + ]
1834 [ + - ]: 867475 : byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
1835 : : else
1836 [ + - ]: 6959884 : byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
1837 : : }
1838 [ + + ]: 260840 : else if( !free_vertices.empty() )
1839 : : {
1840 : : // sort and remove duplicates from free_vertices list.
1841 [ + - ]: 2 : std::sort( free_vertices.begin(), free_vertices.end() );
1842 [ + - ][ + - ]: 2 : free_vertices.erase( std::unique( free_vertices.begin(), free_vertices.end() ), free_vertices.end() );
1843 : :
1844 [ + + ]: 1429 : for( size_t i = 0; i < vertexHandlesArray.size(); ++i )
1845 : : {
1846 [ + - ][ + - ]: 2854 : if( ( byteArray[i] & MsqVertex::MSQ_DEPENDENT ) ||
[ + - ][ + - ]
1847 [ + - ][ + - ]: 1427 : std::binary_search( free_vertices.begin(), free_vertices.end(), vertexHandlesArray[i] ) )
1848 [ + - ]: 1427 : byteArray[i] &= ~MsqVertex::MSQ_PATCH_FIXED;
1849 : : else
1850 [ # # ]: 0 : byteArray[i] |= MsqVertex::MSQ_PATCH_FIXED;
1851 : : }
1852 : : }
1853 : :
1854 : : // set element types
1855 [ + - ]: 1128315 : elementArray.resize( elementHandlesArray.size() );
1856 [ + + ]: 5069195 : for( size_t i = 0; i < elementHandlesArray.size(); ++i )
1857 [ + - ][ + - ]: 3940880 : elementArray[i].set_element_type( elem_topologies[i] );
[ + - ]
1858 : :
1859 : : // get element connectivity, group vertices by free/slave/fixed state
1860 [ + - ][ + - ]: 1128315 : initialize_data( arrptr( offsetArray ), arrptr( byteArray ), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
1861 : :
1862 : : // get vertex coordinates
1863 [ + - ]: 1128315 : vertexArray.resize( vertexHandlesArray.size() );
1864 [ + - ][ + - ]: 2256630 : get_mesh()->vertices_get_coordinates( arrptr( vertexHandlesArray ), arrptr( vertexArray ),
1865 [ + - ][ + - ]: 2256630 : vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
1866 : :
1867 : : // set vertex flags
1868 [ + + ][ + - ]: 10156385 : for( size_t i = 0; i < vertexArray.size(); ++i )
1869 [ + - ][ + - ]: 10156385 : vertexArray[i].flags() = byteArray[i];
[ + - ]
1870 : : }
1871 : :
1872 : 0 : void PatchData::get_sample_location( size_t element_index, Sample sample, Vector3D& result, MsqError& err ) const
1873 : : {
1874 [ # # ]: 0 : const MsqMeshEntity& elem = element_by_index( element_index );
1875 [ # # ]: 0 : const NodeSet ho_bits = non_slave_node_set( element_index );
1876 [ # # ][ # # ]: 0 : const MappingFunction* const f = get_mapping_function( elem.get_element_type() );
1877 [ # # ]: 0 : if( !f )
1878 : : {
1879 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No mapping function", MsqError::UNSUPPORTED_ELEMENT );
1880 : 0 : return;
1881 : : }
1882 : :
1883 : : double coeff[27];
1884 : : size_t num_coeff, index[27];
1885 [ # # ][ # # ]: 0 : f->coefficients( sample, ho_bits, coeff, index, num_coeff, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
1886 [ # # ][ # # ]: 0 : f->convert_connectivity_indices( elem.node_count(), index, num_coeff, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1887 : :
1888 [ # # ]: 0 : const size_t* const conn = elem.get_vertex_index_array();
1889 [ # # ]: 0 : assert( num_coeff > 0 );
1890 [ # # ][ # # ]: 0 : result = coeff[0] * vertex_by_index( conn[index[0]] );
[ # # ]
1891 [ # # ]: 0 : for( unsigned i = 1; i < num_coeff; ++i )
1892 [ # # ][ # # ]: 0 : result += coeff[i] * vertex_by_index( conn[index[i]] );
[ # # ]
1893 : : }
1894 : :
1895 : 26228601 : NodeSet PatchData::non_slave_node_set( size_t element_index ) const
1896 : : {
1897 [ + - ]: 26228601 : const MsqMeshEntity& elem = element_by_index( element_index );
1898 [ + - ]: 26228601 : EntityTopology type = elem.get_element_type();
1899 [ + - ]: 26228601 : const size_t* conn = elem.get_vertex_index_array();
1900 [ + - ]: 26228601 : const size_t n = elem.node_count();
1901 : :
1902 [ + - ]: 26228601 : MsqError err;
1903 : : bool have_midedge, have_midface, have_midelem;
1904 [ + - ]: 26228601 : unsigned num_edge = 0, num_face = 0, num_corner = TopologyInfo::corners( type );
1905 [ + - ]: 26228601 : TopologyInfo::higher_order( type, n, have_midedge, have_midface, have_midelem, err );
1906 [ + - ]: 26228601 : num_edge = TopologyInfo::edges( type );
1907 [ + - ][ + + ]: 26228601 : if( TopologyInfo::dimension( type ) == 2 )
1908 : 16235706 : num_face = 1;
1909 : : else
1910 [ + - ]: 9992895 : num_face = TopologyInfo::faces( type );
1911 : :
1912 [ + - ]: 26228601 : NodeSet result;
1913 [ + - ]: 26228601 : result.set_all_corner_nodes( type );
1914 [ + + ]: 26228601 : if( have_midedge )
1915 : : {
1916 [ + + ]: 407212 : for( unsigned i = 0; i < num_edge; ++i )
1917 : : {
1918 [ + - ][ + - ]: 310132 : if( !( vertex_by_index( conn[num_corner + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
[ + + ]
1919 [ + - ]: 251694 : result.set_mid_edge_node( i );
1920 : : }
1921 : : }
1922 [ - + ]: 26228601 : if( have_midface )
1923 : : {
1924 [ # # ]: 0 : for( unsigned i = 0; i < num_face; ++i )
1925 : : {
1926 [ # # ][ # # ]: 0 : if( !( vertex_by_index( conn[num_corner + num_edge + i] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
[ # # ]
1927 [ # # ]: 0 : result.set_mid_face_node( i );
1928 : : }
1929 : : }
1930 [ - + ][ # # ]: 26228601 : if( have_midelem &&
[ - + ]
1931 [ # # ][ # # ]: 0 : !( vertex_by_index( conn[num_corner + num_edge + num_face] ).get_flags() & MsqVertex::MSQ_DEPENDENT ) )
1932 [ # # ]: 0 : result.set_mid_region_node();
1933 : :
1934 : 26228601 : return result;
1935 : : }
1936 : :
1937 : 450 : void PatchData::get_samples( size_t element, std::vector< Sample >& samples, MsqError& ) const
1938 : : {
1939 [ + - ]: 450 : NodeSet ns = get_samples( element );
1940 [ + - ][ + - ]: 450 : samples.resize( ns.num_nodes() );
1941 : 450 : std::vector< Sample >::iterator i = samples.begin();
1942 : :
1943 : : unsigned j;
1944 [ + - ][ + - ]: 450 : EntityTopology type = element_by_index( element ).get_element_type();
1945 [ + - ][ + + ]: 1850 : for( j = 0; j < TopologyInfo::corners( type ); ++j )
1946 [ + - ][ + + ]: 1400 : if( ns.corner_node( j ) ) *( i++ ) = Sample( 0, j );
[ + - ][ + - ]
[ + - ]
1947 [ + - ][ + + ]: 1850 : for( j = 0; j < TopologyInfo::edges( type ); ++j )
1948 [ + - ][ - + ]: 1400 : if( ns.mid_edge_node( j ) ) *( i++ ) = Sample( 1, j );
[ # # ][ # # ]
[ # # ]
1949 [ + - ][ - + ]: 450 : if( TopologyInfo::dimension( type ) == 3 )
1950 : : {
1951 [ # # ][ # # ]: 0 : for( j = 0; j < TopologyInfo::faces( type ); ++j )
1952 [ # # ][ # # ]: 0 : if( ns.mid_face_node( j ) ) *( i++ ) = Sample( 2, j );
[ # # ][ # # ]
[ # # ]
1953 [ # # ][ # # ]: 0 : if( ns.mid_region_node() ) *( i++ ) = Sample( 3, 0 );
[ # # ][ # # ]
[ # # ]
1954 : : }
1955 [ + - ][ + + ]: 450 : else if( ns.mid_face_node( 0 ) )
1956 [ + - ][ + - ]: 400 : *( i++ ) = Sample( 2, 0 );
[ + - ]
1957 : :
1958 [ + - ][ - + ]: 450 : assert( i == samples.end() );
1959 : 450 : }
1960 : :
1961 : 2 : bool PatchData::attach_extra_data( ExtraData* data )
1962 : : {
1963 [ - + ]: 2 : if( data->patchNext ) { return false; }
1964 : :
1965 [ - + ]: 2 : if( !data->patchPtr )
1966 : 0 : data->patchPtr = this;
1967 [ - + ]: 2 : else if( data->patchPtr != this )
1968 : 0 : return false;
1969 : :
1970 : 2 : data->patchNext = dataList;
1971 : 2 : dataList = data;
1972 : 2 : return true;
1973 : : }
1974 : :
1975 : 0 : bool PatchData::remove_extra_data( ExtraData* data )
1976 : : {
1977 [ # # ]: 0 : if( data->patchPtr != this ) return false;
1978 : :
1979 [ # # ]: 0 : if( dataList == data )
1980 : : {
1981 : 0 : dataList = data->patchNext;
1982 : 0 : data->patchNext = 0;
1983 : 0 : data->patchPtr = 0;
1984 : 0 : return true;
1985 : : }
1986 : :
1987 [ # # ]: 0 : for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
1988 [ # # ]: 0 : if( iter->patchNext == data )
1989 : : {
1990 : 0 : iter->patchNext = data->patchNext;
1991 : 0 : data->patchNext = 0;
1992 : 0 : data->patchPtr = 0;
1993 : 0 : return true;
1994 : : }
1995 : :
1996 : 0 : return false;
1997 : : }
1998 : :
1999 : 1128371 : void PatchData::notify_new_patch()
2000 : : {
2001 [ - + ]: 1128371 : for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
2002 : 0 : iter->notify_new_patch();
2003 : 1128371 : }
2004 : :
2005 : 0 : void PatchData::notify_sub_patch( PatchData& sub_patch, const size_t* vertex_map, const size_t* element_map,
2006 : : MsqError& err )
2007 : : {
2008 [ # # ]: 0 : for( ExtraData* iter = dataList; iter; iter = iter->patchNext )
2009 : : {
2010 [ # # ][ # # ]: 0 : iter->notify_sub_patch( sub_patch, vertex_map, element_map, err );MSQ_ERRRTN( err );
[ # # ]
2011 : : }
2012 : : }
2013 : :
2014 : 1871 : void PatchData::notify_patch_destroyed()
2015 : : {
2016 : : // Remove all ExtraDatas from list and notify them
2017 : : // that they are being destroyed.
2018 [ + + ]: 1873 : while( dataList )
2019 : : {
2020 : 2 : ExtraData* dead_data = dataList;
2021 : 2 : dataList = dataList->patchNext;
2022 : 2 : dead_data->patchNext = 0;
2023 : 2 : dead_data->patchPtr = 0;
2024 : 2 : dead_data->notify_patch_destroyed();
2025 : : }
2026 : 1871 : }
2027 : :
2028 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
