Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
smooth_medial.C
Go to the documentation of this file.
1 /* *******************************************************************
2  * Rocstar Simulation Suite *
3  * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4  * *
5  * Illinois Rocstar LLC *
6  * Champaign, IL *
7  * www.illinoisrocstar.com *
8  * sales@illinoisrocstar.com *
9  * *
10  * License: See LICENSE file in top level of distribution package or *
11  * http://opensource.org/licenses/NCSA *
12  *********************************************************************/
13 /* *******************************************************************
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17  * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18  * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21  * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22  *********************************************************************/
23 #include "Rocmop.h"
24 #include "roccom.h"
25 #include "Pane.h"
26 #include "Rocblas.h"
27 #include "Rocmap.h"
28 #include "geometry.h"
29 #include "Rocsurf.h"
30 #include "Generic_element_2.h"
31 
32 MOP_BEGIN_NAMESPACE
33 
34 using MAP::Rocmap;
35 
36 //============== Debug features =======
37 // Export debugging information to user window, if user has registered
38 // "back-door" attributes to get debug info such as normals.
39 // Should be set to true for development purpose.
40 #define EXPORT_DEBUG_INFO 1
41 
42 static const double half_pi = 1.5707963267949;
43 static const double pi = 3.14159265358979;
44 
45 inline double SQR( double x) { return (x*x); }
46 
47 // ----------------------------------------------------------------------------
48 void dsytrd3(double A[3][3], double Q[3][3], double d[3], double e[2])
49 // ----------------------------------------------------------------------------
50 // Reduces a symmetric 3x3 matrix to tridiagonal form by applying
51 // (unitary) Householder transformations:
52 // [ d[0] e[0] ]
53 // A = Q . [ e[0] d[1] e[1] ] . Q^T
54 // [ e[1] d[2] ]
55 // The function accesses only the diagonal and upper triangular parts of
56 // A. The access is read-only.
57 // ----------------------------------------------------------------------------
58 // Note: Based on implementat of Joachim Kopp (www.mpi-hd.mpg.de/~jkopp/3x3/)
59 // ---------------------------------------------------------------------------
60 {
61  double u[3], q[3];
62  double omega, f;
63  double K, h, g;
64 
65  // Initialize Q to the identitity matrix
66  for (int i=0; i < 3; i++)
67  {
68  Q[i][i] = 1.0;
69  for (int j=0; j < i; j++)
70  Q[i][j] = Q[j][i] = 0.0;
71  }
72 
73  // Bring first row and column to the desired form
74  h = SQR(A[0][1]) + SQR(A[0][2]);
75  if (A[0][1] > 0)
76  g = -sqrt(h);
77  else
78  g = sqrt(h);
79  e[0] = g;
80  f = g * A[0][1];
81  u[1] = A[0][1] - g;
82  u[2] = A[0][2];
83 
84  omega = h - f;
85  if (omega > 0.0)
86  {
87  omega = 1.0 / omega;
88  K = 0.0;
89  for (int i=1; i < 3; i++)
90  {
91  f = A[1][i] * u[1] + A[i][2] * u[2];
92  q[i] = omega * f; // p
93  K += u[i] * f; // u* A u
94  }
95  K *= 0.5 * SQR(omega);
96 
97  for (int i=1; i < 3; i++)
98  q[i] = q[i] - K * u[i];
99 
100  d[0] = A[0][0];
101  d[1] = A[1][1] - 2.0*q[1]*u[1];
102  d[2] = A[2][2] - 2.0*q[2]*u[2];
103 
104  // Store inverse Householder transformation in Q
105  for (int j=1; j < 3; j++)
106  {
107  f = omega * u[j];
108  for (int i=1; i < 3; i++)
109  Q[i][j] = Q[i][j] - f*u[i];
110  }
111 
112  // Calculate updated A[1][2] and store it in e[1]
113  e[1] = A[1][2] - q[1]*u[2] - u[1]*q[2];
114  }
115  else
116  {
117  for (int i=0; i < 3; i++)
118  d[i] = A[i][i];
119  e[1] = A[1][2];
120  }
121 }
122 
123 // ----------------------------------------------------------------------------
124 int dsyevq3(double A[3][3], double Q[3][3], double w[3])
125 // ----------------------------------------------------------------------------
126 // Calculates the eigenvalues and normalized eigenvectors of a symmetric 3x3
127 // matrix A using the QL algorithm with implicit shifts, preceded by a
128 // Householder reduction to tridiagonal form.
129 // The function accesses only the diagonal and upper triangular parts of A.
130 // The access is read-only.
131 // ----------------------------------------------------------------------------
132 // Parameters:
133 // A: The symmetric input matrix
134 // Q: Storage buffer for eigenvectors
135 // w: Storage buffer for eigenvalues
136 // ----------------------------------------------------------------------------
137 // Return value:
138 // 0: Success
139 // -1: Error (no convergence)
140 // ----------------------------------------------------------------------------
141 // Dependencies:
142 // dsytrd3()
143 // ----------------------------------------------------------------------------
144 // Note: Based on implementat of Joachim Kopp (www.mpi-hd.mpg.de/~jkopp/3x3/)
145 // ----------------------------------------------------------------------------
146 {
147  double e[3]; // The third element is used only as temporary workspace
148  double g, r, p, f, b, s, c, t; // Intermediate storage
149  int nIter;
150  int m;
151 
152  // Transform A to real tridiagonal form by the Householder method
153  dsytrd3(A, Q, w, e);
154 
155  // Calculate eigensystem of the remaining real symmetric tridiagonal matrix
156  // with the QL method
157  //
158  // Loop over all off-diagonal elements
159  for (int l=0; l < 2; l++)
160  {
161  nIter = 0;
162  for (;;)
163  {
164  // Check for convergence and exit iteration loop if off-diagonal
165  // element e(l) is zero
166  for (m=l; m < 2; m++)
167  {
168  g = fabs(w[m])+fabs(w[m+1]);
169  if (fabs(e[m]) + g == g)
170  break;
171  }
172  if (m == l)
173  break;
174 
175  if (nIter++ >= 30)
176  return -1;
177 
178  // Calculate g = d_m - k
179  g = (w[l+1] - w[l]) / (e[l] + e[l]);
180  r = sqrt(SQR(g) + 1.0);
181  if (g > 0)
182  g = w[m] - w[l] + e[l]/(g + r);
183  else
184  g = w[m] - w[l] + e[l]/(g - r);
185 
186  s = c = 1.0;
187  p = 0.0;
188  for (int i=m-1; i >= l; i--)
189  {
190  f = s * e[i];
191  b = c * e[i];
192  if (fabs(f) > fabs(g))
193  {
194  c = g / f;
195  r = sqrt(SQR(c) + 1.0);
196  e[i+1] = f * r;
197  c *= (s = 1.0/r);
198  }
199  else
200  {
201  s = f / g;
202  r = sqrt(SQR(s) + 1.0);
203  e[i+1] = g * r;
204  s *= (c = 1.0/r);
205  }
206 
207  g = w[i+1] - p;
208  r = (w[i] - g)*s + 2.0*c*b;
209  p = s * r;
210  w[i+1] = g + p;
211  g = c*r - b;
212 
213  // Form eigenvectors
214  for (int k=0; k < 3; k++)
215  {
216  t = Q[k][i+1];
217  Q[k][i+1] = s*Q[k][i] + c*t;
218  Q[k][i] = c*Q[k][i] - s*t;
219  }
220  }
221  w[l] -= p;
222  e[l] = g;
223  e[m] = 0.0;
224  }
225  }
226 
227  return 0;
228 }
229 
230 void Rocmop::smooth_surf_medial(){
231  if(_verb > 1)
232  std::cout << "Entering Rocmop::smooth_medial" << std::endl;
233 
234  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");
235 
236  // Calculate the face normals
237  Rocmop::evaluate_face_normals();
238 
239  // Compute the medial quadric
240  Rocmop::compute_medial_quadric();
241 
242 #if EXPORT_DEBUG_INFO
243  COM::Attribute *vnormals = const_cast<COM::Attribute*>(_usr_window->attribute("vnormals"));
244  COM::Attribute *vnrms = _wrk_window->attribute("vnormals");
245  COM::Attribute *awnormals = const_cast<COM::Attribute*>(_usr_window->attribute("angleweighted"));
246  COM::Attribute *awnrms = _wrk_window->attribute("awnormals");
247  COM::Attribute *uwnormals = const_cast<COM::Attribute*>(_usr_window->attribute("unitweighted"));
248  COM::Attribute *uwnrms = _wrk_window->attribute("uwnormals");
249  COM::Attribute *tangranks = const_cast<COM::Attribute*>(_usr_window->attribute("tangranks"));
250  COM::Attribute *tranks = _wrk_window->attribute("tangranks");
251  COM::Attribute *evals = _wrk_window->attribute("lambda");
252  COM::Attribute *eigenvalues = const_cast<COM::Attribute*>(_usr_window->attribute("eigenvalues"));
253  COM::Attribute *ws = _wrk_window->attribute("weights");
254  COM::Attribute *weights = const_cast<COM::Attribute*>(_usr_window->attribute("adefects"));
255 
256  if ( vnormals) Rocblas::copy( vnrms, vnormals);
257  if ( awnormals) Rocblas::copy( awnrms, awnormals);
258  if ( uwnormals) Rocblas::copy( uwnrms, uwnormals);
259  if ( tangranks) Rocblas::copy( tranks, tangranks);
260  if ( eigenvalues) Rocblas::copy(evals, eigenvalues);
261  if ( weights) Rocblas::copy(ws,weights);
262 #endif
263 
264  // Identify ridge edges
265  //identify_ridge_edges();
266 
267  // Redistribute vertices within their tangent spaces
268  // FIXME.... MODIFY THESE TO USE PN TRIANGLES
269  for ( int i=0; i<_rediter; ++i) {
270  redistribute_vertices_ridge();
271  redistribute_vertices_smooth();
272  }
273 
274  // Add the displacements onto the nodal coords
275  // FIXME
276  const std::string att1("disps");
277  COM::Attribute *disps = _wrk_window->attribute(att1);
278  COM::Attribute *nc = _wrk_window->attribute( COM::COM_NC);
279  // Rocblas::add(disps,nc,nc);
280 
281  if(_verb > 2)
282  std::cout << "Exiting Rocmop::smooth_medial" << std::endl;
283 }
284 
285 void Rocmop::
286 compute_medial_quadric() {
287  if(_verb > 1)
288  std::cout << "Entering Rocmop::compute_medial_quadric" << std::endl;
289  // Use the following buffer spaces: _eigvecs for matrix A;
290  // _eigvalues to store bm; disps to store bo, _weights to store r
291  const std::string att1("eigvecs"), att2("lambda"), att3("weights"),
292  att4("facenormals"), att5("facecenters"), att6("tangranks"),
293  att7("cntnranks"), att8("cntnvecs"), att9("vnormals"), att10("disps"),
294  att11("awnormals"), att12("uwnormals");
295  COM::Attribute *A_attr = _wrk_window->attribute(att1);
296  COM::Attribute *bm_attr = _wrk_window->attribute(att2);
297  COM::Attribute *aw_attr = _wrk_window->attribute(att11);
298  COM::Attribute *uw_attr = _wrk_window->attribute(att12);
299  COM::Attribute *r_attr = _wrk_window->attribute(att3);
300 
301  int facenormals_id = _wrk_window->attribute(att4)->id();
302  int facecenters_id = _wrk_window->attribute(att5)->id();
303  int tangranks_id = _wrk_window->attribute(att6)->id();
304  int cntnranks_id = _wrk_window->attribute(att7)->id();
305  int cntnvecs_id = _wrk_window->attribute(att8)->id();
306  int vnormals_id = _wrk_window->attribute(att9)->id();
307  int awnormals_id = _wrk_window->attribute(att11)->id();
308  int uwnormals_id = _wrk_window->attribute(att12)->id();
309  int disps_id = _wrk_window->attribute(att10)->id();
310 
311  double zero = 0.;
312  Rocblas::copy_scalar( &zero, A_attr);
313  Rocblas::copy_scalar( &zero, bm_attr);
314  Rocblas::copy_scalar( &zero, r_attr);
315 
316  // Initialize disps to 0
317  std::vector< COM::Pane*> allpanes;
318  _wrk_window->panes(allpanes);
319  for(int i =0, local_npanes = allpanes.size();
320  i<local_npanes; ++i){
321  Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
322  ( allpanes[i]->attribute(disps_id)->pointer());
323  for(int i =0, ni = allpanes[i]->size_of_real_nodes(); i<ni; ++i){
324  ds[i] = Vector_3<double>(0.0,0.0,0.0);
325  }
326  }
327 
328  // Loop through the panes and its real faces, calculate A, b, and c = sum_i(wi)
329  std::vector< COM::Pane*>::iterator it = allpanes.begin();
330  SURF::Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();
331  for ( int i=0, local_npanes = allpanes.size();
332  i<local_npanes; ++i, ++it, ++pm_it) {
333  COM::Pane *pane = *it;
334 
335  // Obtain nodal coordinates of current pane, assuming contiguous layout
336  const Vector_3<double> *pnts = reinterpret_cast<Vector_3<double>*>
337  (pane->attribute(COM_NC)->pointer());
338  Vector_3<double> *nrms = reinterpret_cast<Vector_3<double>*>
339  ( pane->attribute(facenormals_id)->pointer());
340 
341  // Obtain pointers to A_attr, bm_attr, and r_attr
342  Vector_3<double> *As = reinterpret_cast<Vector_3<double>*>
343  ( pane->attribute(A_attr->id())->pointer());
344  Vector_3<double> *bs_m = reinterpret_cast<Vector_3<double>*>
345  ( pane->attribute(bm_attr->id())->pointer());
346  Vector_3<double> *aws = reinterpret_cast<Vector_3<double>*>
347  ( pane->attribute(aw_attr->id())->pointer());
348  Vector_3<double> *uws = reinterpret_cast<Vector_3<double>*>
349  ( pane->attribute(uw_attr->id())->pointer());
350  double *rs = reinterpret_cast<double*>
351  ( pane->attribute(r_attr->id())->pointer());
352  Vector_3<double> *cnts = reinterpret_cast<Vector_3<double>*>
353  ( pane->attribute(facecenters_id)->pointer());
354 
355  // Loop through real elements of the current pane
356  Element_node_enumerator ene( pane, 1);
357  for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
358  // If speed is uniform within a face, then unit normal is the same
359  // after propagation.
360  const Vector_3<double> &ns_nz = nrms[j];
361  Vector_3<double> nrm = ns_nz;
362 
363  int ne = ene.size_of_edges();
364  int uindex=ene[ne-1]-1, vindex=ene[0]-1;
365 
366  // Loop through all vertices of current face.
367  for ( int k=0; k<ne; ++k) {
368  int windex = ene[(k+1==ne)?0:k+1]-1;
369 
370  Vector_3<double> cnt = cnts[j];
371 
372  // Calculate the weight of the current face for this node.
373  double w=1;
374 
375 
376  Vector_3<double> ts[2] = { pnts[uindex]-pnts[vindex],
377  pnts[windex]-pnts[vindex]};
378 
379  double s = std::sqrt((ts[0]*ts[0])*(ts[1]*ts[1]));
380  double angle=0;
381  if ( s>=0) {
382  double cosw = ts[0]*ts[1]/s;
383  if (cosw>1) cosw=1; else if ( cosw<-1) cosw=-1;
384  angle = std::acos( cosw);
385  }
386 
387  if ( _wght_scheme!=SURF::E2N_ONE) {
388 
389  switch ( _wght_scheme) {
390  case SURF::E2N_ANGLE:
391  case SURF::E2N_SPHERE: {
392 
393  if ( _wght_scheme==SURF::E2N_SPHERE)
394  w = std::sin(angle)/s;
395  else if ( w>half_pi)
396  w = std::max(0.,pi-angle);
397  else
398  w = angle;
399  break;
400  }
401  case SURF::E2N_AREA: {
402  w = Vector_3<double>::cross_product(ts[0],ts[1]).norm()/2.; break;
403  }
404  }
405  nrm = w*ns_nz;
406  }
407 
408  // Obtain the addresses of the linear system and weight
409  Vector_3<double> *A_v = &As[3*vindex];
410  Vector_3<double> &b_mv = bs_m[vindex];
411  Vector_3<double> &aw = aws[vindex];
412  Vector_3<double> &uw = uws[vindex];
413  double &r_v = rs[vindex];
414 
415  // Update A, b_m, and r for the current vertex
416  A_v[0] += ns_nz[0]*nrm;
417  A_v[1] += ns_nz[1]*nrm;
418  A_v[2] += ns_nz[2]*nrm;
419  b_mv -= nrm;
420  aw += (std::sin(angle)/s)*ns_nz;
421  uw += ns_nz;
422  r_v += angle;
423 
424  // Update indices for next node
425  uindex=vindex; vindex=windex;
426  }
427 
428 #if defined(BOUNDARY_TREATMENT) && BOUNDARY_TREATMENT
429  // Loop through the halfedges of the face to handle border edges
430  Halfedge h( &*pm_it, Edge_ID(j+1,0), SURF::REAL_PANE), h0=h;
431  do {
432  if ( h.is_border()) {
433  // Incorporate normal plane passing through the border edge into
434  // its incident vertices with weight pi/2.
435  Vector_3<double> nrm_nz = Vector_3<double>::cross_product( ns_nz, h.tangent()).normalize();
436  Vector_3<double> nrm = half_pi*nrm_nz;
437 
438  // Update the origin and destination of h with the orthogonal pane
439  for ( int k=0; k<2; ++k) {
440  int vindex = ene[ (k?h:h.next()).id().lid()]-1;
441  // Obtain the addresses of the linear system and weight
442  Vector_3<double> *A_v = &As[3*vindex];
443  Vector_3<double> &b_mv = bs_m[vindex];
444  double &r_v = rs[vindex];
445 
446  // Update A, b_m, and r for the current vertex.
447  // No need to update as sigma is 0 for it.
448  A_v[0] += nrm_nz[0]*nrm;
449  A_v[1] += nrm_nz[1]*nrm;
450  A_v[2] += nrm_nz[2]*nrm;
451  b_mv -= nrm;
452  r_v += half_pi;
453  }
454  }
455  } while ( (h=h.next()) != h0);
456 #endif
457  }
458  }
459 
460  // Reduce on shared nodes for A, b_m, and r
461  _wm->reduce_on_shared_nodes( A_attr, SURF::Window_manifold_2::OP_SUM);
462  _wm->reduce_on_shared_nodes( bm_attr, SURF::Window_manifold_2::OP_SUM);
463  _wm->reduce_on_shared_nodes( aw_attr, SURF::Window_manifold_2::OP_SUM);
464  _wm->reduce_on_shared_nodes( uw_attr, SURF::Window_manifold_2::OP_SUM);
465  _wm->reduce_on_shared_nodes( r_attr, SURF::Window_manifold_2::OP_SUM);
466 
467  // Loop through the panes and its real vertices
468  it = allpanes.begin();
469  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
470 
471  COM::Pane *pane = *it;
472 
473  // Obtain pointers to A_attr, bo_attr, bm_attr, and r_attr
474  Vector_3<double> *As = reinterpret_cast<Vector_3<double>*>
475  ( pane->attribute(A_attr->id())->pointer());
476  Vector_3<double> *bs_m = reinterpret_cast<Vector_3<double>*>
477  ( pane->attribute(bm_attr->id())->pointer());
478  double *rs = reinterpret_cast<double*>
479  ( pane->attribute(r_attr->id())->pointer());
480  int *tranks = reinterpret_cast<int*>
481  ( pane->attribute(tangranks_id)->pointer());
482  int *cranks = reinterpret_cast<int*>
483  ( pane->attribute(cntnranks_id)->pointer());
484  Vector_3<double> *cvecs = reinterpret_cast<Vector_3<double>*>
485  ( pane->attribute(cntnvecs_id)->pointer());
486  int *cnstrs = _cnstr_types?reinterpret_cast<int*>
487  ( pane->attribute(_cnstr_types->id())->pointer()):NULL;
488  Vector_3<double> *cnstr_dirs = _cnstr_dirs?reinterpret_cast<Vector_3<double>*>
489  ( pane->attribute(_cnstr_dirs->id())->pointer()):NULL;
490  Vector_3<double> *vnrms = reinterpret_cast<Vector_3<double>*>
491  ( pane->attribute(vnormals_id)->pointer());
492  Vector_3<double> *awnrms = reinterpret_cast<Vector_3<double>*>
493  ( pane->attribute(awnormals_id)->pointer());
494  Vector_3<double> *uwnrms = reinterpret_cast<Vector_3<double>*>
495  ( pane->attribute(uwnormals_id)->pointer());
496 
497 
498  // Loop through all real nodes of the pane
499  //std::cout << "Number of nodes = " << pane->size_of_real_nodes() << std::endl;
500  for ( int j=0, jn=pane->size_of_real_nodes(); j<jn; ++j) {
501 
502  // Skip vertices with zero weight, such as those at
503  // edge- and face-centers for quadratic elements
504  if ( rs[j]==0) continue;
505  Vector_3<double> *es = &As[3*j]; // eigen-vectors of current vertex.
506  Vector_3<double> *cs = &cvecs[2*j]; // cnstr-tangent of current vertex.
507  // Make copy of matrix A as it will be overwritten by eigenvectors
508  Vector_3<double> A_v[3] = {es[0], es[1], es[2]};
509  // Likewise, b will be overwritten by eigenvalues
510  Vector_3<double> b_v = bs_m[j];
511  awnrms[j].normalize();
512  uwnrms[j].normalize();
513 
514 
515  // Perform eigen-analysis for the vertex using medial quadric
516  // orank is the rank of the orthogonal (primar) space of A.
517  // The non-constrained solution is stored in vnrm.
518  int orank = eigen_analyze_vertex( es, bs_m[j], &vnrms[j], 2*pi-rs[j]);
519  // EDIT
520  //orank = (orank == 3) ? 2 : orank;
521 
522  //COM_assertion_msg((orank <=3) && (orank >= 1), "Orank is invalid\n");
523 
524  // Store the rank of the tangent space.
525  tranks[j] = 3-orank;
526 
527  // Copy eigenvectors into cvecs.
528  for ( int k=tranks[j]-1; k>=0; --k)
529  cs[k] = es[orank+k];
530 
531  // Compute normal estimate using medial quadric,
532  // and save it into vnrms (overwrite the non-constrained solution)
533  Vector_3<double> d(0,0,0);
534  // Set the tolerance to avoid perturbation due to small eigenvalues
535  double tol = bs_m[j][0]*5.e-3;
536  for ( int k=0; k<orank; ++k)
537  d += (es[k]*b_v/-std::max(bs_m[j][k],tol))*es[k];
538  vnrms[j] = d;
539  }
540  }
541  if(_verb > 1)
542  std::cout << "Exiting Rocmop::compute_medial_quadric" << std::endl;
543 }
544 
545 // Regularize all nodes in their tangent spaces
546 void Rocmop::
547 redistribute_vertices_ridge() {
548  // Use the following buffer spaces: _eigvalues for vector c (only first one)
549 
550  const std::string att1("disps"), att2("lambda"), att3("weights"),
551  att4("cntnvecs"), att5("cntnranks");
552  int disps_id = _wrk_window->attribute(att1)->id();
553  // c_attr are the edge centers
554  // w_attr is the nodal safe factor
555  COM::Attribute *c_attr = _wrk_window->attribute(att2);
556  COM::Attribute *w_attr = _wrk_window->attribute(att3);
557  int cntnvecs_id = _wrk_window->attribute(att4)->id();
558  int cntnranks_id = _wrk_window->attribute(att5)->id();
559 
560  double zero = 0.;
561  Rocblas::copy_scalar( &zero, c_attr);
562 
563  // Loop through the panes and its real faces
564  std::vector< COM::Pane*> allpanes;
565  _wrk_window->panes(allpanes);
566  std::vector< COM::Pane*>::iterator it = allpanes.begin();
567  Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();
568  for ( int i=0, local_npanes = allpanes.size();
569  i<local_npanes; ++i, ++it, ++pm_it) {
570  COM::Pane *pane = *it;
571 
572  Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
573  ( pane->attribute(cntnvecs_id)->pointer());
574  const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
575  (pane->attribute(COM_NC)->pointer());
576  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
577  ( pane->attribute(c_attr->id())->pointer());
578  int *cranks = reinterpret_cast<int*>
579  ( pane->attribute(cntnranks_id)->pointer());
580 
581  std::set< Edge_ID> &eset_pn = _edges[i];
582  for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
583  eend=eset_pn.end(); eit!=eend; ++eit) {
584  Edge_ID eid = *eit;
585  bool is_border = eid.is_border();
586  if ( is_border) eid = pm_it->get_opposite_real_edge( eid);
587 
588  // Loop through real elements of the current pane
589  Element_node_enumerator ene( pane, eid.eid());
590  int ne = ene.size_of_edges();
591 
592  int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
593  if ( is_border) std::swap( vindex, windex);
594 
595  if ( cranks[vindex]>0) {
596  Vector_3<double> diff_wv = pnts[windex]-pnts[vindex];
597  // OLD:: Vector_3<double> diff_wv = pnts[windex]-pnts[vindex] + ds[windex] - ds[vindex];
598 
599  // get average of edge centers
600  cs[vindex][0] += 0.5*(es[2*vindex]*diff_wv);
601  }
602  }
603  }
604 
605  _wm->reduce_on_shared_nodes( c_attr, Window_manifold_2::OP_SUM);
606 
607  // Loop through all vertices to computer center
608  it = allpanes.begin(); pm_it=_wm->pm_begin();
609  for ( int i=0, local_npanes = allpanes.size();
610  i<local_npanes; ++i, ++it, ++pm_it) {
611  COM::Pane *pane = *it;
612 
613  Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
614  ( pane->attribute(cntnvecs_id)->pointer());
615  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
616  ( pane->attribute( c_attr->id())->pointer());
617  int *cranks = reinterpret_cast<int*>
618  ( pane->attribute(cntnranks_id)->pointer());
619 
620  // Loop through all real nodes of the pane
621  std::set< Edge_ID> &eset_pn = _edges[i];
622  for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
623  eend=eset_pn.end(); eit!=eend; ++eit) {
624  Edge_ID eid = *eit;
625  bool is_border = eid.is_border();
626  if ( is_border) eid = pm_it->get_opposite_real_edge( eid);
627 
628  // Loop through real elements of the current pane
629  Element_node_enumerator ene( pane, eid.eid());
630  int ne = ene.size_of_edges();
631 
632  int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
633  if ( is_border) std::swap( vindex, windex);
634 
635  if ( cranks[vindex]>0)
636  cs[vindex] = 0.5*cs[vindex][0]*es[2*vindex];
637  }
638  }
639 
640  // Loop through the panes and its real faces to compute safe factors.
641  get_redist_safe_factor( c_attr, w_attr, 1);
642 
643  // Loop through all vertices to update displacement vector
644  it = allpanes.begin(); pm_it=_wm->pm_begin();
645  for ( int i=0, local_npanes = allpanes.size();
646  i<local_npanes; ++i, ++it, ++pm_it) {
647  COM::Pane *pane = *it;
648 
649  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
650  ( pane->attribute( c_attr->id())->pointer());
651  Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
652  ( pane->attribute( disps_id)->pointer());
653  double *ws = reinterpret_cast<double*>
654  ( pane->attribute(w_attr->id())->pointer());
655 
656  // Loop through all real nodes of the pane
657  std::set< Edge_ID> &eset_pn = _edges[i];
658  for ( std::set< Edge_ID>::const_iterator eit=eset_pn.begin(),
659  eend=eset_pn.end(); eit!=eend; ++eit) {
660  Edge_ID eid = *eit;
661  bool is_border = eid.is_border();
662  if ( is_border) eid = pm_it->get_opposite_real_edge( eid);
663 
664  // Loop through real elements of the current pane
665  Element_node_enumerator ene( pane, eid.eid());
666  int ne = ene.size_of_edges();
667 
668  int vindex=ene[eid.lid()]-1, windex=ene[(eid.lid()+1)%ne]-1;
669  if ( is_border) std::swap( vindex, windex);
670 
671  double w=0.5;
672  if ( ws[vindex]>0) w = std::min( w, 1./ws[vindex]);
673  // displacement = w * average of two edge segments
674  ds[vindex] += w*cs[vindex];
675  // do PN triangle stuff here
676  }
677  }
678 }
679 
680 // Regularize all nodes in their tangent spaces
681 void Rocmop::
682 redistribute_vertices_smooth() {
683 
684  // Use the following buffer spaces: _eigvalues for vector c
685  // (only lsast two components); _weights to store weight
686  const std::string att1("disps"), att2("lambda"), att3("weights"),
687  att4("cntnvecs"), att5("cntnranks"), att6("tangranks");
688  int disps_id = _wrk_window->attribute(att1)->id();
689  COM::Attribute *c_attr = _wrk_window->attribute(att2);
690  COM::Attribute *w_attr = _wrk_window->attribute(att3);
691  int cntnvecs_id = _wrk_window->attribute(att4)->id();
692  int cntnranks_id = _wrk_window->attribute(att5)->id();
693  int tangranks_id = _wrk_window->attribute("tangranks")->id();
694 
695  double zero = 0.;
696  Rocblas::copy_scalar( &zero, c_attr);
697  Rocblas::copy_scalar( &zero, w_attr);
698 
699  std::vector< COM::Pane*> allpanes;
700  _wrk_window->panes(allpanes);
701  // Loop through the panes and its real faces
702  std::vector< COM::Pane*>::iterator it = allpanes.begin();
703  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
704  COM::Pane *pane = *it;
705 
706  Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
707  ( pane->attribute(cntnvecs_id)->pointer());
708  const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
709  (pane->attribute(COM_NC)->pointer());
710  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
711  ( pane->attribute(c_attr->id())->pointer());
712  double *ws = reinterpret_cast<double*>
713  ( pane->attribute(w_attr->id())->pointer());
714  Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
715  ( pane->attribute(disps_id)->pointer());
716  int *tranks = reinterpret_cast<int*>
717  ( pane->attribute(tangranks_id)->pointer());
718  int *cranks = reinterpret_cast<int*>
719  ( pane->attribute(cntnranks_id)->pointer());
720 
721  // Loop through real elements of the current pane
722  Element_node_enumerator ene( pane, 1);
723  Element_node_vectors_k_const<Point_3<double> > ps; ps.set( pnts, ene, 1);
724 
725  for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
726  int ne = ene.size_of_edges();
727  int uindex=ene[ne-1]-1, vindex=ene[0]-1;
728 
729  // Loop through all vertices of current face.
730  for ( int k=0; k<ne; ++k) {
731  int windex = ene[(k+1==ne)?0:k+1]-1;
732 
733  // Only redistribute smooth vertices
734  if ( tranks[vindex]==2) {
735  std::cout << "cranks[vindex]-1 = " << cranks[vindex]-1 << std::endl;
736  Vector_3<double> tngs[2] = {pnts[uindex]-pnts[vindex],
737  pnts[windex]-pnts[vindex]};
738  Vector_3<double> disp_v = ds[vindex];
739  Vector_3<double> diff_uv = ds[uindex]-disp_v+tngs[0];
740  Vector_3<double> diff_wv = ds[windex]-disp_v+tngs[1];
741 
742  // Compute angle at vertex
743  double theta = std::acos
744  ( diff_uv*diff_wv/std::sqrt((diff_wv*diff_wv)*(diff_uv*diff_uv)));
745  std::cout << "theta = " << theta << std::endl;
746 
747  if ( theta>half_pi) theta = std::max(0.,pi-theta);
748  std::cout << "theta 2 = " << theta << std::endl;
749 
750  for ( int k=cranks[vindex]-1; k>=0; --k) {
751  cs[vindex][k] += theta*(es[2*vindex+k]*(diff_uv+diff_wv));
752  // std::cout << "cs[" << vindex << "][" << k << "] += "
753  // << theta << "*(" << es[2*vindex+k] << " * ("
754  // << diff_uv << " + " << diff_wv << "))\n";
755  }
756  ws[vindex] += theta;
757  // std::cout << "ws[" << vindex << "] = " << ws[vindex] << std::endl;
758  }
759 
760  uindex=vindex; vindex=windex;
761  }
762  }
763  }
764  _wm->reduce_on_shared_nodes( c_attr, Window_manifold_2::OP_SUM);
765  _wm->reduce_on_shared_nodes( w_attr, Window_manifold_2::OP_SUM);
766 
767  // Loop through all vertices to compute center of stars.
768  it = allpanes.begin();
769  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
770  COM::Pane *pane = *it;
771 
772  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
773  ( pane->attribute( c_attr->id())->pointer());
774  double *ws = reinterpret_cast<double*>
775  ( pane->attribute( w_attr->id())->pointer());
776  int *tranks = reinterpret_cast<int*>
777  ( pane->attribute(tangranks_id)->pointer());
778  Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
779  ( pane->attribute(cntnvecs_id)->pointer());
780  int *cranks = reinterpret_cast<int*>
781  ( pane->attribute(cntnranks_id)->pointer());
782 
783  // Loop through all real nodes of the pane
784  for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
785  if ( tranks[j]==2) {
786  Vector_3<double> s(0,0,0);
787  for ( int k=cranks[j]-1; k>=0; --k) {
788  s += (cs[j][k]*es[2*j+k])/ws[j];
789  }
790  cs[j] = s;
791  }
792  }
793  }
794 
795  // Loop through the panes and its real faces to compute safe factors.
796  get_redist_safe_factor( c_attr, w_attr, 2);
797 
798  // Loop through all vertices to update displacement vector
799  it = allpanes.begin();
800  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
801  COM::Pane *pane = *it;
802 
803  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
804  ( pane->attribute( c_attr->id())->pointer());
805  double *ws = reinterpret_cast<double*>
806  ( pane->attribute( w_attr->id())->pointer());
807  Vector_3<double> *ds = reinterpret_cast<Vector_3<double>*>
808  ( pane->attribute( disps_id)->pointer());
809  int *tranks = reinterpret_cast<int*>
810  ( pane->attribute( tangranks_id)->pointer());
811 
812  // Loop through all real nodes of the pane
813  for ( int j=0, nj=pane->size_of_real_nodes(); j<nj; ++j) {
814  if ( tranks[j]==2) {
815  double w=0.5;
816  if ( ws[j]>0) w = std::min( w, 1./ws[j]);
817  ds[j] += w*cs[j];
818  std::cout << "ds[" << j << "] += " << w << " * " << cs[j] << std::endl;
819  }
820  }
821  }
822 }
823 
824 void Rocmop::
825 get_redist_safe_factor( COM::Attribute *c_attr, COM::Attribute *w_attr,
826  int rank) {
827  double zero=0;
828  Rocblas::copy_scalar( &zero, w_attr);
829 
830  const std::string att1("tangranks"), att2("weights2"),
831  att3("barycrds"), att4("PNelemids");
832  int tangranks_id = _wrk_window->attribute(att1)->id(),
833  weights2_id = _wrk_window->attribute(att2)->id(),
834  barycrds_id = _wrk_window->attribute(att3)->id(),
835  PNelemid_id = _wrk_window->attribute(att4)->id();
836 
837  std::vector< COM::Pane*> allpanes;
838  _wrk_window->panes(allpanes);
839  std::vector< COM::Pane*>::iterator it = allpanes.begin();
840  for (int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
841  COM::Pane *pane = *it;
842 
843  const Point_3<double> *pnts = reinterpret_cast<Point_3<double>*>
844  (pane->attribute(COM_NC)->pointer());
845  Vector_3<double> *cs = reinterpret_cast<Vector_3<double>*>
846  ( pane->attribute(c_attr->id())->pointer());
847  double *ws = reinterpret_cast<double*>
848  ( pane->attribute(w_attr->id())->pointer());
849  int *tranks = reinterpret_cast<int*>
850  ( pane->attribute(tangranks_id)->pointer());
851  double *ws2 = reinterpret_cast<double*>
852  ( pane->attribute(weights2_id)->pointer());
853  double *bcs = reinterpret_cast<double*>
854  ( pane->attribute(barycrds_id)->pointer());
855  int *PNids = reinterpret_cast<int*>
856  ( pane->attribute(PNelemid_id)->pointer());
857 
858  // Loop through real elements of the current pane
859  Element_node_enumerator ene( pane, 1);
860 
861  for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
862  int ne = ene.size_of_edges();
863  int uindex=ene[ne-1]-1, vindex=ene[0]-1;
864  Vector_3<double> nrm_nz;
865 
866  // Loop through all vertices of current face.
867  for ( int k=0; k<ne; ++k) {
868  int windex = ene[(k+1==ne)?0:k+1]-1;
869 
870  // Only redistribute correct classification of vertices
871  if ( tranks[vindex]==rank) {
872  Vector_3<double> tngs[2] = {pnts[uindex]-pnts[vindex],
873  pnts[windex]-pnts[vindex]};
874 
875  // Compute normal of face offset
876  if ( k==0 || ne>3) {
877  nrm_nz = Vector_3<double>::cross_product( tngs[1], tngs[0]);
878  nrm_nz.normalize();
879  }
880 
881  // Solve for vector x, to determine maximum time step
882  Vector_3<double> S[3] = { tngs[1], tngs[0], nrm_nz};
883  Vector_3<double> xs; solve( S, cs[windex], xs);
884 
885  double temp = ws[windex];
886  ws[windex] = std::max( (6-ne)*std::max( xs[0], xs[1]), ws[windex]);
887  if(temp != ws[windex]){
888  ws2[windex] = ws[windex];
889  bcs[windex*2]= xs[0];
890  bcs[windex*2+1] = xs[1];
891  PNids[windex] = j+1;
892  }
893  }
894 
895  uindex=vindex; vindex=windex;
896  }
897  }
898  }
899 
900  _wm->reduce_on_shared_nodes( w_attr, Window_manifold_2::OP_MAX);
901 }
902 
903 void Rocmop::evaluate_face_normals() {
904 
905  if(_verb > 1)
906  std::cout << "Entering Rocmop::evaluate_face_normals" << std::endl;
907 
908  COM_assertion_msg(_wrk_window, "Unexpected NULL pointer encountered.");
909 
910  if(_verb > 3)
911  std::cout << "Getting attribute ids" << std::endl;
912 
913  const std::string att1("facenormals");
914  const std::string att2("facecenters");
915  int facenormals_id = _wrk_window->attribute(att1)->id();
916  int facecenters_id = _wrk_window->attribute(att2)->id();
917 
918  // Loop through the panes and its real faces
919  std::vector< COM::Pane*> allpanes;
920  _wrk_window->panes(allpanes);
921 
922  std::vector< COM::Pane*>::iterator it = allpanes.begin();
923  for ( int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
924  COM::Pane *pane = *it;
925 
926  // Obtain nodal coordinates of current pane, assuming contiguous layout
927  COM_assertion_msg( pane->attribute( COM_NC)->stride()==3 ||
928  pane->size_of_real_nodes()==0,
929  "Coordinates must be stored in contiguous layout.");
930 
931  const Vector_3<double> *pnts = reinterpret_cast< Vector_3<double>* >
932  (pane->attribute(COM_NC)->pointer());
933  Vector_3<double> *nrms = reinterpret_cast< Vector_3<double>* >
934  ( pane->attribute(facenormals_id)->pointer());
935  Vector_3<double> *cnts = reinterpret_cast< Vector_3<double>* >
936  ( pane->attribute(facecenters_id)->pointer());
937 
938  COM_assertion_msg(pnts, "NULL pointer to nodal coords.");
939  COM_assertion_msg(nrms, "NULL pointer to face normals.");
940  COM_assertion_msg(cnts, "NULL pointer to face centers.");
941 
942  // Loop through real elements of the current pane
943  Element_node_enumerator ene( pane, 1);
944  for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
945 
946  // Obtain current face normal
947  Element_node_vectors_k_const<Vector_3<double> > ps; ps.set( pnts, ene, 1);
948  SURF::Generic_element_2 e(ene.size_of_edges(), ene.size_of_nodes());
949 
950  Vector_2<double> nc(0.5,0.5);
951  Vector_3<double> J[2];
952 
953  e.Jacobian( ps, nc, J);
954 
955  nrms[j] = Vector_3<double>::cross_product( J[0], J[1]);
956 
957  nrms[j].normalize();
958 
959  e.interpolate_to_center( ps, &cnts[j]);
960 
961  }
962  }
963 
964  if(_verb > 1)
965  std::cout << "Exiting Rocmop::evaluate_face_normals" << std::endl;
966 }
967 
968 void Rocmop::
969 identify_ridge_edges()
970 {
971  // Allocate buffer space for maximum and minedgev
972  COM::Attribute *maxedgev =
973  _wrk_window->new_attribute( "maxedgev", 'n', COM_DOUBLE, 1, "");
974  _wrk_window->resize_array( maxedgev, 0);
975 
976  COM::Attribute *minedgev =
977  _wrk_window->new_attribute( "minedgev", 'n', COM_DOUBLE, 1, "");
978  _wrk_window->resize_array( minedgev, 0);
979  _wrk_window->init_done();
980 
981  // Loop through the panes and its real faces
982  std::vector< COM::Pane*> allpanes;
983  _wrk_window->panes(allpanes);
984  std::vector< COM::Pane*>::iterator it = allpanes.begin();
985  Window_manifold_2::PM_iterator pm_it=_wm->pm_begin();
986 
987  const std::string att1("eigvecs"), att2("tangranks");
988 
989  int eigvecs_id = _wrk_window->attribute(att1)->id();
990  int tangranks_id = _wrk_window->attribute(att2)->id();
991 
992  std::vector<std::map<Edge_ID, double> > edge_maps( allpanes.size());
993  for ( int i=0, local_npanes = allpanes.size();
994  i<local_npanes; ++i, ++it, ++pm_it) {
995  COM::Pane *pane = *it;
996  std::map< Edge_ID, double> &edges_pn = edge_maps[i];
997 
998  // Obtain pointers to A_attr, bm_attr, and r_attr
999  Vector_3<double> *es = reinterpret_cast<Vector_3<double>*>
1000  ( pane->attribute(eigvecs_id)->pointer());
1001  int *tranks = reinterpret_cast<int*>
1002  ( pane->attribute(tangranks_id)->pointer());
1003  double *minvs = reinterpret_cast<double*>
1004  ( pane->attribute(minedgev->id())->pointer());
1005  double *maxvs = reinterpret_cast<double*>
1006  ( pane->attribute(maxedgev->id())->pointer());
1007 
1008  // Loop through real elements of the current pane
1009  Element_node_enumerator ene( pane, 1);
1010  for ( int j=0, nj=pane->size_of_real_elements(); j<nj; ++j, ene.next()) {
1011  Halfedge h( &*pm_it, Edge_ID(j+1,0), SURF::REAL_PANE), h0=h;
1012 
1013  int vindex = ene[ h.id().lid()]-1;
1014  do {
1015  int windex = ene[ h.next().id().lid()]-1;
1016 
1017  // If vertex is on a ridge
1018  if ( tranks[ vindex]==1) {
1019  int m = 1 + (tranks[windex]<2);
1020  Vector_3<double> tng = h.tangent(); tng.normalize();
1021  double feg = m * (es[3*vindex+2]* tng);
1022 
1023  edges_pn[ h.id()] = feg;
1024  if ( feg > maxvs[vindex]) maxvs[vindex] = feg;
1025  if ( feg < minvs[vindex]) minvs[vindex] = feg;
1026  }
1027 
1028  // Process border edge
1029  if ( h.is_border() && tranks[windex]==1) {
1030  int m = 1 + (tranks[vindex]<2);
1031  Vector_3<double> tng = -h.tangent(); tng.normalize();
1032  double feg = m * (es[3*windex+2]* tng);
1033 
1034  edges_pn[ h.opposite().id()] = feg;
1035  if ( feg > maxvs[windex]) maxvs[windex] = feg;
1036  if ( feg < minvs[windex]) minvs[windex] = feg;
1037  }
1038 
1039  // Update vindex
1040  vindex = windex;
1041  } while ( (h=h.next()) != h0);
1042  }
1043  }
1044 
1045  // Reduce on shared nodes for minedgev and maxedgev
1046  _wm->reduce_on_shared_nodes( minedgev, Window_manifold_2::OP_MIN);
1047  _wm->reduce_on_shared_nodes( maxedgev, Window_manifold_2::OP_MAX);
1048 
1049  it = allpanes.begin();
1050  // Insert ridge edges onto _edges
1051  _edges.clear(); _edges.resize( allpanes.size());
1052  for ( int i=0, local_npanes = allpanes.size(); i<local_npanes; ++i, ++it) {
1053  COM::Pane *pane = *it;
1054  std::map< Edge_ID, double> &emap_pn = edge_maps[i];
1055  std::set< Edge_ID> &eset_pn = _edges[i];
1056 
1057  const double *minvs = reinterpret_cast<double*>
1058  ( pane->attribute(minedgev->id())->pointer());
1059  const double *maxvs = reinterpret_cast<double*>
1060  ( pane->attribute(maxedgev->id())->pointer());
1061 
1062  // Loop throug ebuf and put ridges edges onto edges_pn
1063  for ( std::map< Edge_ID, double>::const_iterator eit=emap_pn.begin(),
1064  eend=emap_pn.end(); eit!=eend; ++eit) {
1065  Element_node_enumerator ene( pane, eit->first.eid());
1066  int vindex = ene[eit->first.lid()]-1;
1067 
1068  // Push edge onto edge
1069  if ( eit->second == minvs[vindex] || eit->second == maxvs[vindex])
1070  eset_pn.insert( eit->first);
1071  }
1072  }
1073 
1074  // Deallocate minedgev and maxedgev
1075  _wrk_window->delete_attribute( minedgev->name());
1076  _wrk_window->delete_attribute( maxedgev->name());
1077  _wrk_window->init_done();
1078 }
1079 
1080 int Rocmop::
1081 eigen_analyze_vertex( Vector_3<double> A_io[3], Vector_3<double> &b_io,
1082  Vector_3<double> *nrm_nz, double ad) {
1083 
1084  // A_io will be replaced by its eigenvectors.
1085  Vector_3<double> *es = A_io;
1086 
1087  // Create a backup of b_io, as b_io will be replaced by eigenvalues.
1088  const Vector_3<double> b = b_io;
1089  Vector_3<double> &lambdas = b_io;
1090 
1091  // Compute the eigenvalues and eigenvectors of the matrix.
1092  // Note that lambdas will be in decreasing order!
1093  compute_eigenvectors( es, lambdas);
1094 
1095  int nrank;
1096  double eps = lambdas[0]*1.e-8;
1097  double gs[] = { b*es[0]/lambdas[0],
1098  b*es[1]/std::max(eps,lambdas[1]),
1099  b*es[2]/std::max(eps,lambdas[2]) };
1100  double gs_abs[] = {std::abs(gs[0]), std::abs(gs[1]), std::abs(gs[2])};
1101 
1102  // Classify offset intersection based on acos(-b.norm()/rv) and lambdas
1103  if ( lambdas[2] >= _saliency_crn *
1104  std::max(lambdas[0]-lambdas[1], lambdas[1]-lambdas[2]) || ad >= .5*pi)
1105  nrank = 3;
1106  else if ( gs_abs[1] >= gs_abs[0] || lambdas[0]*_dir_thres<lambdas[1]) {
1107  if ( lambdas[1] < lambdas[0]*_eig_thres) {
1108  lambdas[1] = lambdas[0]*_eig_thres;
1109  std::cerr << "Rocprop Warning: Mesh contains near cusp."
1110  << "lambdas[1]=" << lambdas[1] << " and lambdas[0]="
1111  << lambdas[0] << std::endl;
1112  }
1113  nrank = 2; // ridge vertex
1114 #if PRINT_FEATURE
1115  if ( gs_abs[1] < gs_abs[0]) {
1116  double theta = (360/pi)*std::atan(std::sqrt(lambdas[1]/lambdas[0]));
1117  if ( theta<fangle_min_eigen) fangle_min_eigen=theta;
1118  if ( theta>fangle_max_eigen) fangle_max_eigen=theta;
1119  }
1120 #endif
1121  }
1122  else
1123  nrank = 1; // smooth vertex
1124 
1125  if ( !_reorthog || nrank==1)
1126  *nrm_nz=( es[0]*b>0)?-es[0]:es[0];
1127  else {
1128  // Solve for x within primary space
1129  Vector_3<double> x(0,0,0);
1130  for ( int k=0; k<nrank; ++k) x -= gs[k]*es[k];
1131  *nrm_nz = x.normalize();
1132  }
1133  return nrank;
1134 }
1135 
1136 #if 0
1137 // Solve the minimization problem (x^T)Ax+2b^Tx with eigen-decomposition.
1138 int Rocmop::
1139 eigen_analyze_vertex( Vector_3<double> A_io[3], Vector_3<double> &b_io,
1140  Vector_3<double> *nrm_nz) {
1141 
1142  // A_io will be replaced by its eigenvectors.
1143  Vector_3<double> *es = A_io;
1144 
1145  // Create a backup of b_io, as b_io will be replaced by eigenvalues.
1146  const Vector_3<double> b = b_io;
1147  Vector_3<double> &lambdas = b_io;
1148 
1149  // Compute the eigenvalues and eigenvectors of the matrix.
1150  // Note that lambdas will be in decreasing order!
1151  compute_eigenvectors( es, lambdas);
1152 
1153  int orank;
1154  double eps = lambdas[0]*1.e-8;
1155  double gs[] = { b*es[0]/lambdas[0],
1156  b*es[1]/std::max(eps,lambdas[1]),
1157  b*es[2]/std::max(eps,lambdas[2]) };
1158  double gs_abs[] = {std::abs(gs[0]), std::abs(gs[1]), std::abs(gs[2])};
1159 
1160  // Classify offset intersection based on acos(-b.norm()/rv) and lambdas
1161  if ( lambdas[2] >= _saliency_crn*
1162  std::max(lambdas[0]-lambdas[1], lambdas[1]-lambdas[2]) ||
1163  gs_abs[2]>=gs_abs[1] && gs_abs[2] >= gs_abs[0])
1164  orank = 3;
1165  else if ( gs_abs[1] >= gs_abs[0] || lambdas[0]*_dir_thres<lambdas[1]) {
1166  if ( lambdas[1] < lambdas[0]*_eig_thres) {
1167  lambdas[1] = lambdas[0]*_eig_thres;
1168  std::cerr << "Rocmop Warning: Mesh contains near cusp." << std::endl;
1169  }
1170  orank = 2; // ridge vertex
1171  }
1172  else
1173  orank = 1; // smooth vertex
1174 
1175  if ( !_reorthog || orank==1)
1176  *nrm_nz=( es[0]*b>0)?-es[0]:es[0];
1177  else {
1178  // Solve for x within primary space
1179  Vector_3<double> x(0,0,0);
1180  for ( int k=0; k<orank; ++k) x -= gs[k]*es[k];
1181  *nrm_nz = x.normalize();
1182  }
1183  return orank;
1184 }
1185 #endif
1186 
1187 // Compute eigenvalues and eigenvectors of a 3x3 matrix A.
1188 // At output, the eigenvalues are saved in lambdas in decending order.
1189 // The columns of A are replaced by the orthonormal eigenvectors.
1190 void Rocmop::
1191 compute_eigenvectors( Vector_3<double> A[3], Vector_3<double> &lambdas) {
1192 
1193  double abuf[3][3] = { {A[0][0], A[0][1], A[0][2]},
1194  {A[1][0], A[1][1], A[1][2]},
1195  {A[2][0], A[2][1], A[2][2]}};
1196  double ebuf[3][3];
1197 
1198  int info = dsyevq3( abuf, ebuf, &lambdas[0]);
1199  COM_assertion_msg( info==0, "Computation of eigenvectos failed");
1200 
1201  std::swap( ebuf[0][1], ebuf[1][0]);
1202  std::swap( ebuf[0][2], ebuf[2][0]);
1203  std::swap( ebuf[1][2], ebuf[2][1]);
1204 
1205  const Vector_3<double> *buf = (const Vector_3<double>*)&ebuf[0][0];
1206 
1207  if ( lambdas[0]<lambdas[1]) {
1208  if ( lambdas[1]<lambdas[2]) { // l2>l1>l0
1209  std::swap( lambdas[0], lambdas[2]);
1210  A[0] = buf[2]; A[1] = buf[1]; A[2] = buf[0];
1211  }
1212  else {
1213  double t = lambdas[0];
1214  lambdas[0] = lambdas[1];
1215  A[0] = buf[1];
1216 
1217  if ( t<lambdas[2]) { // l1>=l2>l0
1218  lambdas[1] = lambdas[2]; lambdas[2] = t;
1219  A[1] = buf[2]; A[2] = buf[0];
1220  }
1221  else { // l1>l0>=l2
1222  lambdas[1] = t;
1223  A[1] = buf[0]; A[2] = buf[2];
1224  }
1225  }
1226  }
1227  else if ( lambdas[0]<lambdas[2]) { // l2>l0>=l1
1228  double t = lambdas[0];
1229  lambdas[0] = lambdas[2]; lambdas[2] = lambdas[1]; lambdas[1] = t;
1230  A[0] = buf[2]; A[1] = buf[0]; A[2] = buf[1];
1231  }
1232  else {
1233  A[0] = buf[0];
1234  if ( lambdas[1]<lambdas[2]) { // l0>=l2>l1
1235  std::swap( lambdas[1], lambdas[2]);
1236  A[1] = buf[2]; A[2] = buf[1];
1237  }
1238  else {// l0>=l1>=l2
1239  A[1] = buf[1]; A[2] = buf[2];
1240  }
1241  }
1242 }
1243 
1244 void Rocmop::
1245 get_constraint_directions( int type, const Vector_3<double> &dir,
1246  int &ndirs, Vector_3<double> dirs[2]) {
1247  if ( type == 1 || type == -1)
1248  COM_assertion_msg( &dir && std::abs(dir.squared_norm()-1)<=1.e-6,
1249  "Constrained direction must be normalized");
1250 
1251  switch ( type) {
1252  case 0:
1253  ndirs = 0;
1254  break;
1255  case 2:
1256  ndirs = 3;
1257  break;
1258  case 1:
1259  ndirs = 1; dirs[0] = dir;
1260  break;
1261  case -1: {
1262  ndirs = 2;
1263 
1264  double sqnrm;
1265  if ( std::abs( dir[2])>0.7) {
1266  sqnrm = dir[1]*dir[1]+dir[2]*dir[2];
1267  dirs[0] = Vector_3<double>( 0, -dir[2], dir[1]);
1268  dirs[1] = Vector_3<double>( sqnrm, -dir[0]*dir[1], -dir[0]*dir[2]);
1269  }
1270  else {
1271  sqnrm = dir[0]*dir[0]+dir[1]*dir[1];
1272  dirs[0] = Vector_3<double>( dir[1], -dir[0], 0);
1273  dirs[1] = Vector_3<double>( -dir[0]*dir[2], -dir[1]*dir[2], sqnrm);
1274  }
1275 
1276  double s = 1./std::sqrt( sqnrm);
1277  dirs[0] *= s; dirs[1] *= s;
1278  break;
1279  }
1280  case 'x':
1281  ndirs = 1;
1282  dirs[0] = Vector_3<double>(1,0,0);
1283  break;
1284  case -'x':
1285  ndirs = 2;
1286  dirs[0] = Vector_3<double>(0,1,0);
1287  dirs[1] = Vector_3<double>(0,0,1);
1288  break;
1289  case 'y':
1290  ndirs = 1;
1291  dirs[0] = Vector_3<double>(0,1,0);
1292  break;
1293  case -'y':
1294  ndirs = 2;
1295  dirs[0] = Vector_3<double>(1,0,0);
1296  dirs[1] = Vector_3<double>(0,0,1);
1297  break;
1298  case 'z':
1299  ndirs = 1;
1300  dirs[0] = Vector_3<double>(0,0,1);
1301  break;
1302  case -'z':
1303  ndirs = 2;
1304  dirs[0] = Vector_3<double>(1,0,0);
1305  dirs[1] = Vector_3<double>(0,1,0);
1306  break;
1307  default: COM_assertion_msg( false, "Unknown type of constraint");
1308  }
1309 }
1310 
1311 MOP_END_NAMESPACE
1312 
1313 
1314 
1315 
1316 
1317 
Rocmop::eigen_analyze_vertex
int eigen_analyze_vertex(Vector_3< double > A_io[3], Vector_3< double > &b_io, Vector_3< double > *nrm_nz, double angle_defect)
Definition: smooth_medial.C:1081
swap
void swap(int &a, int &b)
Definition: buildface.cpp:88
Window_manifold_2::OP_SUM
static MPI_Op OP_SUM
Definition: Manifold_2.h:353
rs
subroutine rs(nm, n, a, w, matz, z, fv1, fv2, ierr)
Definition: arruda_boyce.f90:259
Element_node_enumerator
An adaptor for enumerating node IDs of an element.
Definition: Element_accessors.h:43
d
const NT & d
Definition: rational_rotation.h:73
Rocmop::_usr_window
const COM::Window * _usr_window
The user&#39;s window.
Definition: Rocmop.h:438
Pane.h
Contains the prototypes for the Pane object.
half_pi
static const double half_pi
Definition: smooth_medial.C:42
k
j indices k indices k
Definition: Indexing.h:6
Rocmop::_wght_scheme
char _wght_scheme
Weighting scheme.
Definition: Rocmop.h:419
COM_assertion_msg
#define COM_assertion_msg(EX, msg)
Definition: roccom_assertion.h:64
s
double s
Definition: blastest.C:80
Rocmop::identify_ridge_edges
void identify_ridge_edges()
Identify ridge edges.
Definition: smooth_medial.C:969
Edge_ID
MAP::Facet_ID Edge_ID
Definition: Manifold_2.h:49
Point_3
Definition: mapbasic.h:147
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
roccom.h
This file contains the prototypes for Roccom API.
Halfedge
This class encapsulate a halfedge over a window manifold.
Definition: Manifold_2.h:446
Rocmop::_reorthog
bool _reorthog
Reorthogonalize?
Definition: Rocmop.h:418
Rocmop::_eig_thres
double _eig_thres
Eigenvalue thresholds.
Definition: Rocmop.h:420
Rocmop::compute_medial_quadric
void compute_medial_quadric()
Compute medial quadric for every vertex.
Definition: smooth_medial.C:286
sqrt
double sqrt(double d)
Definition: double.h:73
Vector_3::normalize
Vector_3 & normalize()
Definition: mapbasic.h:114
Rocmop::compute_eigenvectors
static void compute_eigenvectors(Vector_3< double > A[3], Vector_3< double > &lambdas)
Definition: smooth_medial.C:1191
Window_manifold_2::PM_iterator
std::vector< Pane_manifold_2 >::iterator PM_iterator
Definition: Manifold_2.h:187
pi
static const double pi
Definition: smooth_medial.C:43
Element_node_vectors_k_const
This is a helper class for accessing nodal data.
Definition: Element_accessors.h:234
Rocmop::evaluate_face_normals
void evaluate_face_normals()
Evaluate face normals (copied from FaceOffset_3.[hC].
Definition: smooth_medial.C:903
A
rational * A
Definition: vinci_lass.c:67
SQR
double SQR(double x)
Definition: smooth_medial.C:45
Rocmop::redistribute_vertices_smooth
void redistribute_vertices_smooth()
Redistribute smooth vertices within their tangent spaces.
Definition: smooth_medial.C:682
REAL_PANE
Definition: Manifold_2.h:52
sin
NT & sin
Definition: rational_rotation.h:78
Window_manifold_2::OP_MIN
static MPI_Op OP_MIN
Definition: Manifold_2.h:353
Rocmop::get_constraint_directions
static void get_constraint_directions(int type, const Vector_3< double > &dir, int &ndirs, Vector_3< double > dirs[2])
Get orthonormals of the constraints.
Definition: smooth_medial.C:1245
Element_node_enumerator::size_of_edges
int size_of_edges() const
Number of edges per element.
Definition: Element_accessors.h:91
i
blockLoc i
Definition: read.cpp:79
angle
double angle(Vector_3< double > v1, Vector_3< double > v2)
Compute the angle between two vectors.
Definition: geometry.C:61
Rocmop::get_redist_safe_factor
void get_redist_safe_factor(COM::Attribute *c_attr, COM::Attribute *w_attr, int rank)
Definition: smooth_medial.C:825
x
void int int REAL * x
Definition: read.cpp:74
Rocmop::_rediter
int _rediter
No.iterations for vertex redistribution.
Definition: Rocmop.h:423
Rocmop::_wrk_window
COM::Window * _wrk_window
The working window.
Definition: Rocmop.h:439
Element_node_enumerator::size_of_nodes
int size_of_nodes() const
Number of nodes per element.
Definition: Element_accessors.h:74
Vector_3::squared_norm
Type squared_norm() const
Definition: mapbasic.h:110
Vector_2
Definition: mapbasic.h:199
MOP_END_NAMESPACE
#define MOP_END_NAMESPACE
Definition: mopbasic.h:29
Vector_3::cross_product
static Vector_3 cross_product(const Vector_3 &v, const Vector_3 &w)
Definition: mapbasic.h:104
cimg_library::atan
CImg< _cimg_Tfloat > atan(const CImg< T > &instance)
Definition: CImg.h:6061
Rocmop::redistribute_vertices_ridge
void redistribute_vertices_ridge()
Redistribute ridge vertices within their tangent spaces.
Definition: smooth_medial.C:547
Rocmop::_saliency_crn
double _saliency_crn
Definition: Rocmop.h:422
Rocblas.h
Definition for Rocblas API.
Generic_element_2.h
Rocmop::_verb
int _verb
Verbose level.
Definition: Rocmop.h:453
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
Rocmop::smooth_surf_medial
void smooth_surf_medial()
Smooths a surface using the medial quadric.
Definition: smooth_medial.C:230
Rocmop::_dir_thres
double _dir_thres
Another threshold.
Definition: Rocmop.h:421
Rocblas::copy_scalar
static void copy_scalar(const void *x, Attribute *y)
Operation wrapper for copy (x is a scalar pointer).
Definition: op2args.C:583
Rocmop::solve
static void solve(const FT &a1, const FT &a2, const FT &b1, const FT &b2, const FT &c1, const FT &c2, FT &x, FT &y)
Definition: Rocmop.h:378
j
j indices j
Definition: Indexing.h:6
Rocmop::_cnstr_dirs
COM::Attribute * _cnstr_dirs
Stores directions of nodal contraints.
Definition: Rocmop.h:425
q
NT q
Definition: rational_rotation.h:105
Rocmop::_cnstr_types
const COM::Attribute * _cnstr_types
Stores types of nodal constraints.
Definition: Rocmop.h:424
MOP_BEGIN_NAMESPACE
#define MOP_BEGIN_NAMESPACE
Definition: mopbasic.h:28
Element_node_vectors_k_const::set
void set(const Value *p, Element_node_enumerator &ene, int strd)
initialize the accessor with a pointer and a specific stride.
Definition: Element_accessors.h:244
NTS::abs
NT abs(const NT &x)
Definition: number_utils.h:130
dsyevq3
int dsyevq3(double A[3][3], double Q[3][3], double w[3])
Definition: smooth_medial.C:124
Rocblas::copy
static void copy(const Attribute *x, Attribute *y)
Wrapper for copy.
Definition: op2args.C:333
cimg_library::acos
CImg< _cimg_Tfloat > acos(const CImg< T > &instance)
Definition: CImg.h:6051
Element_node_enumerator::next
void next()
Go to the next element within the connectivity tables of a pane.
Definition: Element_accessors.C:79
ni
void int int REAL REAL REAL *z blockDim dim * ni
Definition: read.cpp:77
nj
void int * nj
Definition: read.cpp:74
Vector_3< double >
rank
static int rank
Definition: advectest.C:66
geometry.h
Geometric helper function header file.
dsytrd3
void dsytrd3(double A[3][3], double Q[3][3], double d[3], double e[2])
Definition: smooth_medial.C:48
Window_manifold_2::OP_MAX
static MPI_Op OP_MAX
Definition: Manifold_2.h:353
COM_DOUBLE
Definition: roccom_basic.h:131
cimg_library::cimg::info
void info()
Print informations about CImg environement variables.
Definition: CImg.h:5702
Vector_3::norm
Type norm() const
Definition: mapbasic.h:112
Rocmop::_edges
std::vector< std::set< Edge_ID > > _edges
ridge edges
Definition: Rocmop.h:426