Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
inks/IdealWeightInverseMeanRatio.cpp
Go to the documentation of this file.
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  diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
24  pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
25 
26  ***************************************************************** */
37 #include "IdealWeightInverseMeanRatio.hpp"
38 #include "MeanRatioFunctions.hpp"
39 #include "MsqTimer.hpp"
40 #include "MsqDebug.hpp"
41 
42 #ifdef MSQ_USE_OLD_STD_HEADERS
43 # include <vector.h>
44 #else
45 # include <vector>
46  using std::vector;
47 #endif
48 
49 #include <math.h>
50 
51 namespace Mesquite {
52 
53 IdealWeightInverseMeanRatio::IdealWeightInverseMeanRatio( MsqError& err, double pow_dbl )
54 {
55  set_metric_type(ELEMENT_BASED);
56  set_element_evaluation_mode(ELEMENT_VERTICES, err); MSQ_ERRRTN(err);
57 
58  set_gradient_type(ANALYTICAL_GRADIENT);
59  set_hessian_type(ANALYTICAL_HESSIAN);
60  avgMethod=QualityMetric::LINEAR;
61  feasible=1;
62  set_name("Inverse Mean Ratio");
63 
64  set_metric_power(pow_dbl, err); MSQ_ERRRTN(err);
65 }
66 
68 void IdealWeightInverseMeanRatio::set_metric_power(double pow_dbl, MsqError& err)
69 {
70  if(fabs(pow_dbl)<=MSQ_MIN){
71  MSQ_SETERR(err)(MsqError::INVALID_ARG);
72  return;
73  }
74  if(pow_dbl<0)
75  set_negate_flag(-1);
76  else
77  set_negate_flag(1);
78  a2Con=pow(.5,pow_dbl);
79  b2Con=pow_dbl;
80  c2Con=-pow_dbl;
81  a3Con=pow(1.0/3.0,pow_dbl);
82  b3Con=pow_dbl;
83  c3Con=-2.0*pow_dbl/3.0;
84 }
85 
86 
87 bool IdealWeightInverseMeanRatio::evaluate_element(PatchData &pd,
88  MsqMeshEntity *e,
89  double &m,
90  MsqError &err)
91 {
92  EntityTopology topo = e->get_element_type();
93 
94  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
95  const size_t *v_i = e->get_vertex_index_array();
96 
97  Vector3D n; // Surface normal for 2D objects
98 
99  // Hex element descriptions
100  static const int locs_hex[8][4] = {{0, 1, 3, 4},
101  {1, 2, 0, 5},
102  {2, 3, 1, 6},
103  {3, 0, 2, 7},
104  {4, 7, 5, 0},
105  {5, 4, 6, 1},
106  {6, 5, 7, 2},
107  {7, 6, 4, 3}};
108 
109  const Vector3D d_con(1.0, 1.0, 1.0);
110 
111  int i;
112 
113  m = 0.0;
114  bool metric_valid = false;
115  switch(topo) {
116  case TRIANGLE:
117  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
118  n = n / n.length(); // Need unit normal
119  mCoords[0] = vertices[v_i[0]];
120  mCoords[1] = vertices[v_i[1]];
121  mCoords[2] = vertices[v_i[2]];
122  metric_valid = m_fcn_2e(m, mCoords, n, a2Con, b2Con, c2Con);
123  if (!metric_valid) return false;
124  break;
125 
126  case QUADRILATERAL:
127  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
128  n = n / n.length(); // Need unit normal
129  for (i = 0; i < 4; ++i) {
130  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
131  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
132  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
133  metric_valid = m_fcn_2i(mMetrics[i], mCoords, n,
134  a2Con, b2Con, c2Con, d_con);
135  if (!metric_valid) return false;
136  }
137  m = average_metrics(mMetrics, 4, err);
138  break;
139 
140  case TETRAHEDRON:
141  mCoords[0] = vertices[v_i[0]];
142  mCoords[1] = vertices[v_i[1]];
143  mCoords[2] = vertices[v_i[2]];
144  mCoords[3] = vertices[v_i[3]];
145  metric_valid = m_fcn_3e(m, mCoords, a3Con, b3Con, c3Con);
146  if (!metric_valid) return false;
147  break;
148 
149  case HEXAHEDRON:
150  for (i = 0; i < 8; ++i) {
151  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
152  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
153  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
154  mCoords[3] = vertices[v_i[locs_hex[i][3]]];
155  metric_valid = m_fcn_3i(mMetrics[i], mCoords,
156  a3Con, b3Con, c3Con, d_con);
157  if (!metric_valid) return false;
158  }
159  m = average_metrics(mMetrics, 8, err); MSQ_ERRZERO(err);
160  break;
161 
162  default:
163  break;
164  } // end switch over element type
165  return true;
166 }
167 
168 bool IdealWeightInverseMeanRatio::compute_element_analytical_gradient(PatchData &pd,
169  MsqMeshEntity *e,
170  MsqVertex *v[],
171  Vector3D g[],
172  int nv,
173  double &m,
174  MsqError &err)
175 {
176  EntityTopology topo = e->get_element_type();
177 
178  if (((topo == QUADRILATERAL) || (topo == HEXAHEDRON)) &&
179  ((avgMethod == MINIMUM) || (avgMethod == MAXIMUM))) {
180  MSQ_DBGOUT(1) <<
181  "Minimum and maximum not continuously differentiable.\n"
182  "Element of subdifferential will be returned.\n";
183  }
184 
185  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
186  const size_t *v_i = e->get_vertex_index_array();
187 
188  Vector3D n; // Surface normal for 2D objects
189 
190  // Hex element descriptions
191  static const int locs_hex[8][4] = {{0, 1, 3, 4},
192  {1, 2, 0, 5},
193  {2, 3, 1, 6},
194  {3, 0, 2, 7},
195  {4, 7, 5, 0},
196  {5, 4, 6, 1},
197  {6, 5, 7, 2},
198  {7, 6, 4, 3}};
199 
200  const Vector3D d_con(1.0, 1.0, 1.0);
201 
202  int i, j;
203 
204  bool metric_valid = false;
205  int vert_per_elem = 0;
206 
207  m = 0.0;
208 
209  switch(topo) {
210  case TRIANGLE:
211  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
212  n /= n.length(); // Need unit normal
213  mCoords[0] = vertices[v_i[0]];
214  mCoords[1] = vertices[v_i[1]];
215  mCoords[2] = vertices[v_i[2]];
216  if (!g_fcn_2e(m, mAccumGrad, mCoords, n, a2Con, b2Con, c2Con)) return false;
217 
218  vert_per_elem = 3;
219  break;
220 
221  case QUADRILATERAL:
222  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
223  n /= n.length(); // Need unit normal
224  for (i = 0; i < 4; ++i) {
225  mAccumGrad[i] = 0.0;
226 
227  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
228  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
229  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
230  if (!g_fcn_2i(mMetrics[i], mGradients+3*i, mCoords, n,
231  a2Con, b2Con, c2Con, d_con)) return false;
232  }
233 
234  m = average_metric_and_weights( mMetrics, 4, err ); MSQ_ERRZERO(err);
235  for (i = 0; i < 4; ++i) {
236  mAccumGrad[locs_hex[i][0]] += mMetrics[i]*mGradients[3*i+0];
237  mAccumGrad[locs_hex[i][1]] += mMetrics[i]*mGradients[3*i+1];
238  mAccumGrad[locs_hex[i][2]] += mMetrics[i]*mGradients[3*i+2];
239  }
240 
241  vert_per_elem = 4;
242  break;
243 
244  case TETRAHEDRON:
245  mCoords[0] = vertices[v_i[0]];
246  mCoords[1] = vertices[v_i[1]];
247  mCoords[2] = vertices[v_i[2]];
248  mCoords[3] = vertices[v_i[3]];
249  metric_valid = g_fcn_3e(m, mAccumGrad, mCoords, a3Con, b3Con, c3Con);
250  if (!metric_valid) return false;
251 
252  vert_per_elem = 4;
253  break;
254 
255  case HEXAHEDRON:
256  for (i = 0; i < 8; ++i) {
257  mAccumGrad[i] = 0.0;
258 
259  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
260  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
261  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
262  mCoords[3] = vertices[v_i[locs_hex[i][3]]];
263  if (!g_fcn_3i(mMetrics[i], mGradients+4*i, mCoords,
264  a3Con, b3Con, c3Con, d_con)) return false;
265  }
266 
267  m = average_metric_and_weights( mMetrics, 8, err ); MSQ_ERRZERO(err);
268  for (i = 0; i < 8; ++i) {
269  mAccumGrad[locs_hex[i][0]] += mMetrics[i]*mGradients[4*i+0];
270  mAccumGrad[locs_hex[i][1]] += mMetrics[i]*mGradients[4*i+1];
271  mAccumGrad[locs_hex[i][2]] += mMetrics[i]*mGradients[4*i+2];
272  mAccumGrad[locs_hex[i][3]] += mMetrics[i]*mGradients[4*i+3];
273  }
274 
275  vert_per_elem = 8;
276  break;
277 
278  default:
279  break;
280  } // end switch over element type
281 
282 
283  // This is not very efficient, but is one way to select correct gradients
284  // For gradients, info is returned only for free vertices, in the order of v[].
285  for (i = 0; i < vert_per_elem; ++i) {
286  for (j = 0; j < nv; ++j) {
287  if (vertices + v_i[i] == v[j]) {
288  g[j] = mAccumGrad[i];
289  }
290  }
291  }
292 
293 
294  return true;
295 }
296 
297 
298 bool IdealWeightInverseMeanRatio::compute_element_analytical_hessian(PatchData &pd,
299  MsqMeshEntity *e,
300  MsqVertex *fv[],
301  Vector3D g[],
302  Matrix3D h[],
303  int nfv,
304  double &m,
305  MsqError &err)
306 {
307  EntityTopology topo = e->get_element_type();
308 
309  if (((topo == QUADRILATERAL) || (topo == HEXAHEDRON)) &&
310  ((avgMethod == MINIMUM) || (avgMethod == MAXIMUM))) {
311  MSQ_DBGOUT(1) <<
312  "Minimum and maximum not continuously differentiable.\n"
313  "Element of subdifferential will be returned.\n"
314  "Who knows what the Hessian is?\n" ;
315  }
316 
317  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
318  const size_t *v_i = e->get_vertex_index_array();
319 
320 
321  Vector3D n; // Surface normal for 2D objects
322  Matrix3D outer;
323  double nm, t=0;
324 
325  // Hex element descriptions
326  static const int locs_hex[8][4] = {{0, 1, 3, 4},
327  {1, 2, 0, 5},
328  {2, 3, 1, 6},
329  {3, 0, 2, 7},
330  {4, 7, 5, 0},
331  {5, 4, 6, 1},
332  {6, 5, 7, 2},
333  {7, 6, 4, 3}};
334 
335  const Vector3D d_con(1.0, 1.0, 1.0);
336 
337  int i, j, k, l, ind;
338  int r, c, loc;
339 
340  bool metric_valid = false;
341 
342  m = 0.0;
343 
344  switch(topo) {
345  case TRIANGLE:
346  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
347  n /= n.length(); // Need unit normal
348  mCoords[0] = vertices[v_i[0]];
349  mCoords[1] = vertices[v_i[1]];
350  mCoords[2] = vertices[v_i[2]];
351  if (!h_fcn_2e(m, g, h, mCoords, n, a2Con, b2Con, c2Con)) return false;
352 
353  // zero out fixed elements of g
354  j = 0;
355  for (i = 0; i < 3; ++i) {
356  // if free vertex, see next
357  if (j<nfv && vertices + v_i[i] == fv[j] )
358  ++j;
359  // else zero gradient and Hessian entries
360  else {
361  g[i] = 0.;
362 
363  switch(i) {
364  case 0:
365  h[0].zero(); h[1].zero(); h[2].zero();
366  break;
367 
368  case 1:
369  h[1].zero(); h[3].zero(); h[4].zero();
370  break;
371 
372  case 2:
373  h[2].zero(); h[4].zero(); h[5].zero();
374  }
375  }
376  }
377  break;
378 
379  case QUADRILATERAL:
380  for (i=0; i < 10; ++i) {
381  h[i].zero();
382  }
383 
384  pd.get_domain_normal_at_element(e, n, err); MSQ_ERRZERO(err);
385  n /= n.length(); // Need unit normal
386  for (i = 0; i < 4; ++i) {
387  g[i] = 0.0;
388 
389  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
390  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
391  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
392  if (!h_fcn_2i(mMetrics[i], mGradients+3*i, mHessians+6*i, mCoords, n,
393  a2Con, b2Con, c2Con, d_con)) return false;
394  }
395 
396  switch(avgMethod) {
397  case MINIMUM:
398  MSQ_SETERR(err)("MINIMUM averaging method does not work.",MsqError::INVALID_STATE);
399  return false;
400 
401  case MAXIMUM:
402  MSQ_SETERR(err)("MAXIMUM averaging method does not work.",MsqError::INVALID_STATE);
403  return false;
404 
405  case SUM:
406  m = 0;
407  for (i = 0; i < 4; ++i) {
408  m += mMetrics[i];
409  }
410 
411  l = 0;
412  for (i = 0; i < 4; ++i) {
413  g[locs_hex[i][0]] += mGradients[3*i+0];
414  g[locs_hex[i][1]] += mGradients[3*i+1];
415  g[locs_hex[i][2]] += mGradients[3*i+2];
416 
417  for (j = 0; j < 3; ++j) {
418  for (k = j; k < 3; ++k) {
419  r = locs_hex[i][j];
420  c = locs_hex[i][k];
421 
422  if (r <= c) {
423  loc = 4*r - (r*(r+1)/2) + c;
424  h[loc] += mHessians[l];
425  }
426  else {
427  loc = 4*c - (c*(c+1)/2) + r;
428  h[loc] += transpose(mHessians[l]);
429  }
430  ++l;
431  }
432  }
433  }
434  break;
435 
436  case SUM_SQUARED:
437  m = 0;
438  for (i = 0; i < 4; ++i) {
439  m += (mMetrics[i]*mMetrics[i]);
440  mMetrics[i] *= 2;
441  }
442 
443  l = 0;
444  for (i = 0; i < 4; ++i) {
445  g[locs_hex[i][0]] += mMetrics[i]*mGradients[3*i+0];
446  g[locs_hex[i][1]] += mMetrics[i]*mGradients[3*i+1];
447  g[locs_hex[i][2]] += mMetrics[i]*mGradients[3*i+2];
448 
449  for (j = 0; j < 3; ++j) {
450  for (k = j; k < 3; ++k) {
451  outer = 2.0*outer.outer_product(mGradients[3*i+j],
452  mGradients[3*i+k]);
453 
454  r = locs_hex[i][j];
455  c = locs_hex[i][k];
456 
457  if (r <= c) {
458  loc = 4*r - (r*(r+1)/2) + c;
459  h[loc] += mMetrics[i]*mHessians[l] + outer;
460  }
461  else {
462  loc = 4*c - (c*(c+1)/2) + r;
463  h[loc] += transpose(mMetrics[i]*mHessians[l] + outer);
464  }
465  ++l;
466  }
467  }
468  }
469  break;
470 
471  case LINEAR:
472  m = 0;
473  for (i = 0; i < 4; ++i) {
474  m += mMetrics[i];
475  }
476  m *= 0.25;
477 
478  l = 0;
479  for (i = 0; i < 4; ++i) {
480  g[locs_hex[i][0]] += 0.25*mGradients[3*i+0];
481  g[locs_hex[i][1]] += 0.25*mGradients[3*i+1];
482  g[locs_hex[i][2]] += 0.25*mGradients[3*i+2];
483 
484  for (j = 0; j < 3; ++j) {
485  for (k = j; k < 3; ++k) {
486  r = locs_hex[i][j];
487  c = locs_hex[i][k];
488 
489  if (r <= c) {
490  loc = 4*r - (r*(r+1)/2) + c;
491  h[loc] += mHessians[l];
492  }
493  else {
494  loc = 4*c - (c*(c+1)/2) + r;
495  h[loc] += transpose(mHessians[l]);
496  }
497  ++l;
498  }
499  }
500  }
501 
502  for (i=0; i<10; ++i) {
503  h[i] *= 0.25;
504  }
505  break;
506 
507  case GEOMETRIC:
508  MSQ_SETERR(err)("GEOMETRIC averaging method does not work.",MsqError::INVALID_STATE);
509  return false;
510 
511  default:
512  switch(avgMethod) {
513  case RMS:
514  t = 2.0;
515  break;
516 
517  case HARMONIC:
518  t = -1.0;
519  break;
520 
521  case HMS:
522  t = -2.0;
523  break;
524 
525  default:
526  MSQ_SETERR(err)("averaging method not available.",MsqError::NOT_IMPLEMENTED);
527  break;
528  }
529 
530  m = 0;
531  for (i = 0; i < 4; ++i) {
532  nm = pow(mMetrics[i], t);
533  m += nm;
534 
535  g_factor[i] = 0.25*t*nm / mMetrics[i];
536  h_factor[i] = (t-1)*g_factor[i] / mMetrics[i];
537  }
538 
539  nm = 0.25 * m;
540 
541  l = 0;
542  for (i = 0; i < 4; ++i) {
543  g[locs_hex[i][0]] += g_factor[i]*mGradients[3*i+0];
544  g[locs_hex[i][1]] += g_factor[i]*mGradients[3*i+1];
545  g[locs_hex[i][2]] += g_factor[i]*mGradients[3*i+2];
546 
547  for (j = 0; j < 3; ++j) {
548  for (k = j; k < 3; ++k) {
549  outer = h_factor[i] * outer.outer_product(mGradients[3*i+j],
550  mGradients[3*i+k]);
551 
552  r = locs_hex[i][j];
553  c = locs_hex[i][k];
554 
555  if (r <= c) {
556  loc = 4*r - (r*(r+1)/2) + c;
557  h[loc] += g_factor[i]*mHessians[l] + outer;
558  }
559  else {
560  loc = 4*c - (c*(c+1)/2) + r;
561  h[loc] += transpose(g_factor[i]*mHessians[l] + outer);
562  }
563  ++l;
564  }
565  }
566  }
567 
568  m = pow(nm, 1.0 / t);
569  g_factor[0] = m / (t*nm);
570  h_factor[0] = (1.0 / t - 1)*g_factor[0] / nm;
571 
572  l = 0;
573  for (i = 0; i < 4; ++i) {
574  for (j = i; j < 4; ++j) {
575  outer = outer.outer_product(g[i], g[j]);
576  h[l] = g_factor[0]*h[l] + h_factor[0]*outer;
577  ++l;
578  }
579  g[i] *= g_factor[0];
580  }
581  break;
582  }
583 
584  // zero out fixed elements of gradient and Hessian
585  ind = 0;
586  for (i=0; i<4; ++i) {
587  // if free vertex, see next
588  if (ind<nfv && vertices+v_i[i] == fv[ind] )
589  ++ind;
590  // else zero gradient entry and hessian entries.
591  else {
592  g[i] = 0.;
593  switch(i) {
594  case 0:
595  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
596  break;
597 
598  case 1:
599  h[1].zero(); h[4].zero(); h[5].zero(); h[6].zero();
600  break;
601 
602  case 2:
603  h[2].zero(); h[5].zero(); h[7].zero(); h[8].zero();
604  break;
605 
606  case 3:
607  h[3].zero(); h[6].zero(); h[8].zero(); h[9].zero();
608  break;
609  }
610  }
611  }
612  break;
613 
614  case TETRAHEDRON:
615  mCoords[0] = vertices[v_i[0]];
616  mCoords[1] = vertices[v_i[1]];
617  mCoords[2] = vertices[v_i[2]];
618  mCoords[3] = vertices[v_i[3]];
619  metric_valid = h_fcn_3e(m, g, h, mCoords, a3Con, b3Con, c3Con);
620  if (!metric_valid) return false;
621 
622  // zero out fixed elements of g
623  j = 0;
624  for (i = 0; i < 4; ++i) {
625  // if free vertex, see next
626  if (j<nfv && vertices + v_i[i] == fv[j] )
627  ++j;
628  // else zero gradient entry
629  else {
630  g[i] = 0.;
631 
632  switch(i) {
633  case 0:
634  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
635  break;
636 
637  case 1:
638  h[1].zero(); h[4].zero(); h[5].zero(); h[6].zero();
639  break;
640 
641  case 2:
642  h[2].zero(); h[5].zero(); h[7].zero(); h[8].zero();
643  break;
644 
645  case 3:
646  h[3].zero(); h[6].zero(); h[8].zero(); h[9].zero();
647  break;
648  }
649  }
650  }
651  break;
652 
653  case HEXAHEDRON:
654  for (i=0; i<36; ++i)
655  h[i].zero();
656 
657  for (i = 0; i < 8; ++i) {
658  g[i] = 0.0;
659 
660  mCoords[0] = vertices[v_i[locs_hex[i][0]]];
661  mCoords[1] = vertices[v_i[locs_hex[i][1]]];
662  mCoords[2] = vertices[v_i[locs_hex[i][2]]];
663  mCoords[3] = vertices[v_i[locs_hex[i][3]]];
664  if (!h_fcn_3i(mMetrics[i], mGradients+4*i, mHessians+10*i, mCoords,
665  a3Con, b3Con, c3Con, d_con)) return false;
666  }
667 
668  switch(avgMethod) {
669  case MINIMUM:
670  MSQ_SETERR(err)("MINIMUM averaging method does not work.",MsqError::NOT_IMPLEMENTED);
671  return false;
672 
673  case MAXIMUM:
674  MSQ_SETERR(err)("MAXIMUM averaging method does not work.",MsqError::NOT_IMPLEMENTED);
675  return false;
676 
677  case SUM:
678  m = 0;
679  for (i = 0; i < 8; ++i) {
680  m += mMetrics[i];
681  }
682 
683  l = 0;
684  for (i = 0; i < 8; ++i) {
685  g[locs_hex[i][0]] += mGradients[4*i+0];
686  g[locs_hex[i][1]] += mGradients[4*i+1];
687  g[locs_hex[i][2]] += mGradients[4*i+2];
688  g[locs_hex[i][3]] += mGradients[4*i+3];
689 
690  for (j = 0; j < 4; ++j) {
691  for (k = j; k < 4; ++k) {
692  r = locs_hex[i][j];
693  c = locs_hex[i][k];
694 
695  if (r <= c) {
696  loc = 8*r - (r*(r+1)/2) + c;
697  h[loc] += mHessians[l];
698  }
699  else {
700  loc = 8*c - (c*(c+1)/2) + r;
701  h[loc] += transpose(mHessians[l]);
702  }
703  ++l;
704  }
705  }
706  }
707  break;
708 
709  case SUM_SQUARED:
710  m = 0;
711  for (i = 0; i < 8; ++i) {
712  m += (mMetrics[i]*mMetrics[i]);
713  mMetrics[i] *= 2.0;
714  }
715 
716  l = 0;
717  for (i = 0; i < 8; ++i) {
718  g[locs_hex[i][0]] += mMetrics[i]*mGradients[4*i+0];
719  g[locs_hex[i][1]] += mMetrics[i]*mGradients[4*i+1];
720  g[locs_hex[i][2]] += mMetrics[i]*mGradients[4*i+2];
721  g[locs_hex[i][3]] += mMetrics[i]*mGradients[4*i+3];
722 
723  for (j = 0; j < 4; ++j) {
724  for (k = j; k < 4; ++k) {
725  outer = 2.0*outer.outer_product(mGradients[4*i+j],
726  mGradients[4*i+k]);
727 
728  r = locs_hex[i][j];
729  c = locs_hex[i][k];
730 
731  if (r <= c) {
732  loc = 8*r - (r*(r+1)/2) + c;
733  h[loc] += mMetrics[i]*mHessians[l] + outer;
734  }
735  else {
736  loc = 8*c - (c*(c+1)/2) + r;
737  h[loc] += transpose(mMetrics[i]*mHessians[l] + outer);
738  }
739  ++l;
740  }
741  }
742  }
743  break;
744 
745  case LINEAR:
746  m = 0;
747  for (i = 0; i < 8; ++i) {
748  m += mMetrics[i];
749  }
750  m *= 0.125;
751 
752  l = 0;
753  for (i = 0; i < 8; ++i) {
754  g[locs_hex[i][0]] += 0.125*mGradients[4*i+0];
755  g[locs_hex[i][1]] += 0.125*mGradients[4*i+1];
756  g[locs_hex[i][2]] += 0.125*mGradients[4*i+2];
757  g[locs_hex[i][3]] += 0.125*mGradients[4*i+3];
758 
759  for (j = 0; j < 4; ++j) {
760  for (k = j; k < 4; ++k) {
761  r = locs_hex[i][j];
762  c = locs_hex[i][k];
763 
764  if (r <= c) {
765  loc = 8*r - (r*(r+1)/2) + c;
766  h[loc] += mHessians[l];
767  }
768  else {
769  loc = 8*c - (c*(c+1)/2) + r;
770  h[loc] += transpose(mHessians[l]);
771  }
772  ++l;
773  }
774  }
775  }
776 
777  for (i=0; i<36; ++i)
778  h[i] *= 0.125;
779  break;
780 
781  case GEOMETRIC:
782  MSQ_SETERR(err)("GEOMETRIC averaging method does not work.",MsqError::NOT_IMPLEMENTED);
783  return false;
784 
785  default:
786  switch(avgMethod) {
787  case RMS:
788  t = 2.0;
789  break;
790 
791  case HARMONIC:
792  t = -1.0;
793  break;
794 
795  case HMS:
796  t = -2.0;
797  break;
798 
799  default:
800  MSQ_SETERR(err)("averaging method not available.",MsqError::NOT_IMPLEMENTED);
801  break;
802  }
803 
804  m = 0;
805  for (i = 0; i < 8; ++i) {
806  nm = pow(mMetrics[i], t);
807  m += nm;
808 
809  g_factor[i] = 0.125*t*nm / mMetrics[i];
810  h_factor[i] = (t-1)*g_factor[i] / mMetrics[i];
811  }
812 
813  nm = 0.125 * m;
814 
815  l = 0;
816  for (i = 0; i < 8; ++i) {
817  g[locs_hex[i][0]] += g_factor[i]*mGradients[4*i+0];
818  g[locs_hex[i][1]] += g_factor[i]*mGradients[4*i+1];
819  g[locs_hex[i][2]] += g_factor[i]*mGradients[4*i+2];
820  g[locs_hex[i][3]] += g_factor[i]*mGradients[4*i+3];
821 
822  for (j = 0; j < 4; ++j) {
823  for (k = j; k < 4; ++k) {
824  outer = h_factor[i]*outer.outer_product(mGradients[4*i+j],
825  mGradients[4*i+k]);
826 
827  r = locs_hex[i][j];
828  c = locs_hex[i][k];
829 
830  if (r <= c) {
831  loc = 8*r - (r*(r+1)/2) + c;
832  h[loc] += g_factor[i]*mHessians[l] + outer;
833  }
834  else {
835  loc = 8*c - (c*(c+1)/2) + r;
836  h[loc] += transpose(g_factor[i]*mHessians[l] + outer);
837  }
838  ++l;
839  }
840  }
841  }
842 
843  m = pow(nm, 1.0 / t);
844  g_factor[0] = m / (t*nm);
845  h_factor[0] = (1.0 / t - 1)*g_factor[0] / nm;
846 
847  l = 0;
848  for (i = 0; i < 8; ++i) {
849  for (j = i; j < 8; ++j) {
850  outer = outer.outer_product(g[i], g[j]);
851  h[l] = g_factor[0]*h[l] + h_factor[0]*outer;
852  ++l;
853  }
854  g[i] *= g_factor[0];
855  }
856  break;
857  }
858 
859  // zero out fixed elements of gradient and Hessian
860  ind = 0;
861  for (i=0; i<8; ++i) {
862  // if free vertex, see next
863  if (ind<nfv && vertices+v_i[i] == fv[ind] )
864  ++ind;
865  // else zero gradient entry and hessian entries.
866  else {
867  g[i] = 0.;
868  switch(i) {
869  case 0:
870  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
871  h[4].zero(); h[5].zero(); h[6].zero(); h[7].zero();
872  break;
873 
874  case 1:
875  h[1].zero(); h[8].zero(); h[9].zero(); h[10].zero();
876  h[11].zero(); h[12].zero(); h[13].zero(); h[14].zero();
877  break;
878 
879  case 2:
880  h[2].zero(); h[9].zero(); h[15].zero(); h[16].zero();
881  h[17].zero(); h[18].zero(); h[19].zero(); h[20].zero();
882  break;
883 
884  case 3:
885  h[3].zero(); h[10].zero(); h[16].zero(); h[21].zero();
886  h[22].zero(); h[23].zero(); h[24].zero(); h[25].zero();
887  break;
888 
889  case 4:
890  h[4].zero(); h[11].zero(); h[17].zero(); h[22].zero();
891  h[26].zero(); h[27].zero(); h[28].zero(); h[29].zero();
892  break;
893 
894  case 5:
895  h[5].zero(); h[12].zero(); h[18].zero(); h[23].zero();
896  h[27].zero(); h[30].zero(); h[31].zero(); h[32].zero();
897  break;
898 
899  case 6:
900  h[6].zero(); h[13].zero(); h[19].zero(); h[24].zero();
901  h[28].zero(); h[31].zero(); h[33].zero(); h[34].zero();
902  break;
903 
904  case 7:
905  h[7].zero(); h[14].zero(); h[20].zero(); h[25].zero();
906  h[29].zero(); h[32].zero(); h[34].zero(); h[35].zero();
907  break;
908  }
909  }
910  }
911  break;
912 
913  default:
914  break;
915  } // end switch over element type
916 
917  return true;
918 }
919 
920 } // namespace Mesquite
Mesquite::h_fcn_3i
int h_fcn_3i(double &obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:1973
Mesquite::g_fcn_2i
bool g_fcn_2i(double &obj, Vector3D g_obj[3], const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:584
Mesquite::m_fcn_3e
bool m_fcn_3e(double &obj, const Vector3D x[4], const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:931
MSQ_ERRZERO
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
Definition: includeLinks/MsqError.hpp:72
Mesquite::g_fcn_2e
bool g_fcn_2e(double &obj, Vector3D g_obj[3], const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:170
Mesquite::g_fcn_3e
bool g_fcn_3e(double &obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:973
k
j indices k indices k
Definition: Indexing.h:6
Mesquite::m_fcn_3i
bool m_fcn_3i(double &obj, const Vector3D x[4], const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:1855
Mesquite::EntityTopology
EntityTopology
Definition: Mesquite.hpp:92
Mesquite::transpose
Matrix3D transpose(const Matrix3D &A)
Definition: includeLinks/Matrix3D.hpp:335
Mesquite::MsqError::NOT_IMPLEMENTED
requested functionality is not (yet) implemented
Definition: includeLinks/MsqError.hpp:128
Mesquite::IdealWeightInverseMeanRatio::evaluate_element
bool evaluate_element(PatchData &pd, MsqMeshEntity *element, double &fval, MsqError &err)
evaluate using mesquite objects
Definition: QualityMetric/Shape/IdealWeightInverseMeanRatio.cpp:87
Mesquite::m_fcn_2i
bool m_fcn_2i(double &obj, const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:543
Mesquite::g_fcn_3i
bool g_fcn_3i(double &obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:1895
v
*********************************************************************Illinois Open Source License ****University of Illinois NCSA **Open Source License University of Illinois All rights reserved ****Developed free of to any person **obtaining a copy of this software and associated documentation to deal with the Software without including without limitation the rights to and or **sell copies of the and to permit persons to whom the **Software is furnished to do subject to the following this list of conditions and the following disclaimers ****Redistributions in binary form must reproduce the above **copyright this list of conditions and the following **disclaimers in the documentation and or other materials **provided with the distribution ****Neither the names of the Center for Simulation of Advanced the University of nor the names of its **contributors may be used to endorse or promote products derived **from this Software without specific prior written permission ****THE SOFTWARE IS PROVIDED AS WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **ARISING OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE v
Definition: roccomf90.h:20
Mesquite::QualityMetric::average_metric_and_weights
double average_metric_and_weights(double metric_values[], int num_metric_values, MsqError &err)
Given a list of metric values, calculate the average metric valude according to the current avgMethod...
Definition: QualityMetric/QualityMetric.cpp:495
Mesquite::QualityMetric::set_hessian_type
void set_hessian_type(HESSIAN_TYPE ht)
Sets hessianType for this metric.
Definition: includeLinks/QualityMetric.hpp:224
Mesquite::IdealWeightInverseMeanRatio::compute_element_analytical_hessian
bool compute_element_analytical_hessian(PatchData &pd, MsqMeshEntity *e, MsqVertex *v[], Vector3D g[], Matrix3D h[], int nv, double &m, MsqError &err)
Definition: QualityMetric/Shape/IdealWeightInverseMeanRatio.cpp:298
Mesquite::MsqError::INVALID_ARG
invalid function argument passed
Definition: includeLinks/MsqError.hpp:119
nvc::Vector3D
NVec< 3, double > Vector3D
Definition: Rocon/include/NVec.h:219
Mesquite::h_fcn_3e
bool h_fcn_3e(double &obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:1062
Mesquite::QualityMetric::set_element_evaluation_mode
void set_element_evaluation_mode(ElementEvaluationMode mode, MsqError &err)
Sets the evaluation mode for the ELEMENT_BASED metrics.
Definition: includeLinks/QualityMetric.hpp:402
Mesquite::m_fcn_2e
bool m_fcn_2e(double &obj, const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:131
MSQ_SETERR
#define MSQ_SETERR(err)
Macro to set error - use err.clear() to clear.
Definition: includeLinks/MsqError.hpp:83
i
blockLoc i
Definition: read.cpp:79
Mesquite::IdealWeightInverseMeanRatio::compute_element_analytical_gradient
bool compute_element_analytical_gradient(PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
Virtual function that computes the gradient of the QualityMetric analytically. The base class impleme...
Definition: QualityMetric/Shape/IdealWeightInverseMeanRatio.cpp:168
n
const NT & n
Definition: rational_rotation.h:72
Mesquite::h_fcn_2e
bool h_fcn_2e(double &obj, Vector3D g_obj[3], Matrix3D h_obj[6], const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c)
Definition: includeLinks/MeanRatioFunctions.hpp:244
j
j indices j
Definition: Indexing.h:6
Mesquite::QualityMetric::set_gradient_type
void set_gradient_type(GRADIENT_TYPE grad)
Sets gradType for this metric.
Definition: includeLinks/QualityMetric.hpp:212
Mesquite::MsqError::INVALID_STATE
object is in an invalid state
Definition: includeLinks/MsqError.hpp:121
Mesquite::QualityMetric::set_metric_type
void set_metric_type(MetricType t)
This function should be used in the constructor of every concrete quality metric. ...
Definition: includeLinks/QualityMetric.hpp:291
Mesquite::QualityMetric::average_metrics
double average_metrics(const double metric_values[], const int &num_values, MsqError &err)
average_metrics takes an array of length num_values and averages the contents using averaging method ...
Definition: includeLinks/QualityMetric.hpp:524
Mesquite::IdealWeightInverseMeanRatio::set_metric_power
void set_metric_power(double pow_dbl, MsqError &err)
Sets the power value in the metric computation.
Definition: QualityMetric/Shape/IdealWeightInverseMeanRatio.cpp:68
MSQ_DBGOUT
#define MSQ_DBGOUT(flag)
Check debug flag and return ostream associated with flag.
Definition: includeLinks/MsqDebug.hpp:158
MSQ_ERRRTN
#define MSQ_ERRRTN(err)
If passed error is true, return from a void function.
Definition: includeLinks/MsqError.hpp:68
Mesquite::QualityMetric::set_name
void set_name(msq_std::string st)
Sets the name of this metric.
Definition: includeLinks/QualityMetric.hpp:181
Mesquite::h_fcn_2i
bool h_fcn_2i(double &obj, Vector3D g_obj[3], Matrix3D h_obj[6], const Vector3D x[3], const Vector3D &n, const double a, const Exponent &b, const Exponent &c, const Vector3D &d)
Definition: includeLinks/MeanRatioFunctions.hpp:656
Mesquite::pow
double pow(double value, const Exponent &exp)
Definition: includeLinks/Exponent.hpp:89
Mesquite::MSQ_MIN
const double MSQ_MIN
Definition: Mesquite.hpp:160