1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : ...
21 // E-MAIL : ...
22
23 // Implementation of P1-P0 FVM-FEM
24
25 // compile and link with ./load.link mat\_dervieux.cpp
26
27 /* clang-format off */
28 //ff-c++-LIBRARY-dep:
29 //ff-c++-cpp-dep:
30 /* clang-format on */
31 #define _USE_MATH_DEFINES
32 #include <iostream>
33 #include <cfloat>
34 #include <cmath>
35 using namespace std;
36 #include "ff++.hpp"
37
38 using namespace std;
39
40 R Min(const R a, const R b);
41 R Max(const R a, const R b);
Exchange(R a,R b)42 inline void Exchange(R a, R b) {
43 R c = a;
44 a = b;
45 b = c;
46 };
47 // const R precision = 1e-10;
48
49 // int LireTaille (const char *NomDuFichier, int &nbnoeuds);
50 // int Lire (const char *NomDuFichier, int n, R2 noeuds []);
51 template< typename T >
from_string(const string & Str,T & Dest)52 bool from_string(const string &Str, T &Dest) {
53 // cr�er un flux � partir de la chaine donn�e
54 istringstream iss(Str);
55
56 // tenter la conversion vers Dest
57 iss >> Dest;
58 return !Dest;
59 };
60
metrique(int nbpoints,R2 * Point,R & A,R & B,R & C,R epsilon)61 void metrique(int nbpoints, R2 *Point, R &A, R &B, R &C, R epsilon) {
62 C = 0.;
63
64 R epsilon0 = 1e-5, precision = 1e-15, delta = 1e-10;
65 R inf = 1e50;
66 R Rmin = inf, Rmax = 0.;
67 int indiceX0 = 0;
68 R2 *PPoint = new R2[nbpoints];
69
70 for (int i = 0; i < nbpoints; i++) {
71 Rmax = Max(Rmax, Point[i].norme( ));
72
73 // ---d�placement des points situ�es sur les axes--------------
74 if (abs(Point[i].x) <= precision) {
75 if (Point[i].y < 0) {
76 Point[i].x = -delta;
77 Point[i].y = -sqrt(pow(Point[i].y, 2) - pow(Point[i].x, 2));
78 }
79
80 if (Point[i].y > 0) {
81 Point[i].x = delta;
82 Point[i].y = sqrt(pow(Point[i].y, 2) - pow(Point[i].x, 2));
83 }
84 }
85
86 if (abs(Point[i].y) <= precision) {
87 if (Point[i].x < 0) {
88 Point[i].y = -delta;
89 Point[i].x = -sqrt(pow(Point[i].x, 2) - pow(Point[i].y, 2));
90 }
91
92 if (Point[i].x > 0) {
93 Point[i].y = delta;
94 Point[i].x = sqrt(pow(Point[i].x, 2) - pow(Point[i].y, 2));
95 }
96 }
97
98 // -----------------------------------------------------------
99 assert(abs(Point[i].x * Point[i].y) >= pow(precision, 2));
100
101 if (Rmin > Point[i].norme( )) {
102 indiceX0 = i;
103 Rmin = Point[i].norme( );
104 }
105
106 // cout<<Point[i]<<endl;
107 }
108
109 // -------permutation des indices de la liste des points :
110 // ranger la liste en commen�ant par le point X0-------
111 for (int k = 0; k < nbpoints - indiceX0; k++) {
112 PPoint[k] = Point[k + indiceX0];
113 }
114
115 for (int k = nbpoints - indiceX0; k < nbpoints; k++) {
116 PPoint[k] = Point[k - nbpoints + indiceX0];
117 }
118
119 for (int i = 0; i < nbpoints; i++) {
120 Point[i] = PPoint[i];
121 }
122
123 // ----------------------------------------------------------------
124
125 R X0;
126 R Y0;
127 R bmin = 0., bmax = inf, b1, b2, aik = 0., bik = 0., cik = 0.;
128 R Xk = 0., Yk = 0., Ck = 0., Ak = 0., Bk = 0., Xi = 0., Yi = 0., ri, detXY = 0., Ri, R0, r0;
129
130 X0 = Point[0].x;
131 Y0 = Point[0].y;
132 r0 = Point[0].norme( );
133 assert(r0 == Rmin);
134
135 R EPS = 0.; // pour recuperer la valeur de epsilon0 optimale
136 R epsilonmax = r0 * (1. - r0 / Rmax) / 20.;
137 R Tabepsilon[20];
138 int neps = 4;
139 // --------- discertisation de epsilon0----------------------------------
140 if (epsilonmax > 1e-2) {
141 neps = 10;
142
143 Tabepsilon[0] = 1e-5;
144 Tabepsilon[1] = 1e-4;
145 Tabepsilon[2] = 1e-3;
146
147 for (int i = 3; i < neps; i++) {
148 Tabepsilon[i] = (i - 3) * (epsilonmax - 1e-2) / (neps - 4.) + 1e-2;
149 }
150 } else {
151 Tabepsilon[0] = 1e-5;
152 Tabepsilon[1] = 1e-4;
153 Tabepsilon[2] = 1e-3;
154 Tabepsilon[3] = 1e-2;
155 }
156
157 // ------------------------------------------------------------------------
158
159 if (r0 <= epsilon0) {
160 epsilon0 = r0 * epsilon0;
161 }
162
163 B = A = 1. / ((r0 - epsilon0) * (r0 - epsilon0));
164 R epsilon0min = epsilon0;
165
166 if (abs(Rmin - Rmax) > 1e-5) {
167 int condition = -1;
168
169 for (int ee = 0; ee < neps - 1; ee++) { // boucle sur epsilon0---------------
170 epsilon0 = Tabepsilon[ee];
171 if (r0 <= epsilon0) {
172 epsilon0 = r0 * epsilon0;
173 }
174
175 assert(r0 > epsilon0);
176 R0 = r0 / (r0 - epsilon0);
177
178 for (int i = 1; i < nbpoints; i++) { // boucle sur chaque noeud
179 Xi = Point[i].x;
180 Yi = Point[i].y;
181 ri = Point[i].norme( );
182 if (ri <= epsilon) {
183 epsilon = ri * epsilon;
184 }
185
186 assert(ri > epsilon);
187 Ri = ri / (ri - epsilon);
188 detXY = Xi * Y0 - Yi * X0;
189
190 // ------deplacement des points align�s avec l'origine et X0-----------
191 if (abs(detXY) <= precision) {
192 Xi += delta;
193
194 if (Yi < 0) {
195 Yi = -sqrt(pow(ri, 2) - pow(Xi, 2));
196 } else {
197 Yi = sqrt(pow(ri, 2) - pow(Xi, 2));
198 }
199
200 Point[i].x = Xi;
201 Point[i].y = Yi;
202 ri = Point[i].norme( );
203 if (ri <= epsilon) {
204 epsilon = ri * epsilon;
205 }
206
207 assert(ri > epsilon);
208 Ri = ri / (ri - epsilon);
209 }
210
211 detXY = Xi * Y0 - Yi * X0;
212
213 assert(abs(detXY) >= precision);
214
215 // -----------------------------------------------------------------
216
217 // -----racines du polynome en b � minimiser----------------------------
218
219 R bb1 =
220 (1. / pow(detXY, 2)) *
221 (pow(X0 * Ri, 2) + pow(Xi * R0, 2) -
222 2. * abs(Xi * X0) * sqrt(pow(R0 * Ri, 2) - pow(detXY / (Rmax * (r0 - epsilon0)), 2)));
223 R bb2 =
224 (1. / pow(detXY, 2)) *
225 (pow(X0 * Ri, 2) + pow(Xi * R0, 2) +
226 2. * abs(Xi * X0) * sqrt(pow(R0 * Ri, 2) - pow(detXY / (Rmax * (r0 - epsilon0)), 2)));
227
228 // --fin----racines du polynome en b � minimiser--------------------
229
230 bmax = Min(bb2, pow(Rmax / pow((r0), 2), 2));
231
232 bmin = Max(1. / (Rmax * Rmax), bb1); // minoration de b
233 R Cte = Max(1e-9, (bmax - bmin) * 1e-9);
234 bmin = bmin * (1. + Cte);
235 bmax = bmax * (1. - Cte);
236
237 // bornes de b-----------------------------------------------------------
238
239 // cas: majoration de c --------------------------------------------
240 R Li = X0 * Xi * (pow(Rmax / pow(r0 - epsilon0min, 2), 2) - 1. / pow(Rmax, 2)) +
241 (pow(Ri * X0, 2) - pow(R0 * Xi, 2)) / detXY;
242 R LiXY = Xi * Y0 + Yi * X0;
243
244 if (abs(LiXY) >= precision) {
245 condition = 1;
246
247 if (Xi * X0 > 0) {
248 if (LiXY > 0) {
249 bmin = Max(bmin, -Li / LiXY);
250 } else {
251 bmax = Min(bmax, -Li / LiXY);
252 }
253 } else {
254 if (LiXY < 0) {
255 bmin = Max(bmin, -Li / LiXY);
256 } else {
257 bmax = Min(bmax, -Li / LiXY);
258 }
259 }
260 } else {
261 if (Li < 0) {
262 condition = 0;
263 } else {
264 condition = 1;
265 }
266 }
267
268 // cas minoration de c --------------------------------------------
269 Li = X0 * Xi * (-pow(Rmax / pow(r0 - epsilon0min, 2), 2) + 1. / pow(Rmax, 2)) +
270 (pow(Ri * X0, 2) - pow(R0 * Xi, 2)) / detXY;
271 LiXY = Xi * Y0 + Yi * X0;
272
273 if (abs(LiXY) >= precision) {
274 condition = 1;
275
276 if (Xi * X0 > 0) {
277 if (LiXY < 0) {
278 bmin = Max(bmin, -Li / LiXY);
279 } else {
280 bmax = Min(bmax, -Li / LiXY);
281 }
282 } else {
283 if (LiXY > 0) {
284 bmin = Max(bmin, -Li / LiXY);
285 } else {
286 bmax = Min(bmax, -Li / LiXY);
287 }
288 }
289 } else {
290 if (Li > 0) {
291 condition = 0;
292 } else {
293 condition = 1;
294 }
295 }
296
297 if (condition == 1) {
298 int test = -1;
299 // --cas : minoration de a-----------------------------------------------
300
301 R Gi =
302 ((Xi * Yi * R0 * R0 - X0 * Y0 * Ri * Ri) / detXY + Xi * X0 / (Rmax * Rmax)) / (Yi * Y0);
303
304 if (Xi * X0 > 0) {
305 if (Yi * Y0 > 0) {
306 bmin = Max(bmin, Gi);
307 } else {
308 bmax = Min(bmax, Gi);
309 }
310 } else {
311 if (Yi * Y0 < 0) {
312 bmin = Max(bmin, Gi);
313 } else {
314 bmax = Min(bmax, Gi);
315 }
316 }
317
318 // cas :majoration de a------------------------------------------------
319
320 R Hi = (Xi * X0 * Rmax * Rmax / pow((r0 - epsilon0min), 4) +
321 (Xi * Yi * R0 * R0 - X0 * Y0 * Ri * Ri) / detXY) /
322 (Yi * Y0);
323 if (Xi * X0 > 0) {
324 if (Yi * Y0 > 0) {
325 bmax = Min(bmax, Hi);
326 } else {
327 bmin = Max(bmin, Hi);
328 }
329 } else {
330 if (Yi * Y0 < 0) {
331 bmax = Min(bmax, Hi);
332 } else {
333 bmin = Max(bmin, Hi);
334 }
335 }
336
337 // ------fin bornes de b------------------------------------------------
338 b2 = bmax;
339 b1 = bmin;
340
341 for (int k = 1; k < nbpoints; k++) { // on balaye les contraintes
342 Xk = Point[k].x;
343 Yk = Point[k].y;
344 Bk = (Yk * Yk * Xi * X0 + Xk * (Xk * Yi * Y0 - Yk * (Yi * X0 + Xi * Y0))) / (Xi * X0);
345 Ck = (X0 * Xi * detXY -
346 Xk * (Xi * R0 * R0 * (Yk * Xi - Yi * Xk) + X0 * Ri * Ri * (-Yk * X0 + Y0 * Xk))) /
347 (Xi * X0 * detXY);
348
349 assert(abs(Xi * X0 * Y0 * Yi * Xk * Yk) >= pow(precision, 5));
350 if (abs(Bk) > precision) { // non nul
351 if (Bk <= 0) {
352 bmax = Min(bmax, Ck / Bk);
353 } else {
354 bmin = Max(bmin, Ck / Bk);
355 }
356
357 if ((bmax < b1) || (bmin > b2) || (bmin > bmax)) {
358 // cout<<" i = "<<i<<" k = "<<k<<endl;
359 test = 0;
360 break;
361 } else {
362 test = 1;
363 }
364 } else {
365 if (Ck > precision) {
366 test = 0;
367 break;
368 } else { // Ck<=0
369 test = -1; // 1 peut etre
370 }
371 }
372 }
373
374 if (test == 1) {
375 R a0 = -pow((detXY / (Xi * X0)), 2);
376 R a1 = 2. * (pow(Ri / Xi, 2) + pow(R0 / X0, 2));
377 if (((a0 * bmax + a1) * bmax) < ((a0 * bmin + a1) * bmin)) {
378 bik = bmax;
379 } else {
380 bik = bmin;
381 }
382
383 aik =
384 (Ri * Ri * Y0 * X0 - R0 * R0 * Yi * Xi + bik * Yi * Y0 * detXY) / (detXY * Xi * X0);
385 cik = (-Ri * Ri * X0 * X0 + R0 * R0 * Xi * Xi - bik * (Yi * X0 + Y0 * Xi) * detXY) /
386 (detXY * Xi * X0);
387
388 assert((4. * aik * bik - cik * cik) >= 0.); // aire positive
389 assert(abs((4. * aik * bik - cik * cik) - pow(2. / (Rmax * (r0 - epsilon0)), 2)) >
390 0); // aire positive
391 if ((4. * aik * bik - cik * cik) <= (4. * A * B - C * C)) {
392 A = aik;
393 B = bik;
394 C = cik;
395 }
396 }
397
398 // -----------------------------------------------------------------
399 }
400 }
401 }
402 } else {
403 A = B = 1. / (Rmin * Rmin);
404 C = 0.;
405 }
406
407 delete[] PPoint;
408 }
409
410 // int LireTaille (const char *NomDuFichier, int &nbnoeuds) { // Lire le maillage sur le fichier
411 // de nom NomDuFichier
412 // // Ouverture
413 // du fichier a partir de son nom ifstream f(NomDuFichier); string buffer;
414 //
415 // nbnoeuds = 0;
416 // if (!f) {
417 // cerr << "Erreur a l'ouverture du fichier " << NomDuFichier << endl;
418 // return 1;
419 // }
420 //
421 // while (getline(f, buffer, '\n')) {
422 // if ((buffer[0] != '#') && (buffer != "")) {
423 // nbnoeuds += 1;
424 // }
425 // }
426 //
427 // return 0;
428 // }
429 //
430 // int Lire (const char *NomDuFichier, int n, R2 noeuds []) {
431 // ifstream f(NomDuFichier);
432 // string buffer;
433 // int i = 0;
434 //
435 // while (i < n) {
436 // f >> buffer;
437 // if (buffer[0] == '#')
438 // {getline(f, buffer);} else {
439 // // Lecture X Y Z de chacun des noeuds
440 //
441 // from_string(buffer, noeuds[i++].x);
442 // f >> noeuds[i - 1].y >> buffer;
443 // }
444 // }
445 //
446 // return 0;
447 // }
448
Min(const R a,const R b)449 R Min(const R a, const R b) { return a < b ? a : b; }
450
Max(const R a,const R b)451 R Max(const R a, const R b) { return a > b ? a : b; }
452
453 // metrixkuate(Th,np,o,err,[m11,m12,m22]);
454 class MetricKuate : public E_F0mps {
455 public:
456 typedef bool Result;
457 Expression expTh;
458 Expression expnp;
459 Expression exphmin;
460 Expression exphmax;
461 Expression experr;
462 Expression m11, m12, m22;
463 Expression px, py;
464
MetricKuate(const basicAC_F0 & args)465 MetricKuate(const basicAC_F0 &args) {
466 args.SetNameParam( );
467 expTh = to< pmesh >(args[0]); // a the expression to get the mesh
468 expnp = to< long >(args[1]); // a the expression to get the mesh
469 exphmin = to< double >(args[2]); // a the expression to get the mesh
470 exphmax = to< double >(args[3]); // a the expression to get the mesh
471 experr = to< double >(args[4]); // a the expression to get the mesh
472 // a array expression [ a, b]
473 const E_Array *ma = dynamic_cast< const E_Array * >((Expression)args[5]);
474 const E_Array *mp = dynamic_cast< const E_Array * >((Expression)args[6]);
475 if (ma->size( ) != 3) {
476 CompileError("syntax: MetricKuate(Th,np,o,err,[m11,m12,m22],[xx,yy])");
477 }
478
479 if (mp->size( ) != 2) {
480 CompileError("syntax: MetricKuate(Th,np,o,err,[m11,m12,m22],[xx,yy])");
481 }
482
483 m11 = CastTo< KN< double > * >((*ma)[0]); // fist exp of the array (must be a double)
484 m12 = CastTo< KN< double > * >((*ma)[1]); // second exp of the array (must be a double)
485 m22 = CastTo< KN< double > * >((*ma)[2]); // second exp of the array (must be a double)
486 px = CastTo< double * >((*mp)[0]); // fist exp of the array (must be a double)
487 py = CastTo< double * >((*mp)[1]); // second exp of the array (must be a double)
488 }
489
~MetricKuate()490 ~MetricKuate( ) {}
491
typeargs()492 static ArrayOfaType typeargs( ) {
493 return ArrayOfaType(atype< pmesh >( ), atype< long >( ), atype< double >( ), atype< double >( ),
494 atype< double >( ), atype< E_Array >( ), atype< E_Array >( ));
495 }
496
f(const basicAC_F0 & args)497 static E_F0 *f(const basicAC_F0 &args) { return new MetricKuate(args); }
498
499 AnyType operator( )(Stack s) const;
500 };
501
502 // the evaluation routine
operator ( )(Stack stack) const503 AnyType MetricKuate::operator( )(Stack stack) const {
504 MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
505 const Mesh *pTh = GetAny< pmesh >((*expTh)(stack));
506 long np = GetAny< long >((*expnp)(stack));
507 double hmin = GetAny< double >((*exphmin)(stack));
508 double hmax = GetAny< double >((*exphmax)(stack));
509
510 KN< double > *pm11, *pm12, *pm22;
511 double *pxx, *pyy;
512 pm11 = GetAny< KN< double > * >((*m11)(stack));
513 pm22 = GetAny< KN< double > * >((*m22)(stack));
514 pm12 = GetAny< KN< double > * >((*m12)(stack));
515 pxx = GetAny< double * >((*px)(stack));
516 pyy = GetAny< double * >((*py)(stack));
517 ffassert(pTh);
518 KN< R2 > Pt(np);
519 const Mesh &Th(*pTh);
520 cout << " MetricKuate " << np << " hmin = " << hmin << " hmax = " << hmax << " nv = " << Th.nv
521 << endl;
522
523 ffassert(pm11->N( ) == Th.nv);
524 ffassert(pm12->N( ) == Th.nv);
525 ffassert(pm22->N( ) == Th.nv);
526 {
527 for (int iv = 0; iv < Th.nv; iv++) {
528 R2 P = Th(iv);
529 MeshPointStack(stack)->set(P.x, P.y);
530 double m11 = 1, m12 = 0, m22 = 1;
531
532 for (int i = 0; i < np; i++) {
533 double t = (M_PI * 2. * i + 0.5) / np;
534 *pxx = cos(t);
535 *pyy = sin(t);
536 double ee = fabs(GetAny< double >((*experr)(stack)));
537 *pxx *= M_E;
538 *pyy *= M_E;
539 double eee = fabs(GetAny< double >((*experr)(stack)));
540 ee = max(ee, 1e-30);
541 eee = max(eee, 1e-30);
542 double p = Min(Max(log(eee) - log(ee), 0.1), 10);
543 double c = pow(1. / ee, 1. / p);
544 c = min(max(c, hmin), hmax);
545 Pt[i].x = *pxx * c / M_E;
546 Pt[i].y = *pyy * c / M_E;
547 if (iv == 0) {
548 cout << Pt[i] << " ++++ " << i << " " << t << " " << p
549 << " c = " << R2(*pxx * c / M_E, *pyy * c / M_E) << "e= " << ee << " " << eee << " "
550 << c << endl;
551 }
552 }
553
554 double epsilon = 1e-5;
555 metrique(np, Pt, m11, m22, m12, epsilon);
556 if (iv == 0) {
557 cout << " ---- 11,12,22 : " << m11 << " " << m12 / 2. << " " << m22 << endl;
558 }
559
560 (*pm11)[iv] = m11;
561 (*pm12)[iv] = m12 / 2.;
562 ;
563 (*pm22)[iv] = m22;
564 }
565 }
566 *mp = mps;
567 return true;
568 }
569
Load_Init()570 static void Load_Init( ) {
571 cout << "\n -- lood: init MetricKuate\n";
572 Global.Add("MetricKuate", "(", new OneOperatorCode< MetricKuate >( ));
573 }
574
575 LOADFUNC(Load_Init)
576