forked from rmerks/VirtualLeaf2021
-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathnCercles.py
More file actions
1347 lines (1070 loc) · 54.9 KB
/
Copy pathnCercles.py
File metadata and controls
1347 lines (1070 loc) · 54.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import math as math
import matplotlib.pyplot as plt
from random import *
# ----------------------------------------------------------------------------------------------------------#
# ------------- CONVERSION POLAIRE-CARTÉSIEN ---------------------------------------------------------------#
# ----------------------------------------------------------------------------------------------------------#
def polaire_vers_cartesien(k_noeuds_n_cercle_polaire):
"""Convertit les coordonnées polaires des noeuds en coordonnées cartésiennes.
Args:
k_noeuds_n_cercle_polaire: Liste des noeuds en coordonnées polaires pour chaque cercle
Returns:
Liste des noeuds en coordonnées cartésiennes pour chaque cercle
"""
return [
{
"x": noeud["r"] * math.cos(noeud["theta"]),
"y": noeud["r"] * math.sin(noeud["theta"]),
"index": noeud["index"],
"initial": noeud["initial"],
"sam": noeud["sam"],
"boundary": noeud["boundary"],
"fixed": noeud["fixed"]
}
for cercle in k_noeuds_n_cercle_polaire
for noeud in cercle
]
#-----------------------------------------------------------------------------------------------------#
#------------ GÉNÉRATION DES RAYONS ------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------#
def n_cercles(n, rayon_0, a, rapport_R_T, debug=False):
"""Calcule les rayons de n cercles concentriques.
Args:
n: Nombre de cercles à générer
rayon_0: Rayon du premier cercle (le plus interne)
a: Liste des nombres de points/cellules par cercle
rapport_R_T: Liste des rapports entre rayon et longueur de corde
debug: Active/désactive l'affichage des informations de débogage
Returns:
Liste des rayons des cercles générés
"""
if n <= 0:
return []
rayon_n_cercles = [rayon_0] # Initialisation avec le premier rayon
if debug:
print(f"cercle 0: rayon = {rayon_0}\n")
for i in range(1, n):
rayon_precedent = rayon_n_cercles[i - 1]
# Sélection de l'indice correct pour a
a_idx = i - 1 if i < n - 1 or len(a) <= i - 1 else i - 2
# Facteur multiplicatif pour calculer l'incrément de rayon
# Calcul du nouveau rayon
sin_term = math.sin(math.pi / a[a_idx])
increment = rapport_R_T[i - 1] * 2 * sin_term * rayon_precedent
nouveau_rayon = rayon_precedent + increment
rayon_n_cercles.append(nouveau_rayon)
if debug:
print(f"cercle {i}:")
print(f"rayon précédent = {rayon_precedent}")
print(f"nouveau rayon = {nouveau_rayon}")
print(f"incrément = {increment}")
if i < n - 1:
ratio_calcule = (nouveau_rayon - rayon_precedent) / (2 * sin_term * rayon_precedent)
print(f"ratio R/T calculé = {ratio_calcule}")
print(f"ratio R/T attendu = {rapport_R_T[i]}")
else:
ratio_calcule = (nouveau_rayon - rayon_precedent) / (sin_term * rayon_precedent)
print(f"ratio R/T = {ratio_calcule}")
print(f"a[{a_idx}] = {a[a_idx]}\n")
return rayon_n_cercles
#----------------------------------------------------------------------------------------------------#
#---------- CRÉATION DES NOEUDS INITIAUX ------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
def creer_noeuds_initiaux(n, rayon_n_cercles, a):
"""Crée les nœuds initiaux sans fusion.
Args:
n: Nombre de cercles
rayon_n_cercles: Liste des rayons de chaque cercle
a: Liste des nombres de points/cellules par cercle
Returns:
Tuple contenant:
- Liste des noeuds en coordonnées polaires pour chaque cercle
- Nombre total de noeuds créés
"""
k_noeuds_n_cercle_polaire = []
index_noeud = 0
# Fonction auxiliaire pour créer un nœud
def creer_noeud(r, theta, boundary=False, fixed=False):
nonlocal index_noeud
noeud = {
"r": r,
"theta": theta,
"index": index_noeud,
"initial": True,
"sam": False,
"boundary": boundary,
"fixed": fixed
}
index_noeud += 1
return noeud
for i in range(n):
rayon = rayon_n_cercles[i]
cercle_noeuds = []
if i == 0:
# Premier cercle - tous les nœuds sont fixes
cercle_noeuds = [
creer_noeud(
r=rayon,
theta=2 * math.pi * j / a[i],
fixed=True
) for j in range(a[i])
]
elif 0 < i < n - 1:
# Cercles intermédiaires
# D'abord les points alignés avec le cercle précédent
cercle_noeuds.extend([
creer_noeud(
r=rayon,
theta=2 * math.pi * j / a[i-1]
) for j in range(a[i-1])
])
# Puis les points propres à ce cercle
cercle_noeuds.extend([
creer_noeud(
r=rayon,
theta=2 * math.pi * j / a[i]
) for j in range(a[i])
])
else:
# Dernier cercle - tous les points sont sur la frontière
cercle_noeuds = [
creer_noeud(
r=rayon,
theta=2 * math.pi * j / a[i-1],
boundary=True
) for j in range(a[i-1])
]
k_noeuds_n_cercle_polaire.append(cercle_noeuds)
return k_noeuds_n_cercle_polaire, index_noeud
#----------------------------------------------------------------------------------------------------#
#------------- FUSION DES NOEUDS PROCHES ------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
def fusionner_noeuds(k_noeuds_n_cercle_polaire, n, index_noeud_depart=0, seuil_distance=5):
""" Fusionne les nœuds qui se chevauchent.
Args:
k_noeuds_n_cercle_polaire: Liste des noeuds en coordonnées polaires
n: Nombre de cercles
index_noeud_depart: Index de départ pour les nouveaux noeuds créés
seuil_distance: Distance minimale entre deux noeuds pour les considérer distincts
Returns:
Tuple contenant:
- Liste mise à jour des noeuds en coordonnées polaires
- Dictionnaire de correspondance des indices fusionnés """
index_noeud = index_noeud_depart
# Trier chaque sous-liste par ordre croissant de theta
for i in range(len(k_noeuds_n_cercle_polaire)):
k_noeuds_n_cercle_polaire[i].sort(key=lambda noeud: noeud["theta"])
# Collecter tous les nœuds avec leur cercle d'origine
tous_noeuds = []
for i, cercle in enumerate(k_noeuds_n_cercle_polaire):
for noeud in cercle:
r = noeud["r"]
theta = noeud["theta"]
idx = noeud["index"]
initial = noeud["initial"]
sam = noeud["sam"]
boundary = noeud["boundary"]
fixed = noeud["fixed"]
x = r * math.cos(theta)
y = r * math.sin(theta)
tous_noeuds.append((x, y, r, theta, idx, i, initial, sam,boundary,fixed))
# Identifier les paires de nœuds qui se chevauchent
map_indices = {}
noeuds_a_fusionner = []
for i in range(len(tous_noeuds)):
x1, y1, r1, theta1, idx1, cercle1, initial1,sam1,boundary1,fixed1 = tous_noeuds[i]
for j in range(i + 1, len(tous_noeuds)):
x2, y2, r2, theta2, idx2, cercle2, initial2,sam2,boundary2,fixed2 = tous_noeuds[j]
# Distance entre les nœuds
distance = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
if distance < seuil_distance:
print(f"Nœuds à fusionner : {idx1} et {idx2} (distance={distance:.3f})")
noeuds_a_fusionner.append((i, j))
# Création des nouveaux nœuds
nouveaux_noeuds = []
indices_a_supprimer = set()
for i, j in noeuds_a_fusionner:
x1, y1, r1, theta1, idx1, cercle1, initial1,sam1,boundary1,fixed1 = tous_noeuds[i]
x2, y2, r2, theta2, idx2, cercle2, initial2,sam2,boundary2,fixed2 = tous_noeuds[j]
# Créer le nœud entre les deux
x_nouveau = (x1 + x2) / 2
y_nouveau = (y1 + y2) / 2
r_nouveau = math.sqrt(x_nouveau ** 2 + y_nouveau ** 2)
theta_nouveau = math.atan2(y_nouveau, x_nouveau)
if theta_nouveau < 0:
theta_nouveau += 2 * math.pi
nouvel_indice = index_noeud
index_noeud += 1
# Un nœud est initial seulement si tous les nœuds fusionnés sont initiaux
initial_nouveau = initial1 and initial2
boundary_nouveau = boundary1 or boundary2
sam_nouveau = False
fixed_nouveau = False
# Création du nouveau nœud
nouveau_noeud = {
"r": r_nouveau,
"theta": theta_nouveau,
"index": nouvel_indice,
"initial": initial_nouveau,
"sam":sam_nouveau,
"boundary":boundary_nouveau,
"fixed": fixed_nouveau
}
nouveaux_noeuds.append((nouveau_noeud, cercle1, cercle2))
map_indices[idx1] = nouvel_indice
map_indices[idx2] = nouvel_indice
indices_a_supprimer.add(i)
indices_a_supprimer.add(j)
# Reconstruction des nœuds sans ceux qui sont supprimés
nouveau_k_noeuds_n_cercle_polaire = [[] for _ in range(n)]
# Ajouter les nœuds non supprimés
for i, (x, y, r, theta, idx, cercle, initial, sam,boundary, fixed) in enumerate(tous_noeuds):
if i not in indices_a_supprimer:
nouveau_k_noeuds_n_cercle_polaire[cercle].append({
"r": r,
"theta": theta,
"index": idx,
"initial": initial,
"sam": sam,
"boundary": boundary,
"fixed": fixed
})
# Ajouter les nouveaux nœuds
for nouveau_noeud, cercle1, cercle2 in nouveaux_noeuds:
nouveau_k_noeuds_n_cercle_polaire[cercle1].append(nouveau_noeud.copy())
if cercle1 != cercle2:
nouveau_k_noeuds_n_cercle_polaire[cercle2].append(nouveau_noeud.copy())
# Trier à nouveau
for i in range(len(nouveau_k_noeuds_n_cercle_polaire)):
nouveau_k_noeuds_n_cercle_polaire[i].sort(key=lambda noeud: noeud["theta"])
print(f"Nombre de nouveaux nœuds créés: {len(nouveaux_noeuds)}")
return nouveau_k_noeuds_n_cercle_polaire, map_indices
# ----------------------------------------------------------------------------------------------------#
# ------------ CRÉATION DES CELLULES AVEC NOEUDS INTERMÉDIAIRES ---------------------------------------#
# ----------------------------------------------------------------------------------------------------#
def creer_cellules(rayon_n_cercles, k_noeuds_n_cercle_polaire, n, a):
"""Crée des cellules à 4 nœuds entre cercles consécutifs.
Args:
rayon_n_cercles: Liste des rayons des cercles
k_noeuds_n_cercle_polaire: Coordonnées polaires des nœuds pour chaque cercle
n: Nombre de cercles
a: Liste des nombres de cellules entre chaque paire de cercles consécutifs
Returns:
Liste des cellules, chaque cellule étant définie par une liste de noeuds (indices)
"""
# Fonction auxiliaire pour la comparaison d'angles
def angles_equivalents(a1, a2, tolerance=0.001):
a1_norm = a1 % (2 * math.pi)
a2_norm = a2 % (2 * math.pi)
diff = min(abs(a1_norm - a2_norm), 2 * math.pi - abs(a1_norm - a2_norm))
return diff < tolerance
# Fonction pour vérifier si un angle est entre deux autres
def est_entre_angles(theta, debut, fin, tolerance=0.001):
theta_norm = theta % (2 * math.pi)
debut_norm = debut % (2 * math.pi)
fin_norm = fin % (2 * math.pi)
if debut_norm <= fin_norm:
return debut_norm - tolerance <= theta_norm <= fin_norm + tolerance
else:
return theta_norm >= debut_norm - tolerance or theta_norm <= fin_norm + tolerance
# Prétraitement : indexer et trier les nœuds par angle pour chaque cercle
noeuds_par_cercle = []
for i in range(n):
if i < len(k_noeuds_n_cercle_polaire):
noeuds_par_cercle.append(sorted(k_noeuds_n_cercle_polaire[i], key=lambda noeud: noeud["theta"]))
else:
noeuds_par_cercle.append([])
cellules_temp = [] # Liste temporaire pour les cellules
# Création de la cellule centrale (cercle 0)
if n > 0 and len(k_noeuds_n_cercle_polaire) > 0:
noeuds_cercle_0 = [noeud["index"] for noeud in noeuds_par_cercle[0]]
cellule_centrale = {
"index": 0,
"noeuds": noeuds_cercle_0,
"boundary": 0,
"cell_type": 1,
"target_area": 100.0,
"lambda_celllength": 0,
"at_boundary": True,
"dead": False,
"target_length": 0,
"stiffness": 1,
"source": False,
"pin_fixed": False,
"area": 0.0,
"fixed": False,
"div_counter": 0,
"cercle": 0,
"angle": 0
}
cellules_temp.append(cellule_centrale)
print(f"Cellule centrale créée avec {len(noeuds_cercle_0)} noeuds")
# Création des cellules entre les cercles
for i in range(n - 1):
nb_cellules = a[i] if i < len(a) else a[-1]
for j in range(nb_cellules):
angle_j = 2 * j * math.pi / nb_cellules
angle_j_moins_1 = 2 * (j - 1) * math.pi / nb_cellules
if angle_j_moins_1 < 0:
angle_j_moins_1 += 2 * math.pi
# Variables pour stocker les informations des nœuds trouvés
noeud1, noeud2, noeud3, noeud4 = None, None, None, None
angle1, angle2, angle3, angle4 = None, None, None, None
r1, r2, r3, r4 = None, None, None, None
tolerance = 0.001
# Recherche optimisée des nœuds aux angles spécifiques
for point in noeuds_par_cercle[i]:
theta = point["theta"]
if noeud1 is None and angles_equivalents(theta, angle_j, tolerance):
noeud1, angle1, r1 = point["index"], theta, point["r"]
elif noeud2 is None and angles_equivalents(theta, angle_j_moins_1, tolerance):
noeud2, angle2, r2 = point["index"], theta, point["r"]
# Si les deux nœuds sont trouvés, on peut passer au cercle suivant
if noeud1 is not None and noeud2 is not None:
break
for point in noeuds_par_cercle[i + 1]:
theta = point["theta"]
if noeud3 is None and angles_equivalents(theta, angle_j_moins_1, tolerance):
noeud3, angle3, r3 = point["index"], theta, point["r"]
elif noeud4 is None and angles_equivalents(theta, angle_j, tolerance):
noeud4, angle4, r4 = point["index"], theta, point["r"]
if noeud3 is not None and noeud4 is not None:
break
# Si tous les nœuds principaux sont trouvés
if noeud1 is not None and noeud2 is not None and noeud3 is not None and noeud4 is not None:
# Calculer le centre de la cellule
x1, y1 = r1 * math.cos(angle1), r1 * math.sin(angle1)
x2, y2 = r2 * math.cos(angle2), r2 * math.sin(angle2)
x3, y3 = r3 * math.cos(angle3), r3 * math.sin(angle3)
x4, y4 = r4 * math.cos(angle4), r4 * math.sin(angle4)
centre_x = (x1 + x2 + x3 + x4) / 4
centre_y = (y1 + y2 + y3 + y4) / 4
angle_centre = math.atan2(centre_y, centre_x)
if angle_centre < 0:
angle_centre += 2 * math.pi
# Recherche des nœuds intermédiaires
noeuds_arc_sup = []
for point in noeuds_par_cercle[i]:
theta, index = point["theta"], point["index"]
if index != noeud1 and index != noeud2:
if (est_entre_angles(theta, angle2, angle1, tolerance) or
angles_equivalents(theta, angle1, tolerance) or
angles_equivalents(theta, angle2, tolerance)):
noeuds_arc_sup.append((theta, index))
noeuds_arc_sup.sort(reverse=True)
noeuds_arc_sup = [idx for _, idx in noeuds_arc_sup]
noeuds_arc_inf = []
for point in noeuds_par_cercle[i + 1]:
theta, index = point["theta"], point["index"]
if index != noeud3 and index != noeud4:
if (est_entre_angles(theta, angle3, angle4, tolerance) or
angles_equivalents(theta, angle3, tolerance) or
angles_equivalents(theta, angle4, tolerance)):
noeuds_arc_inf.append((theta, index))
noeuds_arc_inf.sort()
noeuds_arc_inf = [idx for _, idx in noeuds_arc_inf]
# Construction des noeuds de la cellule dans l'ordre
tous_noeuds = [noeud1] + noeuds_arc_sup + [noeud2, noeud3] + noeuds_arc_inf + [noeud4]
cellule = {
"index": len(cellules_temp),
"noeuds": tous_noeuds,
"boundary": 0,
"cell_type": 0 if i == n - 2 else (3 if i == 0 else 0),
"target_area": 100.0,
"lambda_celllength": 0,
"at_boundary": True if i == n - 2 else False,
"dead": False,
"target_length": 0,
"stiffness": 1,
"source": False,
"pin_fixed": False,
"area": 0.0,
"fixed": False,
"div_counter": 0,
"cercle": i + 1,
"angle": angle_centre
}
cellules_temp.append(cellule)
else:
print(f"Impossible de créer la cellule à i={i}, j={j}: nœuds manquants")
print(f"Nœuds trouvés: {noeud1}, {noeud2}, {noeud3}, {noeud4}")
# Organiser les cellules par cercle puis par angle
cellules_par_cercle = {}
for cellule in cellules_temp:
cercle = cellule.get("cercle", 0)
if cercle not in cellules_par_cercle:
cellules_par_cercle[cercle] = []
cellules_par_cercle[cercle].append(cellule)
# Trier et réassembler les cellules
cellules_triees = []
for cercle in sorted(cellules_par_cercle.keys()):
cellules_par_cercle[cercle].sort(key=lambda c: c.get("angle", 0))
cellules_triees.extend(cellules_par_cercle[cercle])
# Réindexer les cellules triées
cellules = []
for i, cellule in enumerate(cellules_triees):
cellule_copy = cellule.copy()
cellule_copy["index"] = i
cellules.append(cellule_copy)
print("cellules", cellules)
return cellules
# ----------------------------------------------------------------------------------------------------#
# ---------- RAFFINNAGE DU MAILLAGE (AJOUT DE NOEUDS) --------------------------------------------------------------#
# ----------------------------------------------------------------------------------------------------#
def raffiner_maillage_separe(k_noeuds_n_cercle_polaire, n, cellules, longueur_seuil_cercle, longueur_seuil_radial,
rayons_n_cercles, seuil_proximite):
# Copier les données initiales
nouveau_k_noeuds_n_cercle_polaire = [[] for _ in range(n)]
for i in range(min(n, len(k_noeuds_n_cercle_polaire))):
for noeud in k_noeuds_n_cercle_polaire[i]:
nouveau_k_noeuds_n_cercle_polaire[i].append(noeud.copy())
nouvelles_cellules = [cellule.copy() for cellule in cellules]
# Tableau pour stocker les nœuds intermédiaires radiaux (qui ne sont pas sur les cercles)
noeuds_radiaux = []
# Identifier les noeuds avec tag spécial
noeuds_speciaux = {}
for i_cercle, cercle in enumerate(nouveau_k_noeuds_n_cercle_polaire):
for noeud in cercle:
if noeud.get("tag_special", False):
r = noeud["r"]
theta = noeud["theta"]
x = r * math.cos(theta)
y = r * math.sin(theta)
noeuds_speciaux[noeud["index"]] = (x, y, r, theta)
# Phase 1: Remaillage des cercles (segments circulaires)
iterations_cercles = 0
segments_cercles_raffines = True
while segments_cercles_raffines and iterations_cercles < 5:
iterations_cercles += 1
segments_cercles_raffines = False
# Convertir en coordonnées cartésiennes pour le calcul des distances
noeuds_cartesien = {}
for i_cercle in range(len(nouveau_k_noeuds_n_cercle_polaire)):
for noeud in nouveau_k_noeuds_n_cercle_polaire[i_cercle]:
r = noeud["r"]
theta = noeud["theta"]
index = noeud["index"]
initial = noeud["initial"]
boundary = noeud["boundary"]
x = r * math.cos(theta)
y = r * math.sin(theta)
noeuds_cartesien[index] = (x, y, r, theta, i_cercle, initial, boundary)
# Ajouter les noeuds radiaux à la liste des noeuds cartésiens
for noeud in noeuds_radiaux:
r = noeud["r"]
theta = noeud["theta"]
index = noeud["index"]
initial = noeud["initial"]
boundary = noeud["boundary"]
x = r * math.cos(theta)
y = r * math.sin(theta)
noeuds_cartesien[index] = (x, y, r, theta, -1, initial, boundary) # -1 indique un noeud radial
# Dictionnaire pour stocker les nouveaux noeuds à ajouter sur les cercles
segments_cercles_a_raffiner = {}
# Vérifier les segments sur un même cercle
for i_cercle in range(n):
if i_cercle >= len(nouveau_k_noeuds_n_cercle_polaire):
continue
cercle_noeuds = nouveau_k_noeuds_n_cercle_polaire[i_cercle]
nb_noeuds = len(cercle_noeuds)
for i in range(nb_noeuds):
idx1 = cercle_noeuds[i]["index"]
idx2 = cercle_noeuds[(i + 1) % nb_noeuds]["index"]
if idx1 in noeuds_cartesien and idx2 in noeuds_cartesien:
x1, y1, r1, theta1, _, _, boundary1 = noeuds_cartesien[idx1]
x2, y2, r2, theta2, _, _, boundary2 = noeuds_cartesien[idx2]
# Distance euclidienne entre les deux noeuds
distance = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
if distance > longueur_seuil_cercle:
# Gestion de l'interpolation angulaire
if theta1 < 0: theta1 += 2 * math.pi
if theta2 < 0: theta2 += 2 * math.pi
# Gérer le cas où l'arc traverse la ligne θ=0
if abs(theta1 - theta2) > math.pi:
if theta1 > theta2:
theta2 += 2 * math.pi
else:
theta1 += 2 * math.pi
# Interpolation angulaire en respectant le rayon du cercle
theta_new = (theta1 + theta2) / 2
r_new = rayons_n_cercles[i_cercle] # Utiliser le rayon exact du cercle
# Normaliser l'angle entre 0 et 2π
if theta_new >= 2 * math.pi:
theta_new -= 2 * math.pi
# Calculer les coordonnées cartésiennes
x_new = r_new * math.cos(theta_new)
y_new = r_new * math.sin(theta_new)
# Vérifier proximité avec noeuds spéciaux
trop_proche = False
for _, (x_special, y_special, _, _) in noeuds_speciaux.items():
dist_special = math.sqrt((x_new - x_special) ** 2 + (y_new - y_special) ** 2)
if dist_special < seuil_proximite:
trop_proche = True
break
if not trop_proche:
# Hériter le statut boundary si au moins un des nœuds est boundary
is_boundary = boundary1 or boundary2
segments_cercles_a_raffiner[tuple(sorted([idx1, idx2]))] = (x_new, y_new, r_new, theta_new, i_cercle, is_boundary)
segments_cercles_raffines = True
# Assigner des indices aux nouveaux noeuds
indices_max = [noeud["index"] for cercle in nouveau_k_noeuds_n_cercle_polaire for noeud in cercle]
indices_max.extend([noeud["index"] for noeud in noeuds_radiaux])
prochain_indice = max(indices_max) + 1 if indices_max else 0
# Ajouter les nouveaux noeuds aux cercles et créer le mappage
map_segments_cercles_indices = {}
for segment_key, (x_new, y_new, r_new, theta_new, cercle, is_boundary) in segments_cercles_a_raffiner.items():
idx1, idx2 = segment_key
nouveau_noeud = {
"r": r_new,
"theta": theta_new,
"index": prochain_indice,
"initial": False, # Nouveau nœud intermédiaire
"sam": False,
"boundary": is_boundary, # Hériter le statut boundary
"fixed": False
}
# Ajouter au cercle correspondant
nouveau_k_noeuds_n_cercle_polaire[cercle].append(nouveau_noeud)
# Mapper le segment au nouvel indice pour mise à jour des cellules
map_segments_cercles_indices[segment_key] = prochain_indice
prochain_indice += 1
# Trier les noeuds de chaque cercle par angle
for i in range(n):
if i < len(nouveau_k_noeuds_n_cercle_polaire):
nouveau_k_noeuds_n_cercle_polaire[i].sort(key=lambda point: point["theta"])
# Mettre à jour les cellules
maj_cellules = []
for cellule in nouvelles_cellules:
nouvelle_liste_noeuds = []
noeuds = cellule["noeuds"]
for i in range(len(noeuds)):
idx1 = noeuds[i]
idx2 = noeuds[(i + 1) % len(noeuds)]
nouvelle_liste_noeuds.append(idx1)
# Vérifier si ce segment a été raffiné (cercle)
segment_key = tuple(sorted([idx1, idx2]))
if segment_key in map_segments_cercles_indices:
nouvelle_liste_noeuds.append(map_segments_cercles_indices[segment_key])
nouvelle_cellule = cellule.copy()
nouvelle_cellule["noeuds"] = nouvelle_liste_noeuds
maj_cellules.append(nouvelle_cellule)
nouvelles_cellules = maj_cellules
# Phase 2: Remaillage des segments radiaux entre nœuds UNIQUEMENT où des segments radiaux existent déjà
iterations_radiaux = 0
segments_radiaux_raffines = True
while segments_radiaux_raffines and iterations_radiaux < 5:
iterations_radiaux += 1
segments_radiaux_raffines = False
# Convertir en coordonnées cartésiennes pour le calcul des distances
noeuds_cartesien = {}
for i_cercle in range(len(nouveau_k_noeuds_n_cercle_polaire)):
for noeud in nouveau_k_noeuds_n_cercle_polaire[i_cercle]:
r = noeud["r"]
theta = noeud["theta"]
index = noeud["index"]
initial = noeud["initial"]
boundary = noeud["boundary"]
x = r * math.cos(theta)
y = r * math.sin(theta)
noeuds_cartesien[index] = (x, y, r, theta, i_cercle, initial, boundary)
# Ajouter les noeuds radiaux à la liste des noeuds cartésiens
for noeud in noeuds_radiaux:
r = noeud["r"]
theta = noeud["theta"]
index = noeud["index"]
initial = noeud["initial"]
boundary = noeud["boundary"]
x = r * math.cos(theta)
y = r * math.sin(theta)
noeuds_cartesien[index] = (x, y, r, theta, -1, initial, boundary)
# Collecter les segments radiaux existants dans les cellules
segments_radiaux_existants = set()
for cellule in nouvelles_cellules:
noeuds = cellule["noeuds"]
for i in range(len(noeuds)):
idx1 = noeuds[i]
idx2 = noeuds[(i + 1) % len(noeuds)]
if idx1 in noeuds_cartesien and idx2 in noeuds_cartesien:
x1, y1, r1, theta1, cercle1, _, _ = noeuds_cartesien[idx1]
x2, y2, r2, theta2, cercle2, _, _ = noeuds_cartesien[idx2]
# Vérifier si c'est un segment radial (rayons différents)
if abs(r1 - r2) > 0.001: # Tolérance pour éviter erreurs d'arrondi
segments_radiaux_existants.add(tuple(sorted([idx1, idx2])))
# Dictionnaire pour stocker les nouveaux nœuds à ajouter
segments_radiaux_a_raffiner = {}
# Pour chaque segment radial existant
for segment in segments_radiaux_existants:
idx1, idx2 = segment
if idx1 in noeuds_cartesien and idx2 in noeuds_cartesien:
x1, y1, r1, theta1, cercle1, initial1, boundary1 = noeuds_cartesien[idx1]
x2, y2, r2, theta2, cercle2, initial2, boundary2 = noeuds_cartesien[idx2]
# Distance euclidienne entre les deux noeuds
distance = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
# Si la distance est supérieure au seuil, on ajoute un nœud intermédiaire
if distance > longueur_seuil_radial:
# Calculer la position du nouveau nœud intermédiaire
x_new = (x1 + x2) / 2
y_new = (y1 + y2) / 2
r_new = math.sqrt(x_new**2 + y_new**2)
theta_new = math.atan2(y_new, x_new)
if theta_new < 0:
theta_new += 2 * math.pi
# Hériter le statut boundary si l'un des nœuds est boundary
is_boundary = boundary1 or boundary2
# Sauvegarder les coordonnées du nouveau nœud
segments_radiaux_a_raffiner[segment] = (x_new, y_new, r_new, theta_new, is_boundary)
segments_radiaux_raffines = True
# Assigner des indices aux nouveaux noeuds radiaux
indices_max = [noeud["index"] for cercle in nouveau_k_noeuds_n_cercle_polaire for noeud in cercle]
indices_max.extend([noeud["index"] for noeud in noeuds_radiaux])
prochain_indice = max(indices_max) + 1 if indices_max else 0
# Créer les nouveaux noeuds radiaux
map_segments_radiaux_indices = {}
for segment_key, (x_new, y_new, r_new, theta_new, is_boundary) in segments_radiaux_a_raffiner.items():
idx1, idx2 = segment_key
nouveau_noeud = {
"r": r_new,
"theta": theta_new,
"index": prochain_indice,
"initial": False, # Nouveau nœud intermédiaire
"sam": False,
"boundary": False, # Hériter le statut boundary
"fixed": False
}
# Ajouter aux nœuds radiaux
noeuds_radiaux.append(nouveau_noeud)
# Mapper le segment au nouvel indice pour mise à jour des cellules
map_segments_radiaux_indices[segment_key] = prochain_indice
prochain_indice += 1
# Mettre à jour les cellules avec les nouveaux noeuds radiaux
maj_cellules = []
for cellule in nouvelles_cellules:
nouvelle_liste_noeuds = []
noeuds = cellule["noeuds"]
for i in range(len(noeuds)):
idx1 = noeuds[i]
idx2 = noeuds[(i + 1) % len(noeuds)]
nouvelle_liste_noeuds.append(idx1)
# Vérifier si ce segment a été raffiné (radialement)
segment_key = tuple(sorted([idx1, idx2]))
if segment_key in map_segments_radiaux_indices:
nouvelle_liste_noeuds.append(map_segments_radiaux_indices[segment_key])
nouvelle_cellule = cellule.copy()
nouvelle_cellule["noeuds"] = nouvelle_liste_noeuds
maj_cellules.append(nouvelle_cellule)
nouvelles_cellules = maj_cellules
k_noeuds_n_cercle_polaire_final = []
for cercle in nouveau_k_noeuds_n_cercle_polaire:
k_noeuds_n_cercle_polaire_final.append(cercle)
# Créer un cercle virtuel pour les noeuds radiaux (permet de rester compatible avec le format attendu)
if noeuds_radiaux:
k_noeuds_n_cercle_polaire_final.append(noeuds_radiaux)
return k_noeuds_n_cercle_polaire_final, nouvelles_cellules
#----------------------------------------------------------------------------------------------------#
#-------------------------------------walls-------------------------------------------#
#----------------------------------------------------------------------------------------------------#
def generer_walls_initiaux(k_noeuds_n_cercle_cartesien, cellules_list):
"""
Génère les walls (parois) entre les cellules adjacentes.
Args:
k_noeuds_n_cercle_cartesien: Liste des noeuds en coordonnées cartésiennes
cellules_list: Liste des cellules
Returns:
Liste des walls générés
"""
# Créer un dictionnaire pour accès rapide aux noeuds
noeuds_dict = {}
for noeud in k_noeuds_n_cercle_cartesien:
noeuds_dict[noeud["index"]] = (noeud["x"], noeud["y"], noeud.get("initial", False),
noeud.get("boundary", False))
# Identifier les segments et les cellules qui les partagent
segments_par_cellule = {}
for cellule in cellules_list:
noeuds = cellule["noeuds"]
cell_index = cellule["index"]
# Pour chaque paire de noeuds consécutifs
for i in range(len(noeuds)):
n1 = noeuds[i]
n2 = noeuds[(i + 1) % len(noeuds)]
segment = tuple(sorted([n1, n2]))
# Ajouter la cellule à ce segment
if segment not in segments_par_cellule:
segments_par_cellule[segment] = []
segments_par_cellule[segment].append(cell_index)
# Générer les walls uniquement pour les segments partagés par deux cellules
walls_liste = []
wall_index = 0
for segment, cellules in segments_par_cellule.items():
# Un wall ne peut exister qu'entre deux cellules
if len(cellules) == 2:
n1, n2 = segment
c1, c2 = sorted(cellules)
# Vérifier que les noeuds existent et sont initiaux
if n1 in noeuds_dict and n2 in noeuds_dict:
_, _, initial1, boundary1 = noeuds_dict[n1]
_, _, initial2, boundary2 = noeuds_dict[n2]
# Vérifier que les deux noeuds sont initiaux
if initial1 and initial2:
# Ne pas créer de wall si les deux noeuds sont sur la frontière
if boundary1 and boundary2:
continue
# Calculer la longueur du wall
(x1, y1, _, _), (x2, y2, _, _) = noeuds_dict[n1], noeuds_dict[n2]
longueur = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
# Créer le wall avec exactement la structure demandée
wall = {
"length": f"{longueur:.4f}",
"c1": str(c1),
"c2": str(c2),
"n1": str(n1),
"n2": str(n2),
"wall_type": "normal",
"viz_flux": "0",
"index": str(wall_index)
}
walls_liste.append(wall)
wall_index += 1
return walls_liste
#----------------------------------------------------------------------------------------------------#
#------------------------- Aire cellules ---------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
def shoelace_formula(cellules_list, noeuds_list, absoluteValue=True):
"""
Calcule l'aire de chaque cellule en utilisant la formule de Shoelace (Gauss).
Args:
cellules_list: Liste des cellules à traiter
noeuds_list: Liste des noeuds avec leurs coordonnées
absoluteValue: Si True, retourne la valeur absolue de l'aire
Returns:
Liste des cellules avec l'aire calculée
"""
# Optimisation par précalcul du dictionnaire de coordonnées
noeuds_dict = {}
for noeud in noeuds_list:
noeuds_dict[noeud["index"]] = (noeud["x"], noeud["y"])
# Préallocation de la liste de résultat
n_cellules = len(cellules_list)
cellules_avec_aire = [None] * n_cellules
# Traitement optimisé par cellule
for i, cellule in enumerate(cellules_list):
indices_noeuds = cellule["noeuds"]
nb_noeuds = len(indices_noeuds)
# Traitement rapide des polygones dégénérés
if nb_noeuds < 3:
cellule_maj = cellule.copy()
cellule_maj["area"] = 0.0
cellules_avec_aire[i] = cellule_maj
continue
# Extraction optimisée des coordonnées avec gestion d'erreur intégrée
try:
# Extraction en une seule fois pour éviter les accès répétés
coords = [noeuds_dict[idx] for idx in indices_noeuds]
# Boucle unifiée pour le calcul de l'aire (sans allocation supplémentaire)
somme = 0.0
for j in range(nb_noeuds):
k = (j + 1) % nb_noeuds
# Formule de Shoelace optimisée
somme += coords[j][0] * coords[k][1] - coords[k][0] * coords[j][1]
# Calcul final avec une seule multiplication
aire = abs(somme * 0.5) if absoluteValue else somme * 0.5
except KeyError:
# En cas d'indice invalide, attribuer une aire nulle
aire = 0.0
# Minimiser les copies en utilisant copy() au lieu de deepcopy()
cellule_maj = cellule.copy()
cellule_maj["area"] = aire
cellules_avec_aire[i] = cellule_maj
return cellules_avec_aire
def calculer_aire_totale_boundary(noeuds_list, absoluteValue=True):
"""
Calcule l'aire totale du tissu en utilisant uniquement les noeuds marqués comme boundary.
Version optimisée.
Args:
noeuds_list: Liste des noeuds avec leurs coordonnées
absoluteValue: Si True, retourne la valeur absolue de l'aire
Returns:
Aire totale du tissu délimité par les noeuds boundary
"""
# Extraction directe des noeuds boundary (plus efficace qu'une compréhension de liste)
noeuds_boundary = []
for noeud in noeuds_list:
if noeud.get("boundary", False):
noeuds_boundary.append(noeud)
# Vérification rapide du nombre de noeuds
n_noeuds = len(noeuds_boundary)
if n_noeuds < 3:
return 0.0
# Calcul du véritable centroïde (dans la version originale, il était initialisé à 0,0)
centre_x = sum(noeud["x"] for noeud in noeuds_boundary) / n_noeuds
centre_y = sum(noeud["y"] for noeud in noeuds_boundary) / n_noeuds
# Pré-calcul des angles pour le tri (évite les calculs répétés)