1 #include <mystdlib.h>
2
3 #include "meshing.hpp"
4 #ifdef SOLIDGEOM
5 #include <csg.hpp>
6 #endif
7 #include <opti.hpp>
8 #include <core/array.hpp>
9 #include <core/taskmanager.hpp>
10
11
12 namespace netgen
13 {
14 using namespace ngcore;
15
16
Func(const Vector & x) const17 double MinFunctionSum :: Func (const Vector & x) const
18 {
19 double retval = 0;
20 for(int i=0; i<functions.Size(); i++)
21 retval += functions[i]->Func(x);
22
23 return retval;
24 }
25
Grad(const Vector & x,Vector & g) const26 void MinFunctionSum :: Grad (const Vector & x, Vector & g) const
27 {
28 g = 0.;
29 VectorMem<3> gi;
30 for(int i=0; i<functions.Size(); i++)
31 {
32 functions[i]->Grad(x,gi);
33 for(int j=0; j<g.Size(); j++)
34 g[j] += gi[j];
35 }
36 }
37
38
FuncGrad(const Vector & x,Vector & g) const39 double MinFunctionSum :: FuncGrad (const Vector & x, Vector & g) const
40 {
41 double retval = 0;
42 g = 0.;
43 VectorMem<3> gi;
44 for(int i=0; i<functions.Size(); i++)
45 {
46 retval += functions[i]->FuncGrad(x,gi);
47 for(int j=0; j<g.Size(); j++)
48 g[j] += gi[j];
49 }
50 return retval;
51 }
52
FuncDeriv(const Vector & x,const Vector & dir,double & deriv) const53 double MinFunctionSum :: FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
54 {
55 double retval = 0;
56 deriv = 0.;
57 double derivi;
58 for(int i=0; i<functions.Size(); i++)
59 {
60 retval += functions[i]->FuncDeriv(x,dir,derivi);
61 deriv += derivi;
62 }
63 return retval;
64 }
65
GradStopping(const Vector & x) const66 double MinFunctionSum :: GradStopping (const Vector & x) const
67 {
68 double minfs(0), mini;
69 for(int i=0; i<functions.Size(); i++)
70 {
71 mini = functions[i]->GradStopping(x);
72 if(i==0 || mini < minfs)
73 minfs = mini;
74 }
75 return minfs;
76 }
77
78
AddFunction(MinFunction & fun)79 void MinFunctionSum :: AddFunction(MinFunction & fun)
80 {
81 functions.Append(&fun);
82 }
83
Function(int i) const84 const MinFunction & MinFunctionSum :: Function(int i) const
85 {
86 return *functions[i];
87 }
Function(int i)88 MinFunction & MinFunctionSum :: Function(int i)
89 {
90 return *functions[i];
91 }
92
PointFunction1(Mesh::T_POINTS & apoints,const NgArray<INDEX_3> & afaces,const MeshingParameters & amp,double ah)93 PointFunction1 :: PointFunction1 (Mesh::T_POINTS & apoints,
94 const NgArray<INDEX_3> & afaces,
95 const MeshingParameters & amp,
96 double ah)
97 : points(apoints), faces(afaces), mp(amp)
98 {
99 h = ah;
100 }
101
102
Func(const Vector & vp) const103 double PointFunction1 :: Func (const Vector & vp) const
104 {
105 double badness = 0;
106 Point<3> pp(vp(0), vp(1), vp(2));
107
108 for (int j = 0; j < faces.Size(); j++)
109 {
110 const INDEX_3 & el = faces[j];
111
112 double bad = CalcTetBadness (points[PointIndex (el.I1())],
113 points[PointIndex (el.I3())],
114 points[PointIndex (el.I2())],
115 pp, 0, mp);
116 badness += bad;
117 }
118
119 return badness;
120 }
121
122
123 double PointFunction1 ::
FuncDeriv(const Vector & x,const Vector & dir,double & deriv) const124 FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
125 {
126 VectorMem<3> hx;
127 const double eps = 1e-6;
128
129 double dirlen = dir.L2Norm();
130 if (dirlen < 1e-14)
131 {
132 deriv = 0;
133 return Func(x);
134 }
135
136 hx.Set(1, x);
137 hx.Add(eps * h / dirlen, dir);
138 double fr = Func (hx);
139 hx.Set(1, x);
140 hx.Add(-eps * h / dirlen, dir);
141 double fl = Func (hx);
142
143 deriv = (fr - fl) / (2 * eps * h) * dirlen;
144
145 return Func(x);
146 }
147
148
FuncGrad(const Vector & x,Vector & g) const149 double PointFunction1 :: FuncGrad (const Vector & x, Vector & g) const
150 {
151 VectorMem<3> hx;
152 double eps = 1e-6;
153
154 hx = x;
155 for (int i = 0; i < 3; i++)
156 {
157 hx(i) = x(i) + eps * h;
158 double fr = Func (hx);
159 hx(i) = x(i) - eps * h;
160 double fl = Func (hx);
161 hx(i) = x(i);
162
163 g(i) = (fr - fl) / (2 * eps * h);
164 }
165
166 return Func(x);
167 }
168
GradStopping(const Vector & x) const169 double PointFunction1 :: GradStopping (const Vector & x) const
170 {
171 double f = Func(x);
172 return 1e-8 * f * f;
173 }
174
175
176 /* Cheap Functional depending of inner point inside triangular surface */
177
178 // is it used ????
179 class CheapPointFunction1 : public MinFunction
180 {
181 Mesh::T_POINTS & points;
182 const NgArray<INDEX_3> & faces;
183 DenseMatrix m;
184 double h;
185 public:
186 CheapPointFunction1 (Mesh::T_POINTS & apoints,
187 const NgArray<INDEX_3> & afaces,
188 double ah);
189
190 virtual double Func (const Vector & x) const;
191 virtual double FuncGrad (const Vector & x, Vector & g) const;
192 };
193
CheapPointFunction1(Mesh::T_POINTS & apoints,const NgArray<INDEX_3> & afaces,double ah)194 CheapPointFunction1 :: CheapPointFunction1 (Mesh::T_POINTS & apoints,
195 const NgArray<INDEX_3> & afaces,
196 double ah)
197 : points(apoints), faces(afaces)
198 {
199 h = ah;
200
201
202 int nf = faces.Size();
203
204 m.SetSize (nf, 4);
205
206 for (int i = 1; i <= nf; i++)
207 {
208 const Point3d & p1 = points[PointIndex(faces.Get(i).I1())];
209 const Point3d & p2 = points[PointIndex(faces.Get(i).I2())];
210 const Point3d & p3 = points[PointIndex(faces.Get(i).I3())];
211 Vec3d v1 (p1, p2);
212 Vec3d v2 (p1, p3);
213 Vec3d n;
214 Cross (v1, v2, n);
215 n /= n.Length();
216
217 m.Elem(i, 1) = n.X();
218 m.Elem(i, 2) = n.Y();
219 m.Elem(i, 3) = n.Z();
220 m.Elem(i, 4) = - (n.X() * p1.X() + n.Y() * p1.Y() + n.Z() * p1.Z());
221 }
222 }
223
Func(const Vector & vp) const224 double CheapPointFunction1 :: Func (const Vector & vp) const
225 {
226
227 /*
228 int j;
229 double badness = 0;
230 Point3d pp(vp.Get(1), vp.Get(2), vp.Get(3));
231
232 for (j = 1; j <= faces.Size(); j++)
233 {
234 const INDEX_3 & el = faces.Get(j);
235
236 double bad = CalcTetBadness (points.Get(el.I1()),
237 points.Get(el.I3()),
238 points.Get(el.I2()),
239 pp, 0);
240 badness += bad;
241 }
242 */
243
244 int i;
245 double badness = 0;
246 VectorMem<4> hv;
247 Vector res(m.Height());
248
249 for (i = 0;i < 3; i++)
250 hv(i) = vp(i);
251 hv(3) = 1;
252 m.Mult (hv, res);
253
254 for (i = 1; i <= res.Size(); i++)
255 {
256 if (res(i-1) < 1e-10)
257 badness += 1e24;
258 else
259 badness += 1 / res(i-1);
260 }
261
262 return badness;
263 }
264
265
FuncGrad(const Vector & x,Vector & g) const266 double CheapPointFunction1 :: FuncGrad (const Vector & x, Vector & g) const
267 {
268 VectorMem<3> hx;
269 double eps = 1e-6;
270
271 hx = x;
272 for (int i = 0; i < 3; i++)
273 {
274 hx(i) = x(i) + eps * h;
275 double fr = Func (hx);
276 hx(i) = x(i) - eps * h;
277 double fl = Func (hx);
278 hx(i) = x(i);
279
280 g(i) = (fr - fl) / (2 * eps * h);
281 }
282
283 return Func(x);
284 }
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299 /* ************* PointFunction **************************** */
300
301
302 class PointFunction
303 {
304 public:
305 Mesh::T_POINTS & points;
306 const Array<Element> & elements;
307 Table<int, PointIndex> &elementsonpoint;
308 bool own_elementsonpoint;
309 const MeshingParameters & mp;
310 PointIndex actpind;
311 double h;
312
313 public:
314 PointFunction (Mesh::T_POINTS & apoints,
315 const Array<Element> & aelements,
316 const MeshingParameters & amp);
317 PointFunction (const PointFunction & pf);
~PointFunction()318 virtual ~PointFunction () { if(own_elementsonpoint) delete &elementsonpoint; }
319 virtual void SetPointIndex (PointIndex aactpind);
SetLocalH(double ah)320 void SetLocalH (double ah) { h = ah; }
GetLocalH() const321 double GetLocalH () const { return h; }
GetPointToElementTable()322 const Table<int, PointIndex> & GetPointToElementTable() { return elementsonpoint; };
323 virtual double PointFunctionValue (const Point<3> & pp) const;
324 virtual double PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const;
325 virtual double PointFunctionValueDeriv (const Point<3> & pp, const Vec<3> & dir, double & deriv) const;
326
327 int MovePointToInner ();
328 };
329
330
PointFunction(const PointFunction & pf)331 PointFunction :: PointFunction (const PointFunction & pf)
332 : points(pf.points), elements(pf.elements), elementsonpoint(pf.elementsonpoint), own_elementsonpoint(false), mp(pf.mp)
333 { }
334
PointFunction(Mesh::T_POINTS & apoints,const Array<Element> & aelements,const MeshingParameters & amp)335 PointFunction :: PointFunction (Mesh::T_POINTS & apoints,
336 const Array<Element> & aelements,
337 const MeshingParameters & amp)
338 : points(apoints), elements(aelements), elementsonpoint(* new Table<int,PointIndex>()), own_elementsonpoint(true), mp(amp)
339 {
340 static Timer tim("PointFunction - build elementsonpoint table"); RegionTimer reg(tim);
341 elementsonpoint = std::move(ngcore::CreateSortedTable<int, PointIndex>( elements.Range(),
342 [&](auto & table, ElementIndex ei)
343 {
344 const auto & el = elements[ei];
345
346 if(el.NP()!=4)
347 return;
348
349 for (PointIndex pi : el.PNums())
350 table.Add (pi, ei);
351 }, points.Size()));
352 }
353
SetPointIndex(PointIndex aactpind)354 void PointFunction :: SetPointIndex (PointIndex aactpind)
355 {
356 actpind = aactpind;
357 }
358
PointFunctionValue(const Point<3> & pp) const359 double PointFunction :: PointFunctionValue (const Point<3> & pp) const
360 {
361 double badness;
362 Point<3> hp;
363
364 badness = 0;
365
366 hp = points[actpind];
367 points[actpind] = Point<3> (pp);
368
369 for (auto ei : elementsonpoint[actpind])
370 {
371 const Element & el = elements[ei];
372 badness += CalcTetBadness (points[el[0]], points[el[1]],
373 points[el[2]], points[el[3]], -1, mp);
374 }
375
376 points[actpind] = Point<3> (hp);
377 return badness;
378 }
379
380
PointFunctionValueGrad(const Point<3> & pp,Vec<3> & grad) const381 double PointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const
382 {
383 double f = 0;
384
385 Point<3> hp = points[actpind];
386 Vec<3> vgradi, vgrad(0,0,0);
387 points[actpind] = Point<3> (pp);
388
389 for (auto ei : elementsonpoint[actpind])
390 {
391 const Element & el = elements[ei];
392 for (int k = 0; k < 4; k++)
393 if (el[k] == actpind)
394 {
395 f += CalcTetBadnessGrad (points[el[0]], points[el[1]],
396 points[el[2]], points[el[3]],
397 -1, k+1, vgradi, mp);
398
399 vgrad += vgradi;
400 }
401 }
402
403 points[actpind] = Point<3> (hp);
404
405 grad = vgrad;
406 return f;
407 }
408
409
PointFunctionValueDeriv(const Point<3> & pp,const Vec<3> & dir,double & deriv) const410 double PointFunction :: PointFunctionValueDeriv (const Point<3> & pp, const Vec<3> & dir,
411 double & deriv) const
412 {
413 Vec<3> vgradi, vgrad(0,0,0);
414
415 Point<3> hp = points[actpind];
416 points[actpind] = pp;
417 double f = 0;
418
419 for (auto ei : elementsonpoint[actpind])
420 {
421 const Element & el = elements[ei];
422
423 for (int k = 1; k <= 4; k++)
424 if (el.PNum(k) == actpind)
425 {
426 f += CalcTetBadnessGrad (points[el.PNum(1)],
427 points[el.PNum(2)],
428 points[el.PNum(3)],
429 points[el.PNum(4)], -1, k, vgradi, mp);
430
431 vgrad += vgradi;
432 }
433 }
434
435 points[actpind] = Point<3> (hp);
436 deriv = dir * vgrad;
437 return f;
438 }
439
MovePointToInner()440 int PointFunction :: MovePointToInner ()
441 {
442 // try point movement
443 NgArray<Element2d> faces;
444
445 for (auto ei : elementsonpoint[actpind])
446 {
447 const Element & el = elements[ei];
448
449 for (int k = 1; k <= 4; k++)
450 if (el.PNum(k) == actpind)
451 {
452 Element2d face(TRIG);
453 el.GetFace (k, face);
454 Swap (face.PNum(2), face.PNum(3));
455 faces.Append (face);
456 }
457 }
458
459 Point3d hp;
460 int hi = FindInnerPoint (points, faces, hp);
461 if (hi)
462 {
463 // cout << "inner point found" << endl;
464 points[actpind] = Point<3> (hp);
465 }
466 else
467 ;
468 // cout << "no inner point found" << endl;
469
470 /*
471 Point3d hp2;
472 int hi2 = FindInnerPoint (points, faces, hp2);
473 if (hi2)
474 {
475 cout << "new: inner point found" << endl;
476 }
477 else
478 cout << "new: no inner point found" << endl;
479
480 (*testout) << "hi(orig) = " << hi << ", hi(new) = " << hi2;
481 if (hi != hi2) (*testout) << "hi different" << endl;
482 */
483
484 return hi;
485 }
486
487
488
489
490
491
492 class CheapPointFunction : public PointFunction
493 {
494 DenseMatrix m;
495 public:
496 CheapPointFunction (Mesh::T_POINTS & apoints,
497 const Array<Element> & aelements,
498 const MeshingParameters & amp);
499 virtual void SetPointIndex (PointIndex aactpind);
500 virtual double PointFunctionValue (const Point<3> & pp) const;
501 virtual double PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const;
502 };
503
504
CheapPointFunction(Mesh::T_POINTS & apoints,const Array<Element> & aelements,const MeshingParameters & amp)505 CheapPointFunction :: CheapPointFunction (Mesh::T_POINTS & apoints,
506 const Array<Element> & aelements,
507 const MeshingParameters & amp)
508 : PointFunction (apoints, aelements, amp)
509 {
510 ;
511 }
512
513
SetPointIndex(PointIndex aactpind)514 void CheapPointFunction :: SetPointIndex (PointIndex aactpind)
515 {
516 actpind = aactpind;
517
518 int ne = elementsonpoint[actpind].Size();
519 int i, j;
520 PointIndex pi1, pi2, pi3;
521
522 m.SetSize (ne, 4);
523
524 for (i = 0; i < ne; i++)
525 {
526 pi1 = 0;
527 pi2 = 0;
528 pi3 = 0;
529
530 const Element & el = elements[elementsonpoint[actpind][i]];
531 for (j = 1; j <= 4; j++)
532 if (el.PNum(j) != actpind)
533 {
534 pi3 = pi2;
535 pi2 = pi1;
536 pi1 = el.PNum(j);
537 }
538
539 const Point3d & p1 = points[pi1];
540 Vec3d v1 (p1, points[pi2]);
541 Vec3d v2 (p1, points[pi3]);
542 Vec3d n;
543 Cross (v1, v2, n);
544 n /= n.Length();
545
546 Vec3d v (p1, points[actpind]);
547 double c = v * n;
548
549 if (c < 0)
550 n *= -1;
551
552 // n is inner normal
553
554 m.Elem(i+1, 1) = n.X();
555 m.Elem(i+1, 2) = n.Y();
556 m.Elem(i+1, 3) = n.Z();
557 m.Elem(i+1, 4) = - (n.X() * p1.X() + n.Y() * p1.Y() + n.Z() * p1.Z());
558 }
559 }
560
PointFunctionValue(const Point<3> & pp) const561 double CheapPointFunction :: PointFunctionValue (const Point<3> & pp) const
562 {
563 VectorMem<4> p4;
564 Vector di;
565 int n = m.Height();
566
567 p4(0) = pp(0);
568 p4(1) = pp(1);
569 p4(2) = pp(2);
570 p4(3) = 1;
571
572 di.SetSize (n);
573 m.Mult (p4, di);
574
575 double sum = 0;
576 for (int i = 0; i < n; i++)
577 {
578 if (di(i) > 0)
579 sum += 1 / di(i);
580 else
581 return 1e16;
582 }
583 return sum;
584 }
585
586
587
588
PointFunctionValueGrad(const Point<3> & pp,Vec<3> & grad) const589 double CheapPointFunction :: PointFunctionValueGrad (const Point<3> & pp, Vec<3> & grad) const
590 {
591 VectorMem<4> p4;
592 Vector di;
593
594 int n = m.Height();
595
596 p4(0) = pp(0);
597 p4(1) = pp(1);
598 p4(2) = pp(2);
599 p4(3) = 1;
600
601 di.SetSize (n);
602 m.Mult (p4, di);
603
604 double sum = 0;
605 grad = 0;
606 for (int i = 0; i < n; i++)
607 {
608 if (di(i) > 0)
609 {
610 double idi = 1 / di(i);
611 sum += idi;
612 grad(0) -= idi * idi * m(i, 0);
613 grad(1) -= idi * idi * m(i, 1);
614 grad(2) -= idi * idi * m(i, 2);
615 }
616 else
617 {
618 return 1e16;
619 }
620 }
621 return sum;
622 }
623
624
625
626
627
628
629
630
631 class Opti3FreeMinFunction : public MinFunction
632 {
633 const PointFunction & pf;
634 Point<3> sp1;
635
636 public:
637 Opti3FreeMinFunction (const PointFunction & apf);
SetPoint(const Point<3> & asp1)638 void SetPoint (const Point<3> & asp1) { sp1 = asp1; }
639 virtual double Func (const Vector & x) const;
640 virtual double FuncGrad (const Vector & x, Vector & g) const;
641 virtual double FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const;
642 virtual double GradStopping (const Vector & x) const;
643 virtual void ApproximateHesse (const Vector & x,
644 DenseMatrix & hesse) const;
645 };
646
Opti3FreeMinFunction(const PointFunction & apf)647 Opti3FreeMinFunction :: Opti3FreeMinFunction (const PointFunction & apf)
648 : pf(apf)
649 {
650 ;
651 }
652
Func(const Vector & x) const653 double Opti3FreeMinFunction :: Func (const Vector & x) const
654 {
655 Point<3> pp;
656 for (int j = 0; j < 3; j++)
657 pp(j) = sp1(j) + x(j);
658 return pf.PointFunctionValue (pp);
659 }
660
FuncGrad(const Vector & x,Vector & grad) const661 double Opti3FreeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
662 {
663 Vec<3> vgrad;
664 Point<3> pp;
665
666 for (int j = 0; j < 3; j++)
667 pp(j) = sp1(j) + x(j);
668
669 double val = pf.PointFunctionValueGrad (pp, vgrad);
670
671 for (int j = 0; j < 3; j++)
672 grad(j) = vgrad(j);
673
674 return val;
675 }
676
FuncDeriv(const Vector & x,const Vector & dir,double & deriv) const677 double Opti3FreeMinFunction :: FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
678 {
679 Point<3> pp;
680
681 for (int j = 0; j < 3; j++)
682 pp(j) = sp1(j) + x(j);
683
684 Vec<3> vdir;
685 for (int j = 0; j < 3; j++)
686 vdir(j) = dir(j);
687
688 return pf.PointFunctionValueDeriv (pp, vdir, deriv);
689 }
690
GradStopping(const Vector & x) const691 double Opti3FreeMinFunction :: GradStopping (const Vector & x) const
692 {
693 double f = Func(x);
694 return 1e-3 * f / pf.GetLocalH();
695 }
696
697
ApproximateHesse(const Vector & x,DenseMatrix & hesse) const698 void Opti3FreeMinFunction :: ApproximateHesse (const Vector & x,
699 DenseMatrix & hesse) const
700 {
701 int n = x.Size();
702
703 Vector hx;
704 hx.SetSize(n);
705
706 double eps = 1e-8;
707 double f, f11, f22; //, f12, f21
708
709 f = Func(x);
710
711 for (int i = 1; i <= n; i++)
712 {
713 for (int j = 1; j < i; j++)
714 {
715 /*
716 hx = x;
717 hx.Elem(i) = x.Get(i) + eps;
718 hx.Elem(j) = x.Get(j) + eps;
719 f11 = Func(hx);
720 hx.Elem(i) = x.Get(i) + eps;
721 hx.Elem(j) = x.Get(j) - eps;
722 f12 = Func(hx);
723 hx.Elem(i) = x.Get(i) - eps;
724 hx.Elem(j) = x.Get(j) + eps;
725 f21 = Func(hx);
726 hx.Elem(i) = x.Get(i) - eps;
727 hx.Elem(j) = x.Get(j) - eps;
728 f22 = Func(hx);
729 */
730 hesse.Elem(i, j) = hesse.Elem(j, i) = 0;
731 // (f11 + f22 - f12 - f21) / (2 * eps * eps);
732 }
733
734 hx = x;
735 hx(i-1) = x(i-1) + eps;
736 f11 = Func(hx);
737 hx(i-1) = x(i-1) - eps;
738 f22 = Func(hx);
739
740 hesse.Elem(i, i) = (f11 + f22 - 2 * f) / (eps * eps) + 1e-12;
741 }
742 }
743
744
745
746
747
748
749 #ifdef SOLIDGEOM
750 class Opti3SurfaceMinFunction : public MinFunction
751 {
752 const PointFunction & pf;
753 Point3d sp1;
754 const Surface * surf;
755 Vec3d t1, t2;
756
757 public:
758 Opti3SurfaceMinFunction (const PointFunction & apf);
759
760 void SetPoint (const Surface * asurf, const Point3d & asp1);
761
762 void CalcNewPoint (const Vector & x, Point3d & np) const;
763 virtual double Func (const Vector & x) const;
764 virtual double FuncGrad (const Vector & x, Vector & g) const;
765 };
766
767
Opti3SurfaceMinFunction(const PointFunction & apf)768 Opti3SurfaceMinFunction :: Opti3SurfaceMinFunction (const PointFunction & apf)
769 : MinFunction(), pf(apf)
770 {
771 ;
772 }
773
SetPoint(const Surface * asurf,const Point3d & asp1)774 void Opti3SurfaceMinFunction :: SetPoint (const Surface * asurf, const Point3d & asp1)
775 {
776 Vec3d n;
777 sp1 = asp1;
778 surf = asurf;
779
780 Vec<3> hn;
781 surf -> GetNormalVector (sp1, hn);
782 n = hn;
783
784 n.GetNormal (t1);
785 t1 /= t1.Length();
786 t2 = Cross (n, t1);
787 }
788
789
CalcNewPoint(const Vector & x,Point3d & np) const790 void Opti3SurfaceMinFunction :: CalcNewPoint (const Vector & x,
791 Point3d & np) const
792 {
793 np.X() = sp1.X() + x.Get(1) * t1.X() + x.Get(2) * t2.X();
794 np.Y() = sp1.Y() + x.Get(1) * t1.Y() + x.Get(2) * t2.Y();
795 np.Z() = sp1.Z() + x.Get(1) * t1.Z() + x.Get(2) * t2.Z();
796
797 Point<3> hnp = np;
798 surf -> Project (hnp);
799 np = hnp;
800 }
801
802
Func(const Vector & x) const803 double Opti3SurfaceMinFunction :: Func (const Vector & x) const
804 {
805 Point3d pp1;
806
807 CalcNewPoint (x, pp1);
808 return pf.PointFunctionValue (pp1);
809 }
810
811
812
FuncGrad(const Vector & x,Vector & grad) const813 double Opti3SurfaceMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
814 {
815 Vec3d n, vgrad;
816 Point3d pp1;
817 VectorMem<3> freegrad;
818
819 CalcNewPoint (x, pp1);
820
821 double badness = pf.PointFunctionValueGrad (pp1, freegrad);
822 vgrad.X() = freegrad.Get(1);
823 vgrad.Y() = freegrad.Get(2);
824 vgrad.Z() = freegrad.Get(3);
825
826 Vec<3> hn;
827 surf -> GetNormalVector (pp1, hn);
828 n = hn;
829
830 vgrad -= (vgrad * n) * n;
831
832 grad.Elem(1) = vgrad * t1;
833 grad.Elem(2) = vgrad * t2;
834
835 return badness;
836 }
837 #endif
838
839
840
841
842
843
844
845
846 #ifdef SOLIDGEOM
847 class Opti3EdgeMinFunction : public MinFunction
848 {
849 const PointFunction & pf;
850 Point3d sp1;
851 const Surface *surf1, *surf2;
852 Vec3d t1;
853
854 public:
855 Opti3EdgeMinFunction (const PointFunction & apf);
856
857 void SetPoint (const Surface * asurf1, const Surface * asurf2,
858 const Point3d & asp1);
859 void CalcNewPoint (const Vector & x, Point3d & np) const;
860 virtual double FuncGrad (const Vector & x, Vector & g) const;
861 virtual double Func (const Vector & x) const;
862 };
863
Opti3EdgeMinFunction(const PointFunction & apf)864 Opti3EdgeMinFunction :: Opti3EdgeMinFunction (const PointFunction & apf)
865 : MinFunction(), pf(apf)
866 {
867 ;
868 }
869
SetPoint(const Surface * asurf1,const Surface * asurf2,const Point3d & asp1)870 void Opti3EdgeMinFunction :: SetPoint (const Surface * asurf1,
871 const Surface * asurf2,
872 const Point3d & asp1)
873 {
874 Vec3d n1, n2;
875 sp1 = asp1;
876 surf1 = asurf1;
877 surf2 = asurf2;
878
879 Vec<3> hn1, hn2;
880 surf1 -> GetNormalVector (sp1, hn1);
881 surf2 -> GetNormalVector (sp1, hn2);
882 n1 = hn1;
883 n2 = hn2;
884 t1 = Cross (n1, n2);
885 }
886
CalcNewPoint(const Vector & x,Point3d & np) const887 void Opti3EdgeMinFunction :: CalcNewPoint (const Vector & x,
888 Point3d & np) const
889 {
890 np.X() = sp1.X() + x.Get(1) * t1.X();
891 np.Y() = sp1.Y() + x.Get(1) * t1.Y();
892 np.Z() = sp1.Z() + x.Get(1) * t1.Z();
893 Point<3> hnp = np;
894 ProjectToEdge (surf1, surf2, hnp);
895 np = hnp;
896 }
897
Func(const Vector & x) const898 double Opti3EdgeMinFunction :: Func (const Vector & x) const
899 {
900 Vector g(x.Size());
901 return FuncGrad (x, g);
902 }
903
904
FuncGrad(const Vector & x,Vector & grad) const905 double Opti3EdgeMinFunction :: FuncGrad (const Vector & x, Vector & grad) const
906 {
907 Vec3d n1, n2, v1, vgrad;
908 Point3d pp1;
909 double badness;
910 VectorMem<3> freegrad;
911
912 CalcNewPoint (x, pp1);
913
914
915 badness = pf.PointFunctionValueGrad (pp1, freegrad);
916
917 vgrad.X() = freegrad.Get(1);
918 vgrad.Y() = freegrad.Get(2);
919 vgrad.Z() = freegrad.Get(3);
920
921 Vec<3> hn1, hn2;
922 surf1 -> GetNormalVector (pp1, hn1);
923 surf2 -> GetNormalVector (pp1, hn2);
924 n1 = hn1;
925 n2 = hn2;
926
927 v1 = Cross (n1, n2);
928 v1 /= v1.Length();
929
930 grad.Elem(1) = (vgrad * v1) * (t1 * v1);
931 return badness;
932 }
933 #endif
934
935
936
937
938
WrongOrientation(const Mesh::T_POINTS & points,const Element & el)939 int WrongOrientation (const Mesh::T_POINTS & points, const Element & el)
940 {
941 const Point3d & p1 = points[el.PNum(1)];
942 const Point3d & p2 = points[el.PNum(2)];
943 const Point3d & p3 = points[el.PNum(3)];
944 const Point3d & p4 = points[el.PNum(4)];
945
946 Vec3d v1(p1, p2);
947 Vec3d v2(p1, p3);
948 Vec3d v3(p1, p4);
949 Vec3d n;
950
951 Cross (v1, v2, n);
952 double vol = n * v3;
953
954 return (vol > 0);
955 }
956
957
958
959
960
961
962
963
964
965
966
967 /* ************* JacobianPointFunction **************************** */
968
969
970
971
972 // class JacobianPointFunction : public MinFunction
973 // {
974 // public:
975 // Mesh::T_POINTS & points;
976 // const NgArray<Element> & elements;
977 // TABLE<INDEX> elementsonpoint;
978 // PointIndex actpind;
979
980 // public:
981 // JacobianPointFunction (Mesh::T_POINTS & apoints,
982 // const NgArray<Element> & aelements);
983
984 // virtual void SetPointIndex (PointIndex aactpind);
985 // virtual double Func (const Vector & x) const;
986 // virtual double FuncGrad (const Vector & x, Vector & g) const;
987 // virtual double FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const;
988 // };
989
990
991 JacobianPointFunction ::
JacobianPointFunction(Mesh::T_POINTS & apoints,const Array<Element> & aelements)992 JacobianPointFunction (Mesh::T_POINTS & apoints,
993 const Array<Element> & aelements)
994 : points(apoints), elements(aelements), elementsonpoint(apoints.Size())
995 {
996 for (int i = 0; i < elements.Size(); i++)
997 for (int j = 1; j <= elements[i].NP(); j++)
998 elementsonpoint.Add1 (elements[i].PNum(j), i+1);
999
1000 onplane = false;
1001 }
1002
SetPointIndex(PointIndex aactpind)1003 void JacobianPointFunction :: SetPointIndex (PointIndex aactpind)
1004 {
1005 actpind = aactpind;
1006 }
1007
1008
Func(const Vector & v) const1009 double JacobianPointFunction :: Func (const Vector & v) const
1010 {
1011 int j;
1012 double badness = 0;
1013
1014 Point<3> hp = points[actpind];
1015
1016 points[actpind] = hp + Vec<3> (v(0), v(1), v(2));
1017
1018 if(onplane)
1019 points[actpind] -= (v(0)*nv(0)+v(1)*nv(1)+v(2)*nv(2)) * nv;
1020
1021
1022 for (auto eli : elementsonpoint[actpind])
1023 badness += elements[eli].CalcJacobianBadness (points);
1024
1025 points[actpind] = hp;
1026
1027 return badness;
1028 }
1029
1030
1031
1032
1033
1034 double JacobianPointFunction ::
FuncGrad(const Vector & x,Vector & g) const1035 FuncGrad (const Vector & x, Vector & g) const
1036 {
1037 int j, k;
1038 int lpi;
1039 double badness = 0;//, hbad;
1040
1041 Point<3> hp = points[actpind];
1042 points[actpind] = hp + Vec<3> (x(0), x(1), x(2));
1043
1044 if(onplane)
1045 points[actpind] -= (x(0)*nv(0)+x(1)*nv(1)+x(2)*nv(2)) * nv;
1046
1047 Vec<3> hderiv;
1048 //Vec3d vdir;
1049 g.SetSize(3);
1050 g = 0;
1051
1052 for (auto ei : elementsonpoint[actpind])
1053 {
1054 const Element & el = elements[ei];
1055
1056 lpi = 0;
1057 for (k = 1; k <= el.GetNP(); k++)
1058 if (el.PNum(k) == actpind)
1059 lpi = k;
1060 if (!lpi) cerr << "loc point not found" << endl;
1061
1062 badness += elements[ei].
1063 CalcJacobianBadnessGradient (points, lpi, hderiv);
1064
1065 for(k=0; k<3; k++)
1066 g(k) += hderiv(k);
1067
1068 /*
1069 for (k = 1; k <= 3; k++)
1070 {
1071 vdir = Vec3d(0,0,0);
1072 vdir.X(k) = 1;
1073
1074 hbad = elements.Get(eli).
1075 CalcJacobianBadnessDirDeriv (points, lpi, vdir, hderiv);
1076 //(*testout) << "hderiv " << k << ": " << hderiv << endl;
1077 g.Elem(k) += hderiv;
1078 if (k == 1)
1079 badness += hbad;
1080 }
1081 */
1082 }
1083
1084 if(onplane)
1085 {
1086 double scal = nv(0)*g(0) + nv(1)*g(1) + nv(2)*g(2);
1087 g(0) -= scal*nv(0);
1088 g(1) -= scal*nv(1);
1089 g(2) -= scal*nv(2);
1090 }
1091
1092 //(*testout) << "g = " << g << endl;
1093
1094
1095 points[actpind] = hp;
1096
1097 return badness;
1098 }
1099
1100
1101 double JacobianPointFunction ::
FuncDeriv(const Vector & x,const Vector & dir,double & deriv) const1102 FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const
1103 {
1104 int j, k;
1105 int lpi;
1106 double badness = 0;
1107
1108 Point<3> hp = points[actpind];
1109 points[actpind] = Point<3> (hp + Vec3d (x(0), x(1), x(2)));
1110
1111 if(onplane)
1112 points[actpind] -= (Vec3d (x(0), x(1), x(2))*nv) * nv;
1113
1114 double hderiv;
1115 deriv = 0;
1116 Vec<3> vdir(dir(0), dir(1), dir(2));
1117
1118 if(onplane)
1119 {
1120 double scal = vdir * nv;
1121 vdir -= scal*nv;
1122 }
1123
1124 for (auto ei : elementsonpoint[actpind])
1125 {
1126 const Element & el = elements[ei];
1127
1128 lpi = 0;
1129 for (k = 1; k <= el.GetNP(); k++)
1130 if (el.PNum(k) == actpind)
1131 lpi = k;
1132 if (!lpi) cerr << "loc point not found" << endl;
1133
1134 badness += elements[ei].
1135 CalcJacobianBadnessDirDeriv (points, lpi, vdir, hderiv);
1136 deriv += hderiv;
1137 }
1138
1139 points[actpind] = hp;
1140
1141 return badness;
1142
1143 }
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154 #ifdef SOLIDGEOMxxxx
ImproveMesh(const CSG eometry & geometry,OPTIMIZEGOAL goal)1155 void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal)
1156 {
1157 INDEX i, eli;
1158 int j;
1159 int typ = 1;
1160
1161 if (!&geometry || geometry.GetNSurf() == 0)
1162 {
1163 ImproveMesh (goal);
1164 return;
1165 }
1166
1167 const char * savetask = multithread.task;
1168 multithread.task = "Optimize Volume: Smooth Mesh";
1169
1170
1171 TABLE<INDEX> surfelementsonpoint(points.Size());
1172 Vector x(3), xsurf(2), xedge(1);
1173 int surf, surf1, surf2, surf3;
1174
1175 int uselocalh = mparam.uselocalh;
1176
1177 (*testout) << setprecision(8);
1178 (*testout) << "Improve Mesh" << "\n";
1179 PrintMessage (3, "ImproveMesh");
1180 // (*mycout) << "Vol = " << CalcVolume (points, volelements) << endl;
1181
1182
1183 for (i = 1; i <= surfelements.Size(); i++)
1184 for (j = 1; j <= 3; j++)
1185 surfelementsonpoint.Add1 (surfelements.Get(i).PNum(j), i);
1186
1187
1188 PointFunction * pf;
1189 if (typ == 1)
1190 pf = new PointFunction(points, volelements);
1191 else
1192 pf = new CheapPointFunction(points, volelements);
1193
1194 // pf->SetLocalH (h);
1195
1196 Opti3FreeMinFunction freeminf(*pf);
1197 Opti3SurfaceMinFunction surfminf(*pf);
1198 Opti3EdgeMinFunction edgeminf(*pf);
1199
1200 OptiParameters par;
1201 par.maxit_linsearch = 20;
1202 par.maxit_bfgs = 20;
1203
1204 int printmod = 1;
1205 char printdot = '.';
1206 if (points.Size() > 1000)
1207 {
1208 printmod = 10;
1209 printdot = '+';
1210 }
1211 if (points.Size() > 10000)
1212 {
1213 printmod = 100;
1214 printdot = '*';
1215 }
1216
1217 for (i = 1; i <= points.Size(); i++)
1218 {
1219 // if (ptyps.Get(i) == FIXEDPOINT) continue;
1220 if (ptyps.Get(i) != INNERPOINT) continue;
1221
1222 if (multithread.terminate)
1223 throw NgException ("Meshing stopped");
1224 /*
1225 if (multithread.terminate)
1226 break;
1227 */
1228 multithread.percent = 100.0 * i /points.Size();
1229
1230 /*
1231 if (points.Size() < 1000)
1232 PrintDot ();
1233 else
1234 if (i % 10 == 0)
1235 PrintDot ('+');
1236 */
1237 if (i % printmod == 0) PrintDot (printdot);
1238
1239 // (*testout) << "Now point " << i << "\n";
1240 // (*testout) << "Old: " << points.Get(i) << "\n";
1241
1242 pf->SetPointIndex (i);
1243
1244 // if (uselocalh)
1245 {
1246 double lh = GetH (points.Get(i));
1247 pf->SetLocalH (GetH (points.Get(i)));
1248 par.typx = lh / 10;
1249 // (*testout) << "lh(" << points.Get(i) << ") = " << lh << "\n";
1250 }
1251
1252 surf1 = surf2 = surf3 = 0;
1253
1254 for (j = 1; j <= surfelementsonpoint.EntrySize(i); j++)
1255 {
1256 eli = surfelementsonpoint.Get(i, j);
1257 int surfi = surfelements.Get(eli).GetIndex();
1258
1259 if (surfi)
1260 {
1261 surf = GetFaceDescriptor(surfi).SurfNr();
1262
1263 if (!surf1)
1264 surf1 = surf;
1265 else if (surf1 != surf)
1266 {
1267 if (!surf2)
1268 surf2 = surf;
1269 else if (surf2 != surf)
1270 surf3 = surf;
1271 }
1272 }
1273 else
1274 {
1275 surf1 = surf2 = surf3 = 1; // simulates corner point
1276 }
1277 }
1278
1279
1280 if (surf2 && !surf3)
1281 {
1282 // (*testout) << "On Edge" << "\n";
1283 /*
1284 xedge = 0;
1285 edgeminf.SetPoint (geometry.GetSurface(surf1),
1286 geometry.GetSurface(surf2),
1287 points.Elem(i));
1288 BFGS (xedge, edgeminf, par);
1289
1290 edgeminf.CalcNewPoint (xedge, points.Elem(i));
1291 */
1292 }
1293
1294 if (surf1 && !surf2)
1295 {
1296 // (*testout) << "In Surface" << "\n";
1297 /*
1298 xsurf = 0;
1299 surfminf.SetPoint (geometry.GetSurface(surf1),
1300 points.Get(i));
1301 BFGS (xsurf, surfminf, par);
1302
1303 surfminf.CalcNewPoint (xsurf, points.Elem(i));
1304 */
1305 }
1306
1307 if (!surf1)
1308 {
1309 // (*testout) << "In Volume" << "\n";
1310 x = 0;
1311 freeminf.SetPoint (points.Elem(i));
1312 // par.typx =
1313 BFGS (x, freeminf, par);
1314
1315 points.Elem(i).X() += x.Get(1);
1316 points.Elem(i).Y() += x.Get(2);
1317 points.Elem(i).Z() += x.Get(3);
1318 }
1319
1320 // (*testout) << "New Point: " << points.Elem(i) << "\n" << "\n";
1321
1322 }
1323 PrintDot ('\n');
1324 // (*mycout) << "Vol = " << CalcVolume (points, volelements) << endl;
1325
1326 multithread.task = savetask;
1327
1328 }
1329 #endif
1330
1331
1332
1333
ImproveMeshSequential(const MeshingParameters & mp,OPTIMIZEGOAL goal)1334 void Mesh :: ImproveMeshSequential (const MeshingParameters & mp, OPTIMIZEGOAL goal)
1335 {
1336 static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t);
1337
1338 (*testout) << "Improve Mesh" << "\n";
1339 PrintMessage (3, "ImproveMesh");
1340
1341 int np = GetNP();
1342 int ne = GetNE();
1343
1344
1345 if (goal == OPT_QUALITY)
1346 {
1347 double bad1 = CalcTotalBad (mp);
1348 (*testout) << "Total badness = " << bad1 << endl;
1349 PrintMessage (5, "Total badness = ", bad1);
1350 }
1351
1352 Vector x(3);
1353
1354 (*testout) << setprecision(8);
1355
1356 //int uselocalh = mparam.uselocalh;
1357
1358
1359 PointFunction pf(points, volelements, mp);
1360
1361 Opti3FreeMinFunction freeminf(pf);
1362
1363 OptiParameters par;
1364 par.maxit_linsearch = 20;
1365 par.maxit_bfgs = 20;
1366
1367 NgArray<double, PointIndex::BASE> pointh (points.Size());
1368
1369 if(lochfunc)
1370 {
1371 for (PointIndex pi : points.Range())
1372 pointh[pi] = GetH(points[pi]);
1373 }
1374 else
1375 {
1376 pointh = 0;
1377 for (Element & el : VolumeElements())
1378 {
1379 double h = pow(el.Volume(points),1./3.);
1380 for (PointIndex pi : el.PNums())
1381 if (h > pointh[pi])
1382 pointh[pi] = h;
1383 }
1384 }
1385
1386
1387 int printmod = 1;
1388 char printdot = '.';
1389 if (points.Size() > 1000)
1390 {
1391 printmod = 10;
1392 printdot = '+';
1393 }
1394 if (points.Size() > 10000)
1395 {
1396 printmod = 100;
1397 printdot = '*';
1398 }
1399
1400
1401 const char * savetask = multithread.task;
1402 multithread.task = "Optimize Volume: Smooth Mesh";
1403
1404 for (PointIndex pi : points.Range())
1405 if ( (*this)[pi].Type() == INNERPOINT )
1406 {
1407 if (multithread.terminate)
1408 throw NgException ("Meshing stopped");
1409
1410 multithread.percent = 100.0 * (pi+1-PointIndex::BASE) / points.Size();
1411
1412 if ( (pi+1-PointIndex::BASE) % printmod == 0) PrintDot (printdot);
1413
1414 double lh = pointh[pi];
1415 pf.SetLocalH (lh);
1416 par.typx = lh;
1417
1418 freeminf.SetPoint (points[pi]);
1419 pf.SetPointIndex (pi);
1420
1421 x = 0;
1422 int pok;
1423 pok = freeminf.Func (x) < 1e10;
1424
1425 if (!pok)
1426 {
1427 pok = pf.MovePointToInner ();
1428
1429 freeminf.SetPoint (points[pi]);
1430 pf.SetPointIndex (pi);
1431 }
1432
1433 if (pok)
1434 {
1435 //*testout << "start BFGS, pok" << endl;
1436 BFGS (x, freeminf, par);
1437 //*testout << "BFGS complete, pok" << endl;
1438 points[pi](0) += x(0);
1439 points[pi](1) += x(1);
1440 points[pi](2) += x(2);
1441 }
1442 }
1443 PrintDot ('\n');
1444
1445 multithread.task = savetask;
1446
1447 if (goal == OPT_QUALITY)
1448 {
1449 double bad1 = CalcTotalBad (mp);
1450 (*testout) << "Total badness = " << bad1 << endl;
1451 PrintMessage (5, "Total badness = ", bad1);
1452 }
1453 }
1454
ImproveMesh(const MeshingParameters & mp,OPTIMIZEGOAL goal)1455 void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal)
1456 {
1457 static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t);
1458 static Timer tcoloring("coloring");
1459 static Timer tcalcbadmax("Calc badmax");
1460 static Timer topt("optimize");
1461 static Timer trange("range");
1462 static Timer tloch("loch");
1463
1464 // return ImproveMeshSequential(mp, goal);
1465 BuildBoundaryEdges(false);
1466
1467 (*testout) << "Improve Mesh" << "\n";
1468 PrintMessage (3, "ImproveMesh");
1469
1470 int np = GetNP();
1471 int ne = GetNE();
1472
1473 PointFunction pf_glob(points, volelements, mp);
1474
1475 auto & elementsonpoint = pf_glob.GetPointToElementTable();
1476
1477 const auto & getDofs = [&] (int i)
1478 {
1479 i += PointIndex::BASE;
1480 return FlatArray<int>(elementsonpoint[i].Size(), elementsonpoint[i].Data());
1481 };
1482
1483 Array<int> colors(points.Size());
1484
1485 tcoloring.Start();
1486 int ncolors = ngcore::ComputeColoring( colors, ne, getDofs );
1487 auto color_table = CreateTable<PointIndex, int>( points.Size(),
1488 [&] ( auto & table, int i )
1489 {
1490 PointIndex pi = i+static_cast<int>(PointIndex::BASE);
1491 table.Add(colors[i], pi);
1492 }, ncolors);
1493
1494 tcoloring.Stop();
1495
1496 if (goal == OPT_QUALITY)
1497 {
1498 double bad1 = CalcTotalBad (mp);
1499 (*testout) << "Total badness = " << bad1 << endl;
1500 PrintMessage (5, "Total badness = ", bad1);
1501 }
1502
1503
1504 (*testout) << setprecision(8);
1505
1506 Array<double, PointIndex> pointh (points.Size());
1507
1508 if(lochfunc)
1509 {
1510 RegionTimer rt(tloch);
1511 ParallelForRange(points.Range(), [&] (auto myrange)
1512 {
1513 for(auto pi : myrange)
1514 pointh[pi] = GetH(points[pi]);
1515 });
1516 }
1517 else
1518 {
1519 pointh = 0;
1520 for (Element & el : VolumeElements())
1521 {
1522 double h = pow(el.Volume(points),1./3.);
1523 for (PointIndex pi : el.PNums())
1524 if (h > pointh[pi])
1525 pointh[pi] = h;
1526 }
1527 }
1528
1529 const char * savetask = multithread.task;
1530 multithread.task = "Optimize Volume: Smooth Mesh";
1531
1532 topt.Start();
1533 int counter = 0;
1534 for (auto icolor : Range(ncolors))
1535 {
1536 if (multithread.terminate)
1537 throw NgException ("Meshing stopped");
1538
1539 ParallelForRange( color_table[icolor].Range(), [&](auto myrange)
1540 {
1541 RegionTracer reg(ngcore::TaskManager::GetThreadId(), trange, myrange.Size());
1542 Vector x(3);
1543
1544 PointFunction pf{pf_glob};
1545
1546 Opti3FreeMinFunction freeminf(pf);
1547
1548 OptiParameters par;
1549 par.maxit_linsearch = 20;
1550 par.maxit_bfgs = 20;
1551
1552 for (auto i : myrange)
1553 {
1554 PointIndex pi = color_table[icolor][i];
1555 if ( (*this)[pi].Type() == INNERPOINT )
1556 {
1557 counter++;
1558
1559 double lh = pointh[pi];
1560 pf.SetLocalH (lh);
1561 par.typx = lh;
1562
1563 freeminf.SetPoint (points[pi]);
1564 pf.SetPointIndex (pi);
1565
1566 x = 0;
1567 int pok;
1568 pok = freeminf.Func (x) < 1e10;
1569
1570 if (!pok)
1571 {
1572 pok = pf.MovePointToInner ();
1573
1574 freeminf.SetPoint (points[pi]);
1575 pf.SetPointIndex (pi);
1576 }
1577
1578 if (pok)
1579 {
1580 //*testout << "start BFGS, pok" << endl;
1581 BFGS (x, freeminf, par);
1582 //*testout << "BFGS complete, pok" << endl;
1583 points[pi](0) += x(0);
1584 points[pi](1) += x(1);
1585 points[pi](2) += x(2);
1586 }
1587 }
1588 }
1589 }, 4*ngcore::TaskManager::GetNumThreads());
1590 }
1591 topt.Stop();
1592
1593 multithread.task = savetask;
1594
1595 if (goal == OPT_QUALITY)
1596 {
1597 double bad1 = CalcTotalBad (mp);
1598 (*testout) << "Total badness = " << bad1 << endl;
1599 PrintMessage (5, "Total badness = ", bad1);
1600 }
1601 }
1602
1603
1604
1605 // Improve Condition number of Jacobian, any elements
ImproveMeshJacobian(const MeshingParameters & mp,OPTIMIZEGOAL goal,const NgBitArray * usepoint)1606 void Mesh :: ImproveMeshJacobian (const MeshingParameters & mp,
1607 OPTIMIZEGOAL goal, const NgBitArray * usepoint)
1608 {
1609 // int i, j;
1610
1611 (*testout) << "Improve Mesh Jacobian" << "\n";
1612 PrintMessage (3, "ImproveMesh Jacobian");
1613
1614 int np = GetNP();
1615 int ne = GetNE();
1616
1617
1618 Vector x(3);
1619
1620 (*testout) << setprecision(8);
1621
1622 JacobianPointFunction pf(points, volelements);
1623
1624
1625 OptiParameters par;
1626 par.maxit_linsearch = 20;
1627 par.maxit_bfgs = 20;
1628
1629 NgBitArray badnodes(np);
1630 badnodes.Clear();
1631
1632 for (int i = 1; i <= ne; i++)
1633 {
1634 const Element & el = VolumeElement(i);
1635 double bad = el.CalcJacobianBadness (Points());
1636 if (bad > 1)
1637 for (int j = 1; j <= el.GetNP(); j++)
1638 badnodes.Set (el.PNum(j));
1639 }
1640
1641 NgArray<double, PointIndex::BASE, PointIndex> pointh (points.Size());
1642
1643 if(lochfunc)
1644 {
1645 // for(i = 1; i<=points.Size(); i++)
1646 for (PointIndex pi : points.Range())
1647 pointh[pi] = GetH(points[pi]);
1648 }
1649 else
1650 {
1651 pointh = 0;
1652 for (int i=0; i<GetNE(); i++)
1653 {
1654 const Element & el = VolumeElement(i+1);
1655 double h = pow(el.Volume(points),1./3.);
1656 for(int j=1; j<=el.GetNV(); j++)
1657 if(h > pointh[el.PNum(j)])
1658 pointh[el.PNum(j)] = h;
1659 }
1660 }
1661
1662
1663
1664 const char * savetask = multithread.task;
1665 multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
1666
1667 // for (PointIndex pi = points.Begin(); i < points.End(); pi++)
1668 for (PointIndex pi : points.Range())
1669 {
1670 if ((*this)[pi].Type() != INNERPOINT)
1671 continue;
1672
1673 if(usepoint && !usepoint->Test(pi))
1674 continue;
1675
1676 //(*testout) << "improvejac, p = " << i << endl;
1677
1678 if (goal == OPT_WORSTCASE && !badnodes.Test(pi))
1679 continue;
1680 // (*testout) << "smooth p " << i << endl;
1681
1682 /*
1683 if (multithread.terminate)
1684 break;
1685 */
1686 if (multithread.terminate)
1687 throw NgException ("Meshing stopped");
1688
1689 multithread.percent = 100.0 * pi / points.Size();
1690
1691 if (points.Size() < 1000)
1692 PrintDot ();
1693 else
1694 if (pi % 10 == 0)
1695 PrintDot ('+');
1696
1697 double lh = pointh[pi];
1698 par.typx = lh;
1699
1700 pf.SetPointIndex (pi);
1701
1702 x = 0;
1703 int pok = (pf.Func (x) < 1e10);
1704
1705 if (pok)
1706 {
1707 //*testout << "start BFGS, Jacobian" << endl;
1708 BFGS (x, pf, par);
1709 //*testout << "end BFGS, Jacobian" << endl;
1710 points[pi](0) += x(0);
1711 points[pi](1) += x(1);
1712 points[pi](2) += x(2);
1713 }
1714 else
1715 {
1716 cout << "el not ok" << endl;
1717 }
1718 }
1719 PrintDot ('\n');
1720
1721
1722 multithread.task = savetask;
1723 }
1724
1725
1726
1727
1728 // Improve Condition number of Jacobian, any elements
ImproveMeshJacobianOnSurface(const MeshingParameters & mp,const NgBitArray & usepoint,const NgArray<Vec<3> * > & nv,OPTIMIZEGOAL goal,const NgArray<NgArray<int,PointIndex::BASE> * > * idmaps)1729 void Mesh :: ImproveMeshJacobianOnSurface (const MeshingParameters & mp,
1730 const NgBitArray & usepoint,
1731 const NgArray< Vec<3>* > & nv,
1732 OPTIMIZEGOAL goal,
1733 const NgArray< NgArray<int,PointIndex::BASE>* > * idmaps)
1734 {
1735 // int i, j;
1736
1737 (*testout) << "Improve Mesh Jacobian" << "\n";
1738 PrintMessage (3, "ImproveMesh Jacobian");
1739
1740 int np = GetNP();
1741 int ne = GetNE();
1742
1743
1744 Vector x(3);
1745
1746 (*testout).precision(8);
1747
1748 JacobianPointFunction pf(points, volelements);
1749
1750 NgArray< NgArray<int,PointIndex::BASE>* > locidmaps;
1751 const NgArray< NgArray<int,PointIndex::BASE>* > * used_idmaps;
1752
1753 if(idmaps)
1754 used_idmaps = idmaps;
1755 else
1756 {
1757 used_idmaps = &locidmaps;
1758
1759 for(int i=1; i<=GetIdentifications().GetMaxNr(); i++)
1760 {
1761 if(GetIdentifications().GetType(i) == Identifications::PERIODIC)
1762 {
1763 locidmaps.Append(new NgArray<int,PointIndex::BASE>);
1764 GetIdentifications().GetMap(i,*locidmaps.Last(),true);
1765 }
1766 }
1767 }
1768
1769
1770 bool usesum = (used_idmaps->Size() > 0);
1771 MinFunctionSum pf_sum;
1772
1773 JacobianPointFunction * pf2ptr = NULL;
1774 if(usesum)
1775 {
1776 pf2ptr = new JacobianPointFunction(points, volelements);
1777 pf_sum.AddFunction(pf);
1778 pf_sum.AddFunction(*pf2ptr);
1779 }
1780
1781
1782 OptiParameters par;
1783 par.maxit_linsearch = 20;
1784 par.maxit_bfgs = 20;
1785
1786 NgBitArray badnodes(np);
1787 badnodes.Clear();
1788
1789 for (int i = 1; i <= ne; i++)
1790 {
1791 const Element & el = VolumeElement(i);
1792 double bad = el.CalcJacobianBadness (Points());
1793 if (bad > 1)
1794 for (int j = 1; j <= el.GetNP(); j++)
1795 badnodes.Set (el.PNum(j));
1796 }
1797
1798 NgArray<double, PointIndex::BASE> pointh (points.Size());
1799
1800 if(lochfunc)
1801 {
1802 // for(i=1; i<=points.Size(); i++)
1803 for (PointIndex pi : points.Range())
1804 pointh[pi] = GetH(points[pi]);
1805 }
1806 else
1807 {
1808 pointh = 0;
1809 for(int i=0; i<GetNE(); i++)
1810 {
1811 const Element & el = VolumeElement(i+1);
1812 double h = pow(el.Volume(points),1./3.);
1813 for(int j=1; j<=el.GetNV(); j++)
1814 if(h > pointh[el.PNum(j)])
1815 pointh[el.PNum(j)] = h;
1816 }
1817 }
1818
1819
1820 const char * savetask = multithread.task;
1821 multithread.task = "Optimize Volume: Smooth Mesh Jacobian";
1822
1823 // for (PointIndex pi = points.Begin(); pi <= points.End(); pi++)
1824 for (PointIndex pi : points.Range())
1825 if ( usepoint.Test(pi) )
1826 {
1827 //(*testout) << "improvejac, p = " << i << endl;
1828
1829 if (goal == OPT_WORSTCASE && !badnodes.Test(pi))
1830 continue;
1831 // (*testout) << "smooth p " << i << endl;
1832
1833 /*
1834 if (multithread.terminate)
1835 break;
1836 */
1837 if (multithread.terminate)
1838 throw NgException ("Meshing stopped");
1839
1840 multithread.percent = 100.0 * pi / points.Size();
1841
1842 if (points.Size() < 1000)
1843 PrintDot ();
1844 else
1845 if (pi % 10 == 0)
1846 PrintDot ('+');
1847
1848 double lh = pointh[pi];//GetH(points.Get(i));
1849 par.typx = lh;
1850
1851 pf.SetPointIndex (pi);
1852
1853 PointIndex brother (-1);
1854 if(usesum)
1855 {
1856 for(int j=0; brother == -1 && j<used_idmaps->Size(); j++)
1857 {
1858 if(pi < (*used_idmaps)[j]->Size() + PointIndex::BASE)
1859 {
1860 brother = (*(*used_idmaps)[j])[pi];
1861 if(brother == pi || brother == 0)
1862 brother = -1;
1863 }
1864 }
1865 if(brother >= pi)
1866 {
1867 pf2ptr->SetPointIndex(brother);
1868 pf2ptr->SetNV(*nv[brother-1]);
1869 }
1870 }
1871
1872 if(usesum && brother < pi)
1873 continue;
1874
1875 //pf.UnSetNV(); x = 0;
1876 //(*testout) << "before " << pf.Func(x);
1877
1878 pf.SetNV(*nv[pi-1]);
1879
1880 x = 0;
1881 int pok = (brother == -1) ? (pf.Func (x) < 1e10) : (pf_sum.Func (x) < 1e10);
1882
1883 if (pok)
1884 {
1885
1886 if(brother == -1)
1887 BFGS (x, pf, par);
1888 else
1889 BFGS (x, pf_sum, par);
1890
1891
1892 for(int j=0; j<3; j++)
1893 points[pi](j) += x(j);// - scal*nv[i-1].X(j);
1894
1895 if(brother != -1)
1896 for(int j=0; j<3; j++)
1897 points[brother](j) += x(j);// - scal*nv[brother-1].X(j);
1898
1899
1900 }
1901 else
1902 {
1903 cout << "el not ok" << endl;
1904 (*testout) << "el not ok" << endl
1905 << " func " << ((brother == -1) ? pf.Func(x) : pf_sum.Func (x)) << endl;
1906 if(brother != -1)
1907 (*testout) << " func1 " << pf.Func(x) << endl
1908 << " func2 " << pf2ptr->Func(x) << endl;
1909 }
1910 }
1911
1912 PrintDot ('\n');
1913
1914 delete pf2ptr;
1915 for(int i=0; i<locidmaps.Size(); i++)
1916 delete locidmaps[i];
1917
1918 multithread.task = savetask;
1919 }
1920
1921
1922
1923
1924 }
1925