1 // used by splitsimplex.cpp
2 // ---- build a the simplex decomposition
3 /*
4 template void SplitSimplex<R1>(int N,int & nv, R1 *& P, int & nk , int *& K);
5 template void SplitSimplex<R2>(int N,int & nv, R2 *& P, int & nk , int *& K);
6 template void SplitSimplex<R3>(int N,int & nv, R3 *& P, int & nk , int *& K);
7 $:
8 */
9 //
10 //
11 // ORIG-DATE: fev 2009
12 // -*- Mode : c++ -*-
13 //
14 // SUMMARY : Model mesh 2d
15 // USAGE : LGPL
16 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
17 // AUTHOR : Frederic Hecht
18 // E-MAIL : frederic.hecht@ann.jussieu.fr
19 //
20
21 /*
22
23 This file is part of Freefem++
24
25 Freefem++ is free software; you can redistribute it and/or modify
26 it under the terms of the GNU Lesser General Public License as published by
27 the Free Software Foundation; either version 2.1 of the License, or
28 (at your option) any later version.
29
30 Freefem++ is distributed in the hope that it will be useful,
31 but WITHOUT ANY WARRANTY; without even the implied warranty of
32 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 GNU Lesser General Public License for more details.
34
35 You should have received a copy of the GNU Lesser General Public License
36 along with Freefem++; if not, write to the Free Software
37 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
38
39 Thank to the ARN () FF2A3 grant
40 ref:ANR-07-CIS7-002-01 */
41
42 #include <cmath>
43 #include <cstdlib>
44 #include <iostream>
45 #include <cassert>
46 #include <fstream>
47 #include "ufunction.hpp"
48 using namespace std;
49 namespace Fem2D {
50 #include "R3.hpp"
51 };
52 using Fem2D::R1;
53 using Fem2D::R2;
54 using Fem2D::R3;
55 extern long verbosity ;
56
57 #include "splitsimplex.hpp"
58
59 /*
60 construction an array of sub simplex for plot ...
61 see last template function
62 SplitSimplex(N,nv,P,nk,K);
63 N > 1 -> classical split un N^d sub simplex
64 nv : number of point on
65 */
66 // d = 1 trivial
SplitSimplex(int N,R1 * P,int * K,int op=0,R1 * AB=0)67 void SplitSimplex(int N,R1 *P,int *K,int op=0,R1 *AB=0)
68 {
69 assert(N>0);
70 double h = 1./ N;
71 for(int i=0;i<=N;++i)
72 if(AB)
73 P[i+op] = R1(i*h).Bary(AB);
74 else
75 P[i+op] = R1(i*h);
76 int l=0;
77 for(int i=0;i<N;++i) {
78 K[l++]=i+op;
79 K[l++]=i+1+op;
80 if(verbosity>200)
81 cout << "l="<< l/2 << " "<< K[l-2] <<" "<< K[l-1] <<" "<< endl;
82 }
83 }
84 /*
85 0124 1234 1345 3456 2346 3567 1 solutions de taille 6 trouvees.
86 0124 1245 1235 2356 2456 3567 Solution courte numero 7:
87 0124 1246 1236 1356 1456 3567 Solution courte numero 12:
88 */
89
90
91
92 typedef int int4 [4] ;
93
94 // InvIntFunc
95
NumSimplex2(int i)96 inline int NumSimplex2(int i) {return ((i)*(i+1))/2;}
NumSimplex2_1(int l)97 inline int NumSimplex2_1(int l) { return sqrt(l*2)+3;}
NumSimplex3(int i)98 inline int NumSimplex3(int i) {return ((i)*(i+1)*(i+2))/6;}
NumSimplex3_1(int l)99 inline int NumSimplex3_1(int l) { return pow(l*6,1./3)+4;}
100
101 #define InvIntFunction invNumSimplex2
102 #define F(i) NumSimplex2(i)
103 #define F_1(i) NumSimplex2_1(i)
104 #include "InvIntFunc.cpp"
105 #undef F
106 #undef F_1
107 #undef InvIntFunction
108
109 #define InvIntFunction invNumSimplex3
110 #define F(i) NumSimplex3(i)
111 #define F_1(i) NumSimplex3_1(i)
112 #include "InvIntFunc.cpp"
113 #undef F
114 #undef F_1
115 #undef InvIntFunction
116
NumSimplex1(int i)117 inline int NumSimplex1(int i) { return i;}
NumSimplex2(int i,int j)118 inline int NumSimplex2(int i,int j) { return j+NumSimplex2(i+j);}
NumSimplex3(int i,int j,int k)119 inline int NumSimplex3(int i,int j,int k) { return NumSimplex3(i+j+k)+NumSimplex2(j+k)+k;}
120
121
invNumSimplex2(int n,int & i,int & j)122 inline void invNumSimplex2(int n,int &i,int &j)
123 {
124 int k= invNumSimplex2(n); //( i+j)
125 j=n-NumSimplex2(k);
126 i= k-j;
127 // cout << n << " " << k << " -> " << i << " " << j << endl;
128 assert( n == NumSimplex2(i,j));
129 }
invNumSimplex3(int n,int & i,int & j,int & k)130 inline void invNumSimplex3(int n,int &i,int &j,int &k)
131 {
132 int l= invNumSimplex3(n); //=( i+j+k)
133 invNumSimplex2(n-NumSimplex3(l),j,k);
134 assert(j>=0 && k>=0);
135 i=l-k-j;
136 // cout << n << " " << l << "-> " << i << " " << j << " " << k <<endl;
137 assert( n == NumSimplex3(i,j,k)) ;
138 }
139
140
141 // d = 2
SplitSimplex(int N,R2 * P,int * K,int op=0,R2 * ABC=0)142 void SplitSimplex(int N,R2 *P,int *K,int op=0,R2 *ABC=0)
143 {
144 // warning PB of oriention if (i+j<N) , the transo is positif => no swap on 1, 2 ve
145 assert(N>0);
146 int nv = (N+1)*(N+2)/2;
147 double h=1./N;
148 // loop sur les diag i+j = k
149 // num ( i+j,j) lexico croissant
150 for(int l=0;l<nv;l++)
151 {
152 int i,j;
153 invNumSimplex2(l,i,j);
154 if(ABC)
155 P[l+op]= R2(i*h,j*h).Bary(ABC);
156 else
157 P[l+op]= R2(i*h,j*h);
158 assert(l<nv);
159 }
160 // generation des trianges
161 // --------
162 int l=0;
163 for (int i=0;i<N;++i)
164 for (int j=0;j<N;++j)
165 if(i+j<N)
166 {
167 K[l++]= op+NumSimplex2(i,j);
168 K[l++]= op+NumSimplex2(i+1,j);
169 K[l++]= op+NumSimplex2(i,j+1);
170 }
171 else
172 {
173 K[l++]= op+NumSimplex2(N-i,N-j);
174 K[l++]= op+NumSimplex2(N-i-1,N-j);// swap vertices 1 and 2 May 2019 FH (Thanks to AF)
175 K[l++]= op+NumSimplex2(N-i,N-j-1);
176 }
177 }
178
179 /*
180 // d = 3 Surfacic
181 void SplitSimplex(int N,R3 *P,int *K,int op=0,R3 *ABC=0)
182 {
183 assert(N>0);
184 int nv = (N+1)*(N+2)/2;
185 double h=1./N;
186 // loop sur les diag i+j = k
187 // num ( i+j,j) lexico croissant
188 for(int l=0;l<nv;l++)
189 {
190 int i,j,k;
191 invNumSimplex3(l,i,j,k);
192 if(ABC)
193 P[l+op]= R3(i*h,j*h,k*h).Bary(ABC);
194 else
195 P[l+op]= R3(i*h,j*h,k*h);
196 assert(l<nv);
197 }
198 // generation des trianges
199 // --------
200 int l=0;
201 for (int i=0;i<N;++i)
202 for (int j=0;j<N;++j)
203 if(i+j<N)
204 {
205 K[l++]= op+NumSimplex2(i,j);
206 K[l++]= op+NumSimplex2(i+1,j);
207 K[l++]= op+NumSimplex2(i,j+1);
208 }
209 else
210 {
211 K[l++]= op+NumSimplex2(N-i,N-j);
212 K[l++]= op+NumSimplex2(N-i,N-j-1);
213 K[l++]= op+NumSimplex2(N-i-1,N-j);
214 }
215 } */
216
217
SplitSimplex(int N,R3 * P,int * tet,int op=0,R3 * Khat=0)218 void SplitSimplex(int N,R3 *P,int *tet,int op=0,R3* Khat=0)
219 {
220 const int n=N;
221 const int n2=n*n;
222 const int n3=n2*n;
223 const int ntc=6;
224 int nv = (N+1)*(N+2)*(N+3)/6;
225
226 int d1[6][4] = { {0,1,2,4} , {1,2,4,3},{1,3,4,5},{3,4,5,6},{ 2,3,6,4}, {3,5,7,6} };
227 int ntt=0;
228 int n8[8];
229 int ptet=0;
230 double h=1./N;
231 for(int l=0;l<nv;l++)
232 {
233 int i,j,k;
234 invNumSimplex3(l,i,j,k);
235 assert( i>=0&& j >=0 && k >= 0);
236 assert( i<=N && j <= N && k <= N);
237 if(Khat)
238 P[l+op]= R3(i*h,j*h,k*h).Bary(Khat);
239 else
240 P[l+op]= R3(i*h,j*h,k*h);
241 assert(l<nv);
242 }
243 // comment n3 i=m[n]
244 // m = i+j*n+k*n2
245 for (int m=0;m<n3;m++)
246 {
247 int i = m % n;
248 int j = (m / n) % n;
249 int k = m /( n2);
250 for (int l=0;l<8;++l)
251 {
252 int ii = ( i + ( (l & 1) != 0) );
253 int jj = ( j + ( (l & 2) != 0) );
254 int kk = ( k + ( (l & 4) != 0) );
255 int ll= NumSimplex3(ii,jj,kk);
256 n8[l]=ll;
257 if(ii+jj+kk>n) n8[l]=-1;
258 }
259 for (int l=0;l<ntc;++l)
260 if(ntt<n3)
261 {
262 int out =0;
263 for(int m=0;m<4;++m)
264 if( (tet[ptet++]= n8[d1[l][m]]+op) < op) out++;
265 if(out == 0 ) ntt++;
266 else ptet -= 4; // remove tet
267 }
268
269 }
270 /* if(verbosity>199)
271 {
272
273 cout << " SplitSimplex " << endl;
274 for (int i=0,l=0;i<n3;i++)
275 for(int m=0;m<4;++m)
276 cout << tet[l++] << (m==3 ? '\n' : ' ' );
277 cout << ptet << " " << tet << endl;
278 }*/
279 assert(ntt==n3);
280 }
281
282 // Add J. Morice (trunc functions)
283
SplitSurfaceSimplex(int N,int & ntri2,int * & tri)284 void SplitSurfaceSimplex(int N,int &ntri2,int *&tri)
285 {
286 const int n=N;
287 const int n2=n*n;
288
289 int ntri=3*ntri2;
290 int op=0;
291
292 tri = new int[ntri];
293 // generation des trianges
294 // --------
295
296 // face i=0
297 int l=0;
298 if(verbosity>200)
299 cout << "face i=0" << endl;
300 for (int i=0;i<N;++i)
301 for (int j=0;j<N;++j){
302 if(i+j<N)
303 {
304 tri[l++]= op+NumSimplex3(0,i,j);
305 tri[l++]= op+NumSimplex3(0,i+1,j);
306 tri[l++]= op+NumSimplex3(0,i,j+1);
307 }
308 else
309 { // correction orinatation ...
310 tri[l++]= op+NumSimplex3(0,N-i,N-j);
311 tri[l++]= op+NumSimplex3(0,N-i-1,N-j);
312 tri[l++]= op+NumSimplex3(0,N-i,N-j-1);
313 }
314 //cout << "i,j " << i << "," << j << endl;
315 if(verbosity>200)
316 cout << "l="<< l/3 << " "<< tri[l-3] <<" "<< tri[l-2] <<" "<< tri[l-1] <<" "<< endl;
317 }
318 // face j=0
319 if(verbosity>200)
320 cout << "face j=0" << endl;
321 for (int i=0;i<N;++i)
322 for (int j=0;j<N;++j){
323 if(i+j<N)
324 {
325 tri[l++]= op+NumSimplex3(i,0,j);
326 tri[l++]= op+NumSimplex3(i,0,j+1); // inverser les deux lignes
327 tri[l++]= op+NumSimplex3(i+1,0,j);
328 }
329 else
330 { // CHANGE ORIENTATION FH MAY
331 tri[l++]= op+NumSimplex3(N-i,0,N-j);
332 tri[l++]= op+NumSimplex3(N-i,0,N-j-1);
333 tri[l++]= op+NumSimplex3(N-i-1,0,N-j); // inverser les deux lignes
334 }
335 //cout << "i,j " << i << "," << j << endl;
336 if(verbosity>200)
337 cout << "l="<< l/3 << " "<< tri[l-3] <<" "<< tri[l-2] <<" "<< tri[l-1] <<" "<< endl;
338 }
339 // face k=0
340 if(verbosity>200)
341 cout << "face k=0" << endl;
342 for (int i=0;i<N;++i)
343 for (int j=0;j<N;++j){
344 if(i+j<N)
345 {
346 tri[l++]= op+NumSimplex3(i,j,0);
347 tri[l++]= op+NumSimplex3(i+1,j,0);
348 tri[l++]= op+NumSimplex3(i,j+1,0);
349 //tri[l++]= op+NumSimplex3(i+1,j,0);
350 }
351 else
352 { // CHANGE ORIENTATION FH MAY
353 tri[l++]= op+NumSimplex3(N-i,N-j,0);
354 tri[l++]= op+NumSimplex3(N-i-1,N-j,0);
355 tri[l++]= op+NumSimplex3(N-i,N-j-1,0);
356 //tri[l++]= op+NumSimplex3(N-i,N-j-1,0);
357 }
358 //cout << "i,j " << i << "," << j << endl;
359 if(verbosity>200)
360 cout << "l="<< l/3 << " "<< tri[l-3] <<" "<< tri[l-2] <<" "<< tri[l-1] <<" "<< endl;
361 }
362 // face i+j+k=1
363 if(verbosity>200)
364 cout << "dernier face " << endl;
365 for (int k=0;k<N;++k)
366 for (int j=0;j<N;++j){
367 if(k+j<N)
368 {
369 int i=N-j-k;
370 tri[l++]= op+NumSimplex3( i, j, k);
371 tri[l++]= op+NumSimplex3( i-1, j, k+1);
372 tri[l++]= op+NumSimplex3( i-1, j+1, k);
373 }
374 else
375 {
376 int i=N-(N-j-1)-(N-k);
377 tri[l++]= op+NumSimplex3( i, N-j-1, N-k);
378 tri[l++]= op+NumSimplex3( i-1, N-j, N-k);
379 tri[l++]= op+NumSimplex3( i, N-j, N-k-1);
380 }
381 if(verbosity>200)
382 cout << "l="<< l/3 << " "<< tri[l-3] <<" "<< tri[l-2] <<" "<< tri[l-1] <<" "<< endl;
383 }
384 if(verbosity>200)
385 cout << "l= " << l << " ntri=" << ntri << endl;
386 assert( l == ntri);
387 }
388
389
390
391
SplitEdgeSimplex(int N,int & nedge2,int * & edge)392 void SplitEdgeSimplex(int N,int &nedge2,int *&edge)
393 {
394 const int n=N;
395 const int n2=n*n;
396
397 int nedge=2*nedge2;
398 int op=0;
399
400 edge = new int[nedge];
401 // generation des edges
402 // --------
403 // (i,j,k) barycentric coordinates
404
405 int l=0;
406
407 for (int i=0;i<N;++i)
408 for (int j=0;j<N;++j) {
409 if(i+j<N) {
410 edge[l++]= op+NumSimplex2(i+1,j);
411 edge[l++]= op+NumSimplex2(i,j+1);
412 }
413 if(verbosity>200)
414 cout << "l="<< l/2 <<" "<< edge[l-2] <<" "<< edge[l-1] <<" "<< endl;
415 }
416
417 for (int i=0;i<N;++i)
418 for (int j=0;j<N;++j) {
419 if(i+j<N) {
420 edge[l++]= op+NumSimplex2(i,j);
421 edge[l++]= op+NumSimplex2(i,j+1);
422 }
423 if(verbosity>200)
424 cout << "l="<< l/2 <<" "<< edge[l-2] <<" "<< edge[l-1] <<" "<< endl;
425 }
426
427 for (int i=0;i<N;++i)
428 for (int j=0;j<N;++j) {
429 if(i+j<N) {
430 edge[l++]= op+NumSimplex2(i,j);
431 edge[l++]= op+NumSimplex2(i+1,j);
432 }
433 if(verbosity>200)
434 cout << i+j << "l="<< l/2 <<" "<< edge[l-2] <<" "<< edge[l-1] <<" "<< endl;
435 }
436
437 if(verbosity>200)
438 cout << "l= " << l << " nedge=" << nedge << endl;
439 assert( l == nedge);
440 }
441
442
443
444 /*
445 void SplitSimplex(int N,int & nv, R1 *& P, int & nk , int *& K)
446 {
447 typedef R1 Rd;
448 const int d = Rd::d;
449 int cas = (N>0) ? 1 : d+1;
450 N=abs(N);
451 assert(N);
452 int nv1=(N+1);
453 int nk1= N;
454 nv = cas*nv1;
455 nk = cas*nk1;
456
457 P = new Rd[nv];
458 K = new int [nk*(d+1)];
459 if( cas ==1)
460 SplitSimplex( N, P,K) ;
461 else
462 {
463 Rd AB1[2]= { Rd(0),Rd(0.5)};
464 SplitSimplex( N, P,K,0,AB1) ;
465 Rd AB2[2]= { Rd(0.5),Rd(1)};
466 SplitSimplex( N, P,K+nk1,nv1,AB2);
467 }
468 }
469
470
471 void SplitSimplex(int N,int & nv, R2 *& P, int & nk , int *& K)
472 {
473 typedef R2 Rd;
474 const int d = Rd::d;
475 int cas = (N>0) ? 1 : d+1;
476 assert(N);
477 N=abs(N);
478 int nv1=N*(N+1)/2;
479 int nk1= N*N;
480 nv = cas*nv1;
481 nk = cas*nk1;
482 P = new Rd[nv];
483 K = new int [nk*(d+1)];
484 if( cas ==1)
485 SplitSimplex( N, P,K);
486 else
487 {
488 Rd G=Rd::diag(1./(d+1));
489 R2 Khat[d+1];
490 for (int i=0;i<d;++i)
491 Khat[i+1][i]=1; // modit 25/2/2009
492 for(int i=0;i<=d;++i)
493 {
494 Rd S=Khat[i];
495 Khat[i]=G;
496 SplitSimplex( N, P,K+nk1*i,nv1*i,Khat);
497 Khat[i]=S;
498 }
499 }
500 }
501
502 */
503
504 template<class Rd>
SplitSimplex(int N,int & nv,Rd * & P,int & nk,int * & K)505 void SplitSimplex(int N,int & nv, Rd *& P, int & nk , int *& K)
506 {
507 const int d = Rd::d,d1 =d+1;
508 int cas = (N>0) ? 1 : d+1;
509 assert(N);
510 N=abs(N);
511 int nv1=N+1; // nb simplexe
512 int nk1=N; // nb vertices
513 for(int i=2;i<=d;++i)
514 {
515 nk1 *=N;
516 nv1 = (nv1)*(N+i)/i;
517 }
518 nv = cas*nv1;
519 nk = cas*nk1;
520 P = new Rd[nv]; // no bug correct jan 2024 FH (thank to OP)
521 K = new int[nk*d1]; // no bug correct jan 2024 FH (thank to OP)
522 if( cas ==1)
523 SplitSimplex( N, P,K);
524 else
525 {
526 Rd G=Rd::diag(1./(d1));
527 for(int i=0;i<=d;++i)
528 {
529 Rd Khat[d+1];
530 for (int j=1;j<=d;++j)
531 Khat[j][j-1]=1;// bug correct jan 2024 FH (thank to OP)
532 Khat[i]=G;
533 SplitSimplex( N, P,K+(nk1*d1)*i,nv1*i,Khat); // FH no recursion here ...
534 }
535 }
536 if(verbosity>99)
537 {
538 cout << "SplitSimplex : nv =" << nv << " nk :" << nk << " " << N << endl ;
539 for(int i=0; i< nv ; ++i)
540 cout << i << " / " << P[i] << endl;
541
542 for(int k=0,kk=0; k < nk; ++k)
543 {
544 cout << k << " " << kk << " : ";
545 for(int m=0;m<d+1;++m)
546 cout << K[kk++] << " " ;
547 cout << endl;
548 }
549
550 }
551 }
552
553 template void SplitSimplex<R1>(int N,int & nv, R1 *& P, int & nk , int *& K);
554 template void SplitSimplex<R2>(int N,int & nv, R2 *& P, int & nk , int *& K);
555 template void SplitSimplex<R3>(int N,int & nv, R3 *& P, int & nk , int *& K);
556 /*
557 int main(int argc,const char ** argv)
558 {
559 R3 *P;
560 int *K,nv,nk;
561 int N=2;
562 if(argc>1) N=atoi(argv[1]);
563 SplitSimplex(N,nv,P,nk,K);
564 cout << P << " " << K << endl;
565 cout << N << " nv " << nv << " nk =" << nk << endl;
566 for(int i=0;i<nv;++i) cout << P[i] << endl;
567 for(int i=0,l=0;i<nk;i++)
568 {
569 for(int j=0;j<4;j++)
570 cout << K[l++] << " ";
571 cout << endl;
572 }
573
574
575 }
576 */
577