gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshMetric.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <algorithm>
7 #include "meshMetric.h"
8 #include "meshGFaceOptimize.h"
9 #include "Context.h"
10 #include "GEntity.h"
11 #include "GModel.h"
12 #include "gmshLevelset.h"
13 #include "MElementOctree.h"
14 #include "OS.h"
15 #include "Context.h"
16 
17 meshMetric::meshMetric(GModel *gm)
18 {
19  hasAnalyticalMetric = false;
20  _dim = gm->getDim();
21  std::map<MElement *, MElement *> newP;
22  std::map<MElement *, MElement *> newD;
23 
24  if(_dim == 2) {
25  for(auto fit = gm->firstFace(); fit != gm->lastFace(); ++fit) {
26  for(std::size_t i = 0; i < (*fit)->getNumMeshElements(); i++) {
27  MElement *e = (*fit)->getMeshElement(i);
28  MElement *copy = e->copy(_vertexMap, newP, newD);
29  _elements.push_back(copy);
30  }
31  }
32  }
33  else if(_dim == 3) {
34  for(auto rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit) {
35  for(std::size_t i = 0; i < (*rit)->getNumMeshElements(); i++) {
36  MElement *e = (*rit)->getMeshElement(i);
37  MElement *copy = e->copy(_vertexMap, newP, newD);
38  _elements.push_back(copy);
39  }
40  }
41  }
42  _octree = new MElementOctree(_elements);
43  buildVertexToElement(_elements, _adj);
44 }
45 
46 meshMetric::meshMetric(std::vector<MElement *> elements)
47 {
48  hasAnalyticalMetric = false;
49 
50  _dim = elements[0]->getDim();
51  std::map<MElement *, MElement *> newP;
52  std::map<MElement *, MElement *> newD;
53 
54  for(std::size_t i = 0; i < elements.size(); i++) {
55  MElement *e = elements[i];
56  MElement *copy = e->copy(_vertexMap, newP, newD);
57  _elements.push_back(copy);
58  }
59 
60  _octree = new MElementOctree(_elements);
61  buildVertexToElement(_elements, _adj);
62 }
63 
64 void meshMetric::addMetric(int technique, simpleFunction<double> *fct,
65  const std::vector<double> &parameters)
66 {
67  needMetricUpdate = true;
68 
69  int metricNumber = setOfMetrics.size();
70  setOfFcts[metricNumber] = fct;
71  setOfParameters[metricNumber] = parameters;
72  setOfTechniques[metricNumber] = technique;
73 
74  if(fct->hasDerivatives()) hasAnalyticalMetric = true;
75 
76  computeMetric(metricNumber);
77 }
78 
79 void meshMetric::updateMetrics()
80 {
81  if(!setOfMetrics.size()) {
82  Msg::Error("Can't intersect metrics, no metric registered");
83  return;
84  }
85 
86  auto it = _adj.begin();
87  for(; it != _adj.end(); it++) {
88  MVertex *ver = it->first;
89  _nodalMetrics[ver] = setOfMetrics[0][ver];
90  _nodalSizes[ver] = setOfSizes[0][ver];
91 
92  if(setOfMetrics.size() > 1)
93  for(std::size_t i = 1; i < setOfMetrics.size(); i++) {
94  _nodalMetrics[ver] =
95  (_dim == 3) ? intersection_conserve_mostaniso(_nodalMetrics[ver],
96  setOfMetrics[i][ver]) :
97  intersection_conserve_mostaniso_2d(
98  _nodalMetrics[ver], setOfMetrics[i][ver]);
99  _nodalSizes[ver] = std::min(_nodalSizes[ver], setOfSizes[i][ver]);
100  }
101  }
102  needMetricUpdate = false;
103 }
104 
105 void meshMetric::exportInfo(const char *fileendname)
106 {
107  if(needMetricUpdate) updateMetrics();
108  std::stringstream sg, sm, sl, sh, shm;
109  sg << "meshmetric_gradients_" << fileendname;
110  sm << "meshmetric_metric_" << fileendname;
111  sl << "meshmetric_levelset_" << fileendname;
112  sh << "meshmetric_hessian_" << fileendname;
113  shm << "meshmetric_hessianmat_" << fileendname;
114  std::ofstream out_grad(sg.str().c_str());
115  std::ofstream out_metric(sm.str().c_str());
116  std::ofstream out_ls(sl.str().c_str());
117  std::ofstream out_hess(sh.str().c_str());
118  std::ofstream out_hessmat(shm.str().c_str());
119  out_grad << "View \"ls_gradient\"{" << std::endl;
120  out_metric << "View \"metric\"{" << std::endl;
121  out_ls << "View \"ls\"{" << std::endl;
122  out_hess << "View \"hessian\"{" << std::endl;
123  out_hessmat << "View \"hessian_mat\"{" << std::endl;
124  auto itelem = _elements.begin();
125  auto itelemen = _elements.end();
126  for(; itelem != itelemen; itelem++) {
127  MElement *e = *itelem;
128  if(e->getDim() == 2) {
129  out_metric << "TT(";
130  out_grad << "VT(";
131  out_ls << "ST(";
132  out_hess << "ST(";
133  out_hessmat << "TT(";
134  }
135  else {
136  out_metric << "TS(";
137  out_grad << "VS(";
138  out_ls << "SS(";
139  out_hess << "SS(";
140  out_hessmat << "TS(";
141  }
142  for(std::size_t i = 0; i < e->getNumVertices(); i++) {
143  MVertex *ver = e->getVertex(i);
144  out_metric << ver->x() << "," << ver->y() << "," << ver->z();
145  out_grad << ver->x() << "," << ver->y() << "," << ver->z();
146  out_ls << ver->x() << "," << ver->y() << "," << ver->z();
147  out_hess << ver->x() << "," << ver->y() << "," << ver->z();
148  out_hessmat << ver->x() << "," << ver->y() << "," << ver->z();
149  if(i != e->getNumVertices() - 1) {
150  out_metric << ",";
151  out_grad << ",";
152  out_ls << ",";
153  out_hess << ",";
154  out_hessmat << ",";
155  }
156  else {
157  out_metric << "){";
158  out_grad << "){";
159  out_ls << "){";
160  out_hess << "){";
161  out_hessmat << "){";
162  }
163  }
164  for(std::size_t i = 0; i < e->getNumVertices(); i++) {
165  MVertex *ver = e->getVertex(i);
166  out_ls << vals[ver];
167  out_hess << (hessians[ver](0, 0) + hessians[ver](1, 1) +
168  hessians[ver](2, 2));
169  if(i == (e->getNumVertices() - 1)) {
170  out_ls << "};" << std::endl;
171  out_hess << "};" << std::endl;
172  }
173  else {
174  out_ls << ",";
175  out_hess << ",";
176  }
177  for(int k = 0; k < 3; k++) {
178  out_grad << grads[ver](k);
179  if((k == 2) && (i == (e->getNumVertices() - 1)))
180  out_grad << "};" << std::endl;
181  else
182  out_grad << ",";
183  for(int l = 0; l < 3; l++) {
184  out_metric << _nodalMetrics[ver](k, l);
185  out_hessmat << hessians[ver](k, l);
186  if((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) {
187  out_metric << "};" << std::endl;
188  out_hessmat << "};" << std::endl;
189  }
190  else {
191  out_metric << ",";
192  out_hessmat << ",";
193  }
194  }
195  }
196  }
197  }
198  out_grad << "};" << std::endl;
199  out_metric << "};" << std::endl;
200  out_ls << "};" << std::endl;
201  out_hess << "};" << std::endl;
202  out_hessmat << "};" << std::endl;
203  out_grad.close();
204  out_metric.close();
205  out_ls.close();
206  out_hess.close();
207  out_hessmat.close();
208 }
209 
210 meshMetric::~meshMetric()
211 {
212  if(_octree) delete _octree;
213  for(std::size_t i = 0; i < _elements.size(); i++) delete _elements[i];
214 }
215 
216 void meshMetric::computeValues()
217 {
218  auto it = _adj.begin();
219  while(it != _adj.end()) {
220  std::vector<MElement *> lt = it->second;
221  MVertex *ver = it->first;
222  vals[ver] = (*_fct)(ver->x(), ver->y(), ver->z());
223  it++;
224  }
225 }
226 
227 // Determines set of vertices to use for least squares
228 std::vector<MVertex *> getLSBlob(std::size_t minNbPt, v2t_cont::iterator it,
229  v2t_cont &adj)
230 {
231  // static const double RADFACT = 3;
232  //
233  // SPoint3 p0 = it->first->point(); // Vertex coordinates (center of circle)
234  //
235  // double rad = 0.;
236  // for (auto itEl = it->second.begin(); itEl != it->second.end(); itEl++)
237  // rad = std::max(rad,(*itEl)->getOuterRadius());
238  // rad *= RADFACT;
239 
240  std::vector<MVertex *> vv(1, it->first),
241  bvv = vv; // Vector of vertices in blob and in boundary of blob
242  do {
243  std::set<MVertex *> nbvv; // Set of vertices in new boundary
244  for(auto itBV = bvv.begin(); itBV != bvv.end();
245  itBV++) { // For each boundary vertex...
246  std::vector<MElement *> &adjBV = adj[*itBV];
247  for(auto itBVEl = adjBV.begin(); itBVEl != adjBV.end(); itBVEl++) {
248  for(std::size_t iV = 0; iV < (*itBVEl)->getNumVertices();
249  iV++) { // ... look for adjacent vertices...
250  MVertex *v = (*itBVEl)->getVertex(iV);
251  // if ((find(vv.begin(),vv.end(),v) == vv.end()) &&
252  // (p0.distance(v->point()) <= rad)) nbvv.insert(v);
253  if(find(vv.begin(), vv.end(), v) == vv.end())
254  nbvv.insert(v); // ... and add them in the new boundary if they are
255  // not already in the blob
256  }
257  }
258  }
259  if(nbvv.empty())
260  bvv.clear();
261  else {
262  bvv.assign(nbvv.begin(), nbvv.end());
263  vv.insert(vv.end(), nbvv.begin(), nbvv.end());
264  }
265  // } while (!bvv.empty());
266  } while(vv.size() < minNbPt); // Repeat until min. number of points is reached
267 
268  return vv;
269 }
270 
271 // Compute derivatives and second order derivatives using least squares
272 // 2D LS system: a_i0*x^2+a_i1*x*y+a_i2*y^2+a_i3*x+a_i4*y+a_i5=b_i
273 // 3D LS system:
274 // a_i0*x^2+a_i1*x*y+a_i2*x*z+a_i3*y^2+a_i4*y*z+a_i5*z^2+a_i6*x+a_i7*y+a_i8*z+a_i9=b_i
275 void meshMetric::computeHessian()
276 {
277  std::size_t sysDim = (_dim == 2) ? 6 : 10;
278  std::size_t minNbPtBlob = 3 * sysDim;
279 
280  for(auto it = _adj.begin(); it != _adj.end(); it++) {
281  MVertex *ver = it->first;
282  std::vector<MVertex *> vv = getLSBlob(minNbPtBlob, it, _adj);
283  fullMatrix<double> A(vv.size(), sysDim), ATA(sysDim, sysDim);
284  fullVector<double> b(vv.size()), ATb(sysDim), coeffs(sysDim);
285  for(std::size_t i = 0; i < vv.size(); i++) {
286  const double &x = vv[i]->x(), &y = vv[i]->y(), &z = vv[i]->z();
287  if(_dim == 2) {
288  A(i, 0) = x * x;
289  A(i, 1) = x * y;
290  A(i, 2) = y * y;
291  A(i, 3) = x;
292  A(i, 4) = y;
293  A(i, 5) = 1.;
294  }
295  else {
296  A(i, 0) = x * x;
297  A(i, 1) = x * y;
298  A(i, 2) = x * z;
299  A(i, 3) = y * y;
300  A(i, 4) = y * z;
301  A(i, 5) = z * z;
302  A(i, 6) = x;
303  A(i, 7) = y;
304  A(i, 8) = z;
305  A(i, 9) = 1.;
306  }
307  b(i) = vals[vv[i]];
308  }
309  ATA.gemm(A, A, 1., 0., true, false);
310  A.multWithATranspose(b, 1., 0., ATb);
311  ATA.luSolve(ATb, coeffs);
312  const double &x = ver->x(), &y = ver->y(), &z = ver->z();
313  double d2udx2, d2udy2, d2udz2, d2udxy, d2udxz, d2udyz, dudx, dudy, dudz;
314  if(_dim == 2) {
315  d2udx2 = 2. * coeffs(0);
316  d2udy2 = 2. * coeffs(2);
317  d2udz2 = 0.;
318  d2udxy = coeffs(1);
319  d2udxz = 0.;
320  d2udyz = 0.;
321  dudx = d2udx2 * x + d2udxy * y + coeffs(3);
322  dudy = d2udxy * x + d2udy2 * y + coeffs(4);
323  dudz = 0.;
324  }
325  else {
326  d2udx2 = 2. * coeffs(0);
327  d2udy2 = 2. * coeffs(3);
328  d2udz2 = 2. * coeffs(5);
329  d2udxy = coeffs(1);
330  d2udxz = coeffs(2);
331  d2udyz = coeffs(4);
332  dudx = d2udx2 * x + d2udxy * y + d2udxz * z + coeffs(6);
333  dudy = d2udxy * x + d2udy2 * y + d2udyz * z + coeffs(7);
334  dudz = d2udxz * x + d2udyz * y + d2udz2 * z + coeffs(8);
335  }
336  double duNorm = sqrt(dudx * dudx + dudy * dudy + dudz * dudz);
337  if(duNorm == 0. || _technique == meshMetric::HESSIAN ||
338  _technique == meshMetric::EIGENDIRECTIONS ||
339  _technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H)
340  duNorm = 1.;
341  grads[ver] = SVector3(dudx / duNorm, dudy / duNorm, dudz / duNorm);
342  hessians[ver](0, 0) = d2udx2;
343  hessians[ver](0, 1) = d2udxy;
344  hessians[ver](0, 2) = d2udxz;
345  hessians[ver](1, 0) = d2udxy;
346  hessians[ver](1, 1) = d2udy2;
347  hessians[ver](1, 2) = d2udyz;
348  hessians[ver](2, 0) = d2udxz;
349  hessians[ver](2, 1) = d2udyz;
350  hessians[ver](2, 2) = d2udz2;
351  }
352 }
353 
354 void meshMetric::computeMetricLevelSet(MVertex *ver, SMetric3 &hessian,
355  SMetric3 &metric, double &size, double x,
356  double y, double z)
357 {
358  double signed_dist;
359  SVector3 gr;
360  if(ver) {
361  signed_dist = vals[ver];
362  gr = grads[ver];
363  hessian = hessians[ver];
364  }
365  else {
366  signed_dist = (*_fct)(x, y, z);
367  _fct->gradient(x, y, z, gr(0), gr(1), gr(2));
368  _fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
369  hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
370  hessian(2, 1), hessian(2, 2));
371  }
372 
373  double dist = fabs(signed_dist);
374 
375  SMetric3 H;
376  double norm = gr(0) * gr(0) + gr(1) * gr(1) + gr(2) * gr(2);
377  if(dist < _e && norm != 0.0) {
378  double h = hmin * (hmax / hmin - 1) * dist / _e + hmin;
379  double C = 1. / (h * h) - 1. / (hmax * hmax);
380  H(0, 0) += C * gr(0) * gr(0) / norm;
381  H(1, 1) += C * gr(1) * gr(1) / norm;
382  H(2, 2) += C * gr(2) * gr(2) / norm;
383  H(1, 0) = H(0, 1) = C * gr(1) * gr(0) / norm;
384  H(2, 0) = H(0, 2) = C * gr(2) * gr(0) / norm;
385  H(2, 1) = H(1, 2) = C * gr(2) * gr(1) / norm;
386  }
387 
388  fullMatrix<double> V(3, 3);
389  fullVector<double> S(3);
390  H.eig(V, S);
391 
392  double lambda1, lambda2, lambda3;
393  lambda1 = S(0);
394  lambda2 = S(1);
395  lambda3 = (_dim == 3) ? S(2) : 1.;
396 
397  SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
398  SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
399  SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
400 
401  size =
402  std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
403  metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
404 }
405 
406 void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian,
407  SMetric3 &metric, double &size, double x,
408  double y, double z)
409 {
410  SVector3 gr;
411  if(ver != nullptr) {
412  gr = grads[ver];
413  hessian = hessians[ver];
414  }
415  else if(ver == nullptr) {
416  _fct->gradient(x, y, z, gr(0), gr(1), gr(2));
417  _fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
418  hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
419  hessian(2, 1), hessian(2, 2));
420  }
421 
422  double _epsilonP = 1.;
423  double hminP = 1.e-12;
424  double hmaxP = 1.e+12;
425 
426  fullMatrix<double> V(3, 3);
427  fullVector<double> S(3);
428  hessian.eig(V, S);
429 
430  double lambda1 =
431  std::min(std::max(fabs(S(0)) / _epsilonP, 1. / (hmaxP * hmaxP)),
432  1. / (hminP * hminP));
433  double lambda2 =
434  std::min(std::max(fabs(S(1)) / _epsilonP, 1. / (hmaxP * hmaxP)),
435  1. / (hminP * hminP));
436  double lambda3 =
437  (_dim == 3) ?
438  std::min(std::max(fabs(S(2)) / _epsilonP, 1. / (hmaxP * hmaxP)),
439  1. / (hminP * hminP)) :
440  1.;
441 
442  SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
443  SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
444  SVector3 t3 =
445  (_dim == 3) ? SVector3(V(0, 2), V(1, 2), V(2, 2)) : SVector3(0., 0., 1.);
446 
447  size =
448  std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
449  metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
450 }
451 
452 void meshMetric::computeMetricFrey(MVertex *ver, SMetric3 &hessian,
453  SMetric3 &metric, double &size, double x,
454  double y, double z)
455 {
456  double signed_dist;
457  SVector3 gr;
458  if(ver) {
459  signed_dist = vals[ver];
460  gr = grads[ver];
461  hessian = hessians[ver];
462  }
463  else {
464  signed_dist = (*_fct)(x, y, z);
465  _fct->gradient(x, y, z, gr(0), gr(1), gr(2));
466  _fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
467  hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
468  hessian(2, 1), hessian(2, 2));
469  }
470 
471  double dist = fabs(signed_dist);
472 
473  SMetric3 H(1. / (hmax * hmax));
474  double norm = gr(0) * gr(0) + gr(1) * gr(1) + gr(2) * gr(2);
475  if(dist < _e && norm != 0.0) {
476  double h = hmin * (hmax / hmin - 1.0) * dist / _e + hmin;
477  double C = 1. / (h * h) - 1. / (hmax * hmax);
478  double kappa = hessian(0, 0) + hessian(1, 1) + hessian(2, 2);
479  double epsGeom = 4.0 * 3.14 * 3.14 / (kappa * _np * _np);
480  H(0, 0) += C * gr(0) * gr(0) / (norm) + hessian(0, 0) / epsGeom;
481  H(1, 1) += C * gr(1) * gr(1) / (norm) + hessian(1, 1) / epsGeom;
482  H(2, 2) += C * gr(2) * gr(2) / (norm) + hessian(2, 2) / epsGeom;
483  H(1, 0) = H(0, 1) = C * gr(1) * gr(0) / (norm) + hessian(1, 0) / epsGeom;
484  H(2, 0) = H(0, 2) = C * gr(2) * gr(0) / (norm) + hessian(2, 0) / epsGeom;
485  H(2, 1) = H(1, 2) = C * gr(2) * gr(1) / (norm) + hessian(2, 1) / epsGeom;
486  }
487 
488  fullMatrix<double> V(3, 3);
489  fullVector<double> S(3);
490  H.eig(V, S);
491 
492  double lambda1, lambda2, lambda3;
493  lambda1 = S(0);
494  lambda2 = S(1);
495  lambda3 = (_dim == 3) ? S(2) : 1.;
496 
497  if(dist < _e) {
498  lambda1 = std::min(std::max(fabs(S(0)) / _epsilon, 1. / (hmax * hmax)),
499  1. / (hmin * hmin));
500  lambda2 = std::min(std::max(fabs(S(1)) / _epsilon, 1. / (hmax * hmax)),
501  1. / (hmin * hmin));
502  lambda3 = (_dim == 3) ?
503  std::min(std::max(fabs(S(2)) / _epsilon, 1. / (hmax * hmax)),
504  1. / (hmin * hmin)) :
505  1.;
506  }
507 
508  SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
509  SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
510  SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
511 
512  size =
513  std::min(std::min(1 / sqrt(lambda1), 1 / sqrt(lambda2)), 1 / sqrt(lambda3));
514  metric = SMetric3(lambda1, lambda2, lambda3, t1, t2, t3);
515 }
516 
517 void meshMetric::computeMetricEigenDir(MVertex *ver, SMetric3 &hessian,
518  SMetric3 &metric, double &size, double x,
519  double y, double z)
520 {
521  double signed_dist;
522  SVector3 gVec;
523  if(ver) {
524  signed_dist = vals[ver];
525  gVec = grads[ver];
526  hessian = hessians[ver];
527  }
528  else {
529  signed_dist = (*_fct)(x, y, z);
530  _fct->gradient(x, y, z, gVec(0), gVec(1), gVec(2));
531  _fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
532  hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
533  hessian(2, 1), hessian(2, 2));
534  }
535 
536  double dist = fabs(signed_dist);
537 
538  const double metric_value_hmax = 1. / (hmax * hmax);
539  const double gMag = gVec.norm(), invGMag = 1. / gMag;
540 
541  if(signed_dist < _e && signed_dist > _e_moins && gMag != 0.0) {
542  const double metric_value_hmin = 1. / (hmin * hmin);
543  const SVector3 nVec = invGMag * gVec; // Unit normal vector
544  double lambda_n =
545  0.; // Eigenvalues of metric for normal & tangential directions
546  if(_technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H) {
547  // const double h_dist = hmin + ((hmax-hmin)/_e)*dist;
548  // // Characteristic element size in the normal direction - linear interp
549  // between hmin and hmax lambda_n = 1./(h_dist*h_dist);
550  double beta = CTX::instance()->mesh.smoothRatio;
551  double h_dista = std::min((hmin + (dist * log(beta))), hmax);
552  lambda_n = 1. / (h_dista * h_dista);
553  }
554  else if(_technique == meshMetric::EIGENDIRECTIONS) {
555  const double maximum_distance =
556  (signed_dist > 0.) ? _e : fabs(_e_moins); // ... or linear interpolation
557  // between 1/h_min^2 and
558  // 1/h_max^2
559  lambda_n =
560  metric_value_hmin +
561  ((metric_value_hmax - metric_value_hmin) / maximum_distance) * dist;
562  }
563  std::vector<SVector3> tVec; // Unit tangential vectors
564  std::vector<double> kappa; // Curvatures
565  if(_dim == 2) { // 2D curvature formula: cf. R. Goldman, "Curvature formulas
566  // for implicit curves and surfaces", Computer Aided
567  // Geometric Design 22 (2005), pp. 632–658
568  kappa.resize(2);
569  kappa[0] =
570  fabs(-gVec(1) * (-gVec(1) * hessian(0, 0) + gVec(0) * hessian(0, 1)) +
571  gVec(0) * (-gVec(1) * hessian(1, 0) + gVec(0) * hessian(1, 1))) *
572  pow(invGMag, 3);
573  kappa[1] = 1.;
574  tVec.resize(2);
575  tVec[0] = SVector3(-nVec(1), nVec(0), 0.);
576  tVec[1] = SVector3(0., 0., 1.);
577  }
578  else { // 3D curvature formula: cf. A.G. Belyaev, A.A. Pasko and T.L. Kunii,
579  // "Ridges and Ravines on Implicit Surfaces," CGI, pp.530-535,
580  // Computer Graphics International 1998 (CGI'98), 1998
581  fullMatrix<double> ImGG(3, 3);
582  ImGG(0, 0) = 1. - gVec(0) * gVec(0);
583  ImGG(0, 1) = -gVec(0) * gVec(1);
584  ImGG(0, 2) = -gVec(0) * gVec(2);
585  ImGG(1, 0) = -gVec(1) * gVec(0);
586  ImGG(1, 1) = 1. - gVec(1) * gVec(1);
587  ImGG(1, 2) = -gVec(1) * gVec(2);
588  ImGG(2, 0) = -gVec(2) * gVec(0);
589  ImGG(2, 1) = -gVec(2) * gVec(1);
590  ImGG(2, 2) = 1. - gVec(2) * gVec(2);
591  fullMatrix<double> hess(3, 3);
592  hessian.getMat(hess);
593  fullMatrix<double> gN(3, 3); // Gradient of unit normal
594  gN.gemm(ImGG, hess, 1., 0.);
595  gN.scale(invGMag);
596  fullMatrix<double> eigVecL(3, 3), eigVecR(3, 3);
597  fullVector<double> eigValRe(3), eigValIm(3);
598  gN.eig(eigValRe, eigValIm, eigVecL, eigVecR,
599  false); // Eigendecomp. of gradient of unit normal
600  kappa.resize(3); // Store abs. val. of eigenvalues (= principal curvatures
601  // only in non-normal directions)
602  kappa[0] = fabs(eigValRe(0));
603  kappa[1] = fabs(eigValRe(1));
604  kappa[2] = fabs(eigValRe(2));
605  tVec.resize(3); // Store normalized eigenvectors (= principal directions
606  // only in non-normal directions)
607  tVec[0] = SVector3(eigVecR(0, 0), eigVecR(1, 0), eigVecR(2, 0));
608  tVec[0].normalize();
609  tVec[1] = SVector3(eigVecR(0, 1), eigVecR(1, 1), eigVecR(2, 1));
610  tVec[1].normalize();
611  tVec[2] = SVector3(eigVecR(0, 2), eigVecR(1, 2), eigVecR(2, 2));
612  tVec[2].normalize();
613  std::vector<double> tVecDotNVec(3); // Store dot products with normal
614  // vector to look for normal direction
615  tVecDotNVec[0] = fabs(dot(tVec[0], nVec));
616  tVecDotNVec[1] = fabs(dot(tVec[1], nVec));
617  tVecDotNVec[2] = fabs(dot(tVec[2], nVec));
618  const int i_N = max_element(tVecDotNVec.begin(), tVecDotNVec.end()) -
619  tVecDotNVec.begin(); // Index of normal dir. = max. dot
620  // products (close to 0. in
621  // tangential dir.)
622  kappa.erase(kappa.begin() + i_N); // Remove normal dir.
623  tVec.erase(tVec.begin() + i_N);
624  }
625  const double invh_t0 = (_np * kappa[0]) / 6.283185307,
626  invh_t1 = (_np * kappa[1]) /
627  6.283185307; // Inverse of tangential size 0
628  const double lambda_t0 = invh_t0 * invh_t0, lambda_t1 = invh_t1 * invh_t1;
629  const double lambdaP_n =
630  std::min(std::max(lambda_n, metric_value_hmax),
631  metric_value_hmin); // Truncate eigenvalues
632  const double lambdaP_t0 =
633  std::min(std::max(lambda_t0, metric_value_hmax), metric_value_hmin);
634  const double lambdaP_t1 =
635  (_dim == 2) ?
636  1. :
637  std::min(std::max(lambda_t1, metric_value_hmax), metric_value_hmin);
638  metric =
639  SMetric3(lambdaP_n, lambdaP_t0, lambdaP_t1, nVec, tVec[0], tVec[1]);
640  const double h_n = 1. / sqrt(lambdaP_n), h_t0 = 1. / sqrt(lambdaP_t0),
641  h_t1 = 1. / sqrt(lambdaP_t1);
642  size = std::min(std::min(h_n, h_t0), h_t1);
643  }
644  else { // isotropic metric !
645  SMetric3 mymetric(metric_value_hmax);
646  metric = mymetric;
647  size = hmax;
648  }
649 }
650 
651 void meshMetric::computeMetricIsoLinInterp(MVertex *ver, SMetric3 &hessian,
652  SMetric3 &metric, double &size,
653  double x, double y, double z)
654 {
655  double signed_dist;
656  SVector3 gr;
657  if(ver) {
658  signed_dist = vals[ver];
659  gr = grads[ver];
660  hessian = hessians[ver];
661  }
662  else {
663  signed_dist = (*_fct)(x, y, z);
664  _fct->gradient(x, y, z, gr(0), gr(1), gr(2));
665  _fct->hessian(x, y, z, hessian(0, 0), hessian(0, 1), hessian(0, 2),
666  hessian(1, 0), hessian(1, 1), hessian(1, 2), hessian(2, 0),
667  hessian(2, 1), hessian(2, 2));
668  }
669 
670  double dist = fabs(signed_dist);
671  double norm = gr.normalize();
672  size = hmax; // the characteristic element size in all directions - linear
673  // interp between hmin and hmax
674  if(norm != 0.) {
675  if((signed_dist >= 0) && (signed_dist < _e))
676  size = hmin + ((hmax - hmin) / _e) * dist;
677  else if((signed_dist < 0) && (signed_dist > _e_moins))
678  size = hmin - ((hmax - hmin) / _e_moins) * dist;
679  }
680 
681  double lambda = 1. / size / size;
682  metric(0, 0) = lambda;
683  metric(0, 1) = 0.;
684  metric(0, 2) = 0.;
685  metric(1, 0) = 0.;
686  metric(1, 1) = lambda;
687  metric(1, 2) = 0.;
688  metric(2, 0) = 0.;
689  metric(2, 1) = 0.;
690  metric(2, 2) = (_dim == 3) ? lambda : 1.;
691 }
692 
693 // this function scales the mesh metric in order
694 // to reach a target number of elements
695 // We know that the number of elements in the final
696 // mesh will be (assuming M_e the metric at centroid of element e)
697 // N = \sum_e \sqrt {\det (M_e)} V_e
698 // where V_e is the volume of e
699 // assuming that N_{target} = K N, we have
700 // K N = K \sum_e \sqrt {\det (M_e)} V_e
701 // = \sum_e \sqrt {K^2 \det (M_e)} V_e
702 // = \sum_e \sqrt {\det (K^{2/d} M_e)} V_e
703 // where d is the dimension of the problem.
704 // This means that the metric should be scaled by K^{2/d} where
705 // K is N_target / N
706 void meshMetric::scaleMetric(int nbElementsTarget, nodalMetricTensor &nmt)
707 {
708  // compute N
709  double N = 0;
710  for(std::size_t i = 0; i < _elements.size(); i++) {
711  MElement *e = _elements[i];
712  SMetric3 m1 = nmt[e->getVertex(0)];
713  SMetric3 m2 = nmt[e->getVertex(1)];
714  SMetric3 m3 = nmt[e->getVertex(2)];
715  if(_dim == 2) {
716  SMetric3 m = interpolation(m1, m2, m3, 0.3333, 0.3333);
717  N += sqrt(m.determinant()) * e->getVolume() * 4. / sqrt(3.0); // 3.0
718  }
719  else {
720  SMetric3 m4 = nmt[e->getVertex(3)];
721  SMetric3 m = interpolation(m1, m2, m3, m4, 0.25, 0.25, 0.25);
722  N += sqrt(m.determinant()) * e->getVolume() * 12. / sqrt(2.0); // 4.0;
723  }
724  }
725  double scale = pow((double)nbElementsTarget / N, 2.0 / _dim);
726  for(auto it = nmt.begin(); it != nmt.end(); ++it) {
727  if(_dim == 3) { it->second *= scale; }
728  else {
729  it->second(0, 0) *= scale;
730  it->second(1, 0) *= scale;
731  it->second(1, 1) *= scale;
732  }
733  SMetric3 &m = it->second;
734  fullMatrix<double> V(3, 3);
735  fullVector<double> S(3);
736  m.eig(V, S);
737  S(0) = std::min(std::max(S(0), 1 / (hmax * hmax)), 1 / (hmin * hmin));
738  S(1) = std::min(std::max(S(1), 1 / (hmax * hmax)), 1 / (hmin * hmin));
739  if(_dim == 3)
740  S(2) = std::min(std::max(S(2), 1 / (hmax * hmax)), 1 / (hmin * hmin));
741  SVector3 t1(V(0, 0), V(1, 0), V(2, 0));
742  SVector3 t2(V(0, 1), V(1, 1), V(2, 1));
743  SVector3 t3(V(0, 2), V(1, 2), V(2, 2));
744  m = SMetric3(S(0), S(1), S(2), t1, t2, t3);
745  }
746 }
747 
748 void meshMetric::computeMetric(int metricNumber)
749 {
750  _fct = setOfFcts[metricNumber];
751  std::vector<double> parameters = setOfParameters[metricNumber];
752  int technique = setOfTechniques[metricNumber];
753 
754  hmin = parameters.size() >= 3 ? parameters[1] : CTX::instance()->mesh.lcMin;
755  hmax = parameters.size() >= 3 ? parameters[2] : CTX::instance()->mesh.lcMax;
756  _e = parameters[0];
757  _e_moins = (parameters.size() >= 5) ? parameters[4] : -parameters[0];
758  if(_e_moins > 0.) _e_moins *= -1.;
759  _epsilon = parameters[0];
760  _technique = (MetricComputationTechnique)technique;
761  _np = (parameters.size() >= 4) ? parameters[3] : 15.;
762 
763  computeValues();
764  computeHessian();
765 
766  for(auto it = _adj.begin(); it != _adj.end(); it++) {
767  MVertex *ver = it->first;
768  SMetric3 hessian, metric;
769  double size;
770  switch(_technique) {
771  case(LEVELSET): computeMetricLevelSet(ver, hessian, metric, size); break;
772  case(HESSIAN): computeMetricHessian(ver, hessian, metric, size); break;
773  case(FREY): computeMetricFrey(ver, hessian, metric, size); break;
774  case(EIGENDIRECTIONS):
775  computeMetricEigenDir(ver, hessian, metric, size);
776  break;
777  case(EIGENDIRECTIONS_LINEARINTERP_H):
778  computeMetricEigenDir(ver, hessian, metric, size);
779  break;
780  case(ISOTROPIC_LINEARINTERP_H):
781  computeMetricIsoLinInterp(ver, hessian, metric, size);
782  break;
783  }
784 
785  setOfSizes[metricNumber].insert(std::make_pair(ver, size));
786  setOfMetrics[metricNumber].insert(std::make_pair(ver, metric));
787  }
788 
789  if(_technique == HESSIAN) scaleMetric(_epsilon, setOfMetrics[metricNumber]);
790 }
791 
792 double meshMetric::operator()(double x, double y, double z, GEntity *ge)
793 {
794  if(needMetricUpdate) updateMetrics();
795  if(!setOfMetrics.size()) {
796  Msg::Error("No metric defined");
797  return 0.;
798  }
799  SPoint3 xyz(x, y, z), uvw;
800  double initialTol = CTX::instance()->mesh.toleranceReferenceElement;
801  CTX::instance()->mesh.toleranceReferenceElement = 1.e-4;
802  MElement *e = _octree->find(x, y, z, _dim);
803  CTX::instance()->mesh.toleranceReferenceElement = initialTol;
804  double value = 0.;
805  if(e) {
806  e->xyz2uvw(xyz, uvw);
807  double *val = new double[e->getNumVertices()];
808  for(std::size_t i = 0; i < e->getNumVertices(); i++) {
809  val[i] = _nodalSizes[e->getVertex(i)];
810  }
811  value = e->interpolate(val, uvw[0], uvw[1], uvw[2]);
812  delete[] val;
813  }
814  else {
815  Msg::Warning("point %g %g %g not found, looking for nearest node", x, y, z);
816  double minDist = 1.e100;
817  for(auto it = _nodalSizes.begin(); it != _nodalSizes.end(); it++) {
818  const double dist = xyz.distance(it->first->point());
819  if(dist <= minDist) {
820  minDist = dist;
821  value = it->second;
822  }
823  }
824  }
825  return value;
826 }
827 
828 void meshMetric::operator()(double x, double y, double z, SMetric3 &metr,
829  GEntity *ge)
830 {
831  if(needMetricUpdate) { updateMetrics(); }
832  if(!setOfMetrics.size()) {
833  Msg::Error("No metric defined");
834  return;
835  }
836  metr = SMetric3(1.e-22);
837 
838  // RECOMPUTE MESH METRIC AT XYZ
839  if(hasAnalyticalMetric) {
840  int nbMetrics = setOfMetrics.size();
841  std::vector<SMetric3> newSetOfMetrics(nbMetrics);
842  for(int iMetric = 0; iMetric < nbMetrics; iMetric++) {
843  _fct = setOfFcts[iMetric];
844  _technique = (MetricComputationTechnique)setOfTechniques[iMetric];
845  if(_fct->hasDerivatives()) {
846  SMetric3 hessian, metric;
847  double size;
848  switch(_technique) {
849  case(LEVELSET):
850  computeMetricLevelSet(nullptr, hessian, metric, size, x, y, z);
851  break;
852  case(HESSIAN):
853  computeMetricHessian(nullptr, hessian, metric, size, x, y, z);
854  break;
855  case(FREY):
856  computeMetricFrey(nullptr, hessian, metric, size, x, y, z);
857  break;
858  case(EIGENDIRECTIONS):
859  computeMetricEigenDir(nullptr, hessian, metric, size, x, y, z);
860  break;
861  case(EIGENDIRECTIONS_LINEARINTERP_H):
862  computeMetricEigenDir(nullptr, hessian, metric, size, x, y, z);
863  break;
864  case(ISOTROPIC_LINEARINTERP_H):
865  computeMetricIsoLinInterp(nullptr, hessian, metric, size, x, y, z);
866  break;
867  }
868  newSetOfMetrics[iMetric] = metric;
869  }
870  else {
871  // find other metrics here
872  SMetric3 metric;
873  SPoint3 xyz(x, y, z), uvw;
874  double initialTol = CTX::instance()->mesh.toleranceReferenceElement;
875  CTX::instance()->mesh.toleranceReferenceElement = 1.e-4;
876  MElement *e = _octree->find(x, y, z, _dim);
877  CTX::instance()->mesh.toleranceReferenceElement = initialTol;
878  if(e) {
879  e->xyz2uvw(xyz, uvw);
880  SMetric3 m1 = setOfMetrics[iMetric][e->getVertex(0)];
881  SMetric3 m2 = setOfMetrics[iMetric][e->getVertex(1)];
882  SMetric3 m3 = setOfMetrics[iMetric][e->getVertex(2)];
883  if(_dim == 2)
884  metric = interpolation(m1, m2, m3, uvw[0], uvw[1]);
885  else {
886  SMetric3 m4 = setOfMetrics[iMetric][e->getVertex(3)];
887  metric = interpolation(m1, m2, m3, m4, uvw[0], uvw[1], uvw[2]);
888  }
889  newSetOfMetrics[iMetric] = metric;
890  }
891  else {
892  Msg::Warning("point %g %g %g not found, looking for nearest node", x,
893  y, z);
894  }
895  }
896  }
897  // intersect metrics here
898  metr = newSetOfMetrics[0];
899  for(int i = 1; i < nbMetrics; i++)
900  metr = intersection_conserve_mostaniso(metr, newSetOfMetrics[i]);
901  }
902  // INTERPOLATE DISCRETE MESH METRIC
903  else {
904  SPoint3 xyz(x, y, z), uvw;
905  double initialTol = CTX::instance()->mesh.toleranceReferenceElement;
906  CTX::instance()->mesh.toleranceReferenceElement = 1.e-4;
907  MElement *e = _octree->find(x, y, z, _dim);
908  CTX::instance()->mesh.toleranceReferenceElement = initialTol;
909 
910  if(e) {
911  e->xyz2uvw(xyz, uvw);
912  SMetric3 m1 = _nodalMetrics[e->getVertex(0)];
913  SMetric3 m2 = _nodalMetrics[e->getVertex(1)];
914  SMetric3 m3 = _nodalMetrics[e->getVertex(2)];
915  if(_dim == 2)
916  metr = interpolation(m1, m2, m3, uvw[0], uvw[1]);
917  else {
918  SMetric3 m4 = _nodalMetrics[e->getVertex(3)];
919  metr = interpolation(m1, m2, m3, m4, uvw[0], uvw[1], uvw[2]);
920  }
921  }
922  else {
923  Msg::Warning("point %g %g %g not found, looking for nearest node", x, y,
924  z);
925  double minDist = 1.e100;
926  for(auto it = _nodalMetrics.begin(); it != _nodalMetrics.end(); it++) {
927  const double dist = xyz.distance(it->first->point());
928  if(dist <= minDist) {
929  minDist = dist;
930  metr = it->second;
931  }
932  }
933  }
934  }
935 }
936 
937 double meshMetric::getLaplacian(MVertex *v)
938 {
939  MVertex *vNew = _vertexMap[v->getNum()];
940  auto it = hessians.find(vNew);
941  SMetric3 h = it->second;
942  return h(0, 0) + h(1, 1) + h(2, 2);
943 }
944 
945 SVector3 meshMetric::getGradient(MVertex *v)
946 {
947  MVertex *vNew = _vertexMap[v->getNum()];
948  auto it = grads.find(vNew);
949  SVector3 gr = it->second;
950  return gr;
951 }
buildVertexToElement
void buildVertexToElement(std::vector< T * > const &elements, v2t_cont &adj)
Definition: meshGFaceOptimize.h:38
meshMetric::setOfFcts
std::map< int, simpleFunction< double > * > setOfFcts
Definition: meshMetric.h:65
MElement::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MElement.cpp:1128
contextMeshOptions::toleranceReferenceElement
double toleranceReferenceElement
Definition: Context.h:47
meshMetric::_octree
MElementOctree * _octree
Definition: meshMetric.h:47
simpleFunction::hessian
virtual void hessian(double x, double y, double z, scalar &dfdxx, scalar &dfdxy, scalar &dfdxz, scalar &dfdyx, scalar &dfdyy, scalar &dfdyz, scalar &dfdzx, scalar &dfdzy, scalar &dfdzz) const
Definition: simpleFunction.h:30
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
meshMetric::_dim
int _dim
Definition: meshMetric.h:37
SMetric3
Definition: STensor3.h:17
meshMetric::_vertexMap
std::map< int, MVertex * > _vertexMap
Definition: meshMetric.h:48
fullVector< double >
meshMetric::ISOTROPIC_LINEARINTERP_H
@ ISOTROPIC_LINEARINTERP_H
Definition: meshMetric.h:30
meshMetric::EIGENDIRECTIONS_LINEARINTERP_H
@ EIGENDIRECTIONS_LINEARINTERP_H
Definition: meshMetric.h:29
meshMetric::exportInfo
void exportInfo(const char *fileendname)
Definition: meshMetric.cpp:105
meshMetric::setOfSizes
std::map< int, nodalField > setOfSizes
Definition: meshMetric.h:63
OS.h
GEntity.h
meshMetric::needMetricUpdate
bool needMetricUpdate
Definition: meshMetric.h:39
meshMetric::_technique
meshMetric::MetricComputationTechnique _technique
Definition: meshMetric.h:41
MVertex
Definition: MVertex.h:24
MElement::getDim
virtual int getDim() const =0
simpleFunction::gradient
virtual void gradient(double x, double y, double z, scalar &dfdx, scalar &dfdy, scalar &dfdz) const
Definition: simpleFunction.h:25
MElementOctree.h
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
meshMetric::setOfMetrics
std::map< int, nodalMetricTensor > setOfMetrics
Definition: meshMetric.h:62
SMetric3::eig
void eig(fullMatrix< double > &V, fullVector< double > &S, bool s=false) const
Definition: STensor3.cpp:104
SPoint3
Definition: SPoint3.h:14
meshMetric::computeMetricEigenDir
void computeMetricEigenDir(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x=0.0, double y=0.0, double z=0.0)
Definition: meshMetric.cpp:517
MVertex::getNum
std::size_t getNum() const
Definition: MVertex.h:86
getLSBlob
std::vector< MVertex * > getLSBlob(std::size_t minNbPt, v2t_cont::iterator it, v2t_cont &adj)
Definition: meshMetric.cpp:228
fullMatrix::scale
void scale(const scalar s)
Definition: fullMatrix.h:447
contextMeshOptions::lcMin
double lcMin
Definition: Context.h:25
SVector3
Definition: SVector3.h:16
meshGFaceOptimize.h
meshMetric::updateMetrics
void updateMetrics()
Definition: meshMetric.cpp:79
meshMetric::_epsilon
double _epsilon
Definition: meshMetric.h:38
meshMetric::setOfParameters
std::map< int, std::vector< double > > setOfParameters
Definition: meshMetric.h:66
v2t_cont
std::map< MVertex *, std::vector< MElement * >, MVertexPtrLessThan > v2t_cont
Definition: meshGFaceOptimize.h:31
meshMetric::_nodalMetrics
nodalMetricTensor _nodalMetrics
Definition: meshMetric.h:59
MElement::copy
virtual MElement * copy(std::map< int, MVertex * > &vertexMap, std::map< MElement *, MElement * > &newParents, std::map< MElement *, MElement * > &newDomains)
Definition: MElement.cpp:2486
GEntity
Definition: GEntity.h:31
MElement::getNumVertices
virtual std::size_t getNumVertices() const =0
meshMetric.h
SMetric3::determinant
double determinant() const
Definition: STensor3.cpp:74
meshMetric::computeMetricLevelSet
void computeMetricLevelSet(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x=0.0, double y=0.0, double z=0.0)
Definition: meshMetric.cpp:354
contextMeshOptions::lcMax
double lcMax
Definition: Context.h:25
meshMetric::_adj
v2t_cont _adj
Definition: meshMetric.h:46
meshMetric::_fct
simpleFunction< double > * _fct
Definition: meshMetric.h:43
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
fullMatrix< double >
GModel::lastFace
fiter lastFace()
Definition: GModel.h:359
contextMeshOptions::smoothRatio
double smoothRatio
Definition: Context.h:26
meshMetric::MetricComputationTechnique
MetricComputationTechnique
Definition: meshMetric.h:24
meshMetric::computeMetricIsoLinInterp
void computeMetricIsoLinInterp(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x=0.0, double y=0.0, double z=0.0)
Definition: meshMetric.cpp:651
meshMetric::HESSIAN
@ HESSIAN
Definition: meshMetric.h:26
SMetric3::getMat
void getMat(fullMatrix< double > &mat) const
Definition: STensor3.cpp:51
simpleFunction< double >
meshMetric::computeMetricFrey
void computeMetricFrey(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x=0.0, double y=0.0, double z=0.0)
Definition: meshMetric.cpp:452
meshMetric::scaleMetric
void scaleMetric(int nbElementsTarget, nodalMetricTensor &nmt)
Definition: meshMetric.cpp:706
meshMetric::hasAnalyticalMetric
bool hasAnalyticalMetric
Definition: meshMetric.h:40
meshMetric::_e
double _e
Definition: meshMetric.h:38
meshMetric::_np
double _np
Definition: meshMetric.h:38
MElement::getVolume
virtual double getVolume()
Definition: MElement.cpp:567
GModel::getDim
int getDim() const
Definition: GModel.cpp:989
norm
void norm(const double *vec, double *norm)
Definition: gmshLevelset.cpp:202
meshMetric::_elements
std::vector< MElement * > _elements
Definition: meshMetric.h:45
GModel
Definition: GModel.h:44
meshMetric::meshMetric
meshMetric(std::vector< MElement * > elements)
Definition: meshMetric.cpp:46
meshMetric::nodalMetricTensor
std::map< MVertex *, SMetric3 > nodalMetricTensor
Definition: meshMetric.h:55
GModel::firstFace
fiter firstFace()
Definition: GModel.h:355
MElementOctree
Definition: MElementOctree.h:15
SVector3::norm
double norm() const
Definition: SVector3.h:33
interpolation
SMetric3 interpolation(const SMetric3 &m1, const SMetric3 &m2, const double t)
Definition: STensor3.cpp:260
intersection_conserve_mostaniso
SMetric3 intersection_conserve_mostaniso(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:216
meshMetric::_e_moins
double _e_moins
Definition: meshMetric.h:38
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
meshMetric::computeMetric
void computeMetric(int metricNumber)
Definition: meshMetric.cpp:748
MElement
Definition: MElement.h:30
meshMetric::hessians
std::map< MVertex *, SMetric3 > hessians
Definition: meshMetric.h:52
fullMatrix::eig
bool eig(fullVector< double > &eigenValReal, fullVector< double > &eigenValImag, fullMatrix< scalar > &leftEigenVect, fullMatrix< scalar > &rightEigenVect, bool sortRealPart=false)
Definition: fullMatrix.h:727
fullMatrix::gemm
void gemm(const fullMatrix< scalar > &a, const fullMatrix< scalar > &b, scalar alpha=1., scalar beta=1., bool transposeA=false, bool transposeB=false)
Definition: fullMatrix.h:580
S
#define S
Definition: DefaultOptions.h:21
MElementOctree::find
MElement * find(double x, double y, double z, int dim=-1, bool strict=false) const
Definition: MElementOctree.cpp:205
meshMetric::computeMetricHessian
void computeMetricHessian(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x=0.0, double y=0.0, double z=0.0)
Definition: meshMetric.cpp:406
meshMetric::getGradient
SVector3 getGradient(MVertex *v)
Definition: meshMetric.cpp:945
GModel::firstRegion
riter firstRegion()
Definition: GModel.h:354
GModel::lastRegion
riter lastRegion()
Definition: GModel.h:358
meshMetric::getLaplacian
double getLaplacian(MVertex *v)
Definition: meshMetric.cpp:937
Context.h
intersection_conserve_mostaniso_2d
SMetric3 intersection_conserve_mostaniso_2d(const SMetric3 &m1, const SMetric3 &m2)
Definition: STensor3.cpp:236
meshMetric::LEVELSET
@ LEVELSET
Definition: meshMetric.h:25
z
const double z
Definition: GaussQuadratureQuad.cpp:56
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
meshMetric::operator()
virtual double operator()(double x, double y, double z, GEntity *ge=nullptr)
Definition: meshMetric.cpp:792
fullMatrix::luSolve
bool luSolve(const fullVector< scalar > &rhs, fullVector< scalar > &result)
Definition: fullMatrix.h:603
meshMetric::FREY
@ FREY
Definition: meshMetric.h:27
meshMetric::~meshMetric
~meshMetric()
Definition: meshMetric.cpp:210
meshMetric::computeHessian
void computeHessian()
Definition: meshMetric.cpp:275
meshMetric::hmin
double hmin
Definition: meshMetric.h:42
GModel.h
MVertex::y
double y() const
Definition: MVertex.h:61
meshMetric::_nodalSizes
nodalField _nodalSizes
Definition: meshMetric.h:60
meshMetric::setOfTechniques
std::map< int, int > setOfTechniques
Definition: meshMetric.h:67
simpleFunction::hasDerivatives
virtual bool hasDerivatives()
Definition: simpleFunction.h:22
MElement::interpolate
double interpolate(double val[], double u, double v, double w, int stride=1, int order=-1)
Definition: MElement.cpp:1216
gmshLevelset.h
meshMetric::grads
std::map< MVertex *, SVector3 > grads
Definition: meshMetric.h:51
meshMetric::hmax
double hmax
Definition: meshMetric.h:42
meshMetric::addMetric
void addMetric(int technique, simpleFunction< double > *fct, const std::vector< double > &parameters)
Definition: meshMetric.cpp:64
meshMetric::computeValues
void computeValues()
Definition: meshMetric.cpp:216
MVertex::x
double x() const
Definition: MVertex.h:60
meshMetric::EIGENDIRECTIONS
@ EIGENDIRECTIONS
Definition: meshMetric.h:28
SVector3::normalize
double normalize()
Definition: SVector3.h:38
meshMetric::vals
std::map< MVertex *, double > vals
Definition: meshMetric.h:50
scale
static void scale(std::vector< double > &x, double s)
Definition: ConjugateGradients.cpp:21