1 // Classes et fonctions XMeshLab
3 #include "MeshCut_Utils.hxx"
4 #include "MeshCut_Maillage.hxx"
6 #include "MeshCut_Carre.hxx"
7 #include "MeshCut_Cube.hxx"
9 #include "MeshCut_Fonctions.hxx"
10 #include "MeshCut_Cas.hxx"
12 #include "MeshCut_Globals.hxx"
19 using namespace MESHCUT;
22 // ================================== DECLARATION DES VARIABLES GLOBALES ==================================================
24 std::map<std::string, int> MESHCUT::intersections;
26 int MESHCUT::indexNouvellesMailles, MESHCUT::indexNouveauxNoeuds, MESHCUT::offsetMailles;
27 std::string MESHCUT::str_id_GMplus, MESHCUT::str_id_GMmoins;
28 Maillage *MESHCUT::MAILLAGE1, *MESHCUT::MAILLAGE2;
30 std::vector<float> MESHCUT::newXX, MESHCUT::newYY, MESHCUT::newZZ;
31 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::newCNX;
32 std::map<TYPE_MAILLE, int> MESHCUT::cptNouvellesMailles;
33 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::GMplus, MESHCUT::GMmoins;
38 std::string MESHCUT::str_id_maillagenew;
40 float MESHCUT::normale[3], MESHCUT::pointPlan[3];
42 float MESHCUT::epsilon;
47 // ================================== PROGRAMME PRINCIPAL ==================================================
49 int main(int argc, char *argv[])
65 throw std::exception();
66 char *ficMEDin0 = argv[1];
67 ficMEDin = (string) ficMEDin0;
68 char *ficMEDout0 = argv[2];
69 ficMEDout = (string) ficMEDout0;
70 char *id_maillagenew = argv[3];
71 str_id_maillagenew = (string) id_maillagenew;
74 char *id_GMplus = argv[4];
75 str_id_GMplus = (string) id_GMplus;
76 char *id_GMmoins = argv[5];
77 str_id_GMmoins = (string) id_GMmoins;
79 // Vecteur normal au plan de coupe
80 char *charxn = argv[6];
81 xNormal = char2float(charxn);
82 char *charyn = argv[7];
83 yNormal = char2float(charyn);
84 char *charzn = argv[8];
85 zNormal = char2float(charzn);
87 // Point du plan de coupe
88 char *charxm = argv[9];
89 xm = char2float(charxm);
90 char *charym = argv[10];
91 ym = char2float(charym);
92 char *charzm = argv[11];
93 zm = char2float(charzm);
95 // Tolérance : epsilon = tolérance * longueur arête moyenne - où epsilon est la tolérance absolue (distance)
96 char *chtolerance = argv[12];
97 tolerance = char2float(chtolerance);
102 cout << " Cut a tetrahedron mesh by a plane" << endl;
103 cout << " ---------------------------------" << endl;
104 cout << "Syntax:" << endl << endl;
105 cout << argv[0] << " input.med output.med resuMeshName aboveGroup belowGroup nx ny nz px py pz T " << endl;
106 cout << endl << "where:" << endl;
107 cout << " input.med = name of the original mesh file in med format" << endl;
108 cout << " output.med = name of the result mesh file in med format" << endl;
109 cout << " resuMeshName = name of the result mesh" << endl;
110 cout << " aboveGroup = name of the group of volumes above the cut plane" << endl;
111 cout << " belowGroups = name of the group of volumes below the cut plane" << endl;
112 cout << " nx ny nz = vector normal to the cut plane" << endl;
113 cout << " px py pz = a point of the cut plane" << endl;
114 cout << " T = 0 < T < 1 : vertices of a tetrahedron are considered as belonging" << endl;
115 cout << " the cut plane if their distance to the plane is inferior to L*T" << endl;
116 cout << " where L is the mean edge size of the tetrahedron" << endl;
117 ERREUR("--> check arguments!");
120 cout << "Découpe plane :" << endl;
121 cout << "\t" << "Maillage source: " << ficMEDin << endl;
122 cout << "\t" << "Maillage cible: " << ficMEDout << endl;
123 cout << "\t" << "ID nouveau maillage: " << str_id_maillagenew << endl;
124 cout << "\t" << "Nom GMplus: " << str_id_GMplus << endl;
125 cout << "\t" << "Nom GMmoins: " << str_id_GMmoins << endl;
126 cout << "\t" << "Vecteur normal du plan de coupe: xn=" << xNormal << " yn=" << yNormal << " zn=" << zNormal << endl;
127 cout << "\t" << "Point du plan de coupe: xm=" << xm << " ym=" << ym << " zm=" << zm << endl;
128 cout << "\t" << "Tolérance: " << tolerance << endl;
131 if (tolerance <= 0.0)
132 ERREUR("L'argument tolérance doit être strictement positif");
134 // Il faut normer la normale
135 float normeNormal = sqrt(xNormal * xNormal + yNormal * yNormal + zNormal * zNormal);
136 if (normeNormal == 0.0)
137 ERREUR("Vecteur normal nul");
138 normale[0] = xNormal / normeNormal;
139 normale[1] = yNormal / normeNormal;
140 normale[2] = zNormal / normeNormal;
146 // Calcul du coefficient d de l'équation du plan xn x + yn y + zn n + d = 0
147 d = -normale[0] * xm - normale[1] * ym - normale[2] * zm;
149 intersections.clear();
151 // Initialisation des compteurs de nouvelles mailles
152 for (int itm = (int) POI1; itm <= (int) HEXA20; itm++)
154 TYPE_MAILLE tm = (TYPE_MAILLE) itm;
155 cptNouvellesMailles[tm] = 0;
159 int S[4]; // Signature du T4 courant
160 int NG[4]; // Num. globaux des sommets
162 // Acquisition maillage initial
163 //cout << chrono() << " - Acquisition du maillage initial" << endl;
164 MAILLAGE1 = new Maillage((string) "TEMP");
165 MAILLAGE1->inputMED(ficMEDin);
166 cout << chrono() << " - Fin d'acquisition maillage" << endl;
167 indexNouveauxNoeuds = MAILLAGE1->nombreNoeudsMaillage;
169 // Le maillage ne contient aucun TETRA4 : on rend le maillage initial sans modification
170 if (!MAILLAGE1->EFFECTIFS_TYPES[TETRA4])
172 cout << "WARNING : le maillage ne contient aucun élément TETRA4, il ne sera donc pas modifié" << endl;
173 MAILLAGE1->ID = str_id_maillagenew;
174 MAILLAGE1->outputMED(ficMEDout);
175 cout << chrono() << " - Terminé" << endl << endl;
178 // A partir de cet instant le maillage contient forcément des TETRA4
181 // Chargement des distances noeud-plan DNP
182 DNP = (float*) malloc(sizeof(float) * MAILLAGE1->nombreNoeudsMaillage);
183 for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
184 DNP[k] = distanceNoeudPlan(k + 1);
185 cout << chrono() << " - Fin de chargement des distances noeud-plan" << endl;
187 // Longueur d'arête moyenne des T4 intersectant le plan de coupe
188 float LONGUEURS = 0.0;
189 int cptLONGUEURS = 0;
190 for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
194 int *offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
195 for (int is = 0; is < 4; is++)
197 int ng = *(offset + is);
198 if (DNP[ng - 1] > 0.0)
200 else if (DNP[ng - 1] < 0.0)
205 // Ce tetra est à cheval sur le plan de coupe: on calcule ses longueurs d'arêtes
206 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 1));
208 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 2));
210 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 3));
212 LONGUEURS += longueurSegment(*(offset + 1), *(offset + 2));
214 LONGUEURS += longueurSegment(*(offset + 1), *(offset + 3));
216 LONGUEURS += longueurSegment(*(offset + 2), *(offset + 3));
221 // Aucun TETRA4 intercepté par le plan de coupe : on rend MAILLAGE1
222 if (cptLONGUEURS == 0)
225 << "WARNING : le plan de coupe n'intersecte aucun élément TETRA4, le maillage initial ne sera donc pas modifié"
227 MAILLAGE1->ID = str_id_maillagenew;
228 MAILLAGE1->outputMED(ficMEDout);
229 cout << chrono() << " - Terminé" << endl << endl;
232 // A partir de cet instant le maillage contient forcément des TETRA4 intersectant le plan de coupe
235 float longueurMoyenne = LONGUEURS / cptLONGUEURS;
236 epsilon = tolerance * longueurMoyenne;
238 int nT4coupe = cptLONGUEURS / 6;
239 cout << chrono() << " - Fin du calcul de la longueur moyenne des arêtes TETRA4 au voisinage du plan de coupe" << endl;
241 cout << "Nombre de TETRA4 impliqués dans la coupe = " << nT4coupe << endl;
242 cout << "Longueur moyenne = " << longueurMoyenne << endl;
243 cout << "Tolérance = " << tolerance << endl;
244 cout << "Epsilon = " << epsilon << endl;
246 // Détermination des positions de noeuds par rapport au plan de coupe - POSN
247 POSN = (int*) malloc(sizeof(int) * MAILLAGE1->nombreNoeudsMaillage);
248 for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
250 if (DNP[k] > epsilon)
252 else if (DNP[k] < -epsilon)
257 cout << chrono() << " - Fin du positionnement des points par rapport au plan de coupe" << endl;
258 cout << "Début de la boucle sur les TETRA4" << endl;
260 for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
263 for (int is = 0; is < 4; is++)
265 int ng = *(MAILLAGE1->CNX[TETRA4] + 4 * it4 + is);
267 S[is] = *(POSN + ng - 1);
270 // -------------------------------------------------------------------
272 if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
273 GMmoins[TETRA4].push_back(it4);
275 else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
276 GMmoins[TETRA4].push_back(it4);
278 else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
279 { // Cas 3 - Arêtes 2 4 5
282 V[2] = intersectionSegmentPlan(it4, 2);
284 V[4] = intersectionSegmentPlan(it4, 4);
285 V[5] = intersectionSegmentPlan(it4, 5);
289 // -------------------------------------------------------------------
291 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
292 GMmoins[TETRA4].push_back(it4);
294 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
295 GMmoins[TETRA4].push_back(it4);
297 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
298 { // Cas 2, arêtes 2 4
301 V[2] = intersectionSegmentPlan(it4, 2);
303 V[4] = intersectionSegmentPlan(it4, 4);
308 // -------------------------------------------------------------------
310 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
311 { // Cas 3, arêtes 1 3 5
313 V[1] = intersectionSegmentPlan(it4, 1);
315 V[3] = intersectionSegmentPlan(it4, 3);
317 V[5] = intersectionSegmentPlan(it4, 5);
321 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
322 { // Cas 2, arêtes 1 3
324 V[1] = intersectionSegmentPlan(it4, 1);
326 V[3] = intersectionSegmentPlan(it4, 3);
332 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
333 { // Cas 4, arêtes 1 2 3 4
335 V[1] = intersectionSegmentPlan(it4, 1);
336 V[2] = intersectionSegmentPlan(it4, 2);
337 V[3] = intersectionSegmentPlan(it4, 3);
338 V[4] = intersectionSegmentPlan(it4, 4);
343 // -------------------------------------------------------------------
345 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
346 GMmoins[TETRA4].push_back(it4);
348 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
349 GMmoins[TETRA4].push_back(it4);
351 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
352 { // Cas 2, arêtes 2 5
355 V[2] = intersectionSegmentPlan(it4, 2);
358 V[5] = intersectionSegmentPlan(it4, 5);
362 // -------------------------------------------------------------------
364 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
365 GMmoins[TETRA4].push_back(it4);
367 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
368 GMmoins[TETRA4].push_back(it4);
370 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
374 V[2] = intersectionSegmentPlan(it4, 2);
381 // -------------------------------------------------------------------
383 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
384 { // Cas 2, arêtes 1 5
386 V[1] = intersectionSegmentPlan(it4, 1);
390 V[5] = intersectionSegmentPlan(it4, 5);
394 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
397 V[1] = intersectionSegmentPlan(it4, 1);
405 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
406 { // Cas 2, arêtes 1 2
408 V[1] = intersectionSegmentPlan(it4, 1);
409 V[2] = intersectionSegmentPlan(it4, 2);
416 // -------------------------------------------------------------------
418 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
419 { // Cas 3, arêtes 0 3 4
420 V[0] = intersectionSegmentPlan(it4, 0);
423 V[3] = intersectionSegmentPlan(it4, 3);
424 V[4] = intersectionSegmentPlan(it4, 4);
429 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
430 { // Cas 2, arêtes 0 3
431 V[0] = intersectionSegmentPlan(it4, 0);
434 V[3] = intersectionSegmentPlan(it4, 3);
440 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
441 { // Cas 4, arêtes 0 2 3 5
442 V[0] = intersectionSegmentPlan(it4, 0);
444 V[2] = intersectionSegmentPlan(it4, 2);
445 V[3] = intersectionSegmentPlan(it4, 3);
447 V[5] = intersectionSegmentPlan(it4, 5);
451 // -------------------------------------------------------------------
453 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
454 { // Cas 2, arêtes 0 4
455 V[0] = intersectionSegmentPlan(it4, 0);
459 V[4] = intersectionSegmentPlan(it4, 4);
464 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
466 V[0] = intersectionSegmentPlan(it4, 0);
475 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
476 { // Cas 2, arêtes 0 2
477 V[0] = intersectionSegmentPlan(it4, 0);
479 V[2] = intersectionSegmentPlan(it4, 2);
486 // -------------------------------------------------------------------
488 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
489 { // Cas 4, arêtes 0 1 4 5
490 V[0] = intersectionSegmentPlan(it4, 0);
491 V[1] = intersectionSegmentPlan(it4, 1);
494 V[4] = intersectionSegmentPlan(it4, 4);
495 V[5] = intersectionSegmentPlan(it4, 5);
499 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
500 { // Cas 2, arêtes 0 1
501 V[0] = intersectionSegmentPlan(it4, 0);
502 V[1] = intersectionSegmentPlan(it4, 1);
510 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
511 { // Cas 3, arêtes 0 1 2
512 V[0] = intersectionSegmentPlan(it4, 0);
513 V[1] = intersectionSegmentPlan(it4, 1);
514 V[2] = intersectionSegmentPlan(it4, 2);
521 // -------------------------------------------------------------------
523 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == -1)
524 GMmoins[TETRA4].push_back(it4);
526 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 0)
527 GMmoins[TETRA4].push_back(it4);
529 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 1)
530 { // Cas 2, arêtes 4 5
535 V[4] = intersectionSegmentPlan(it4, 4);
536 V[5] = intersectionSegmentPlan(it4, 5);
540 // -------------------------------------------------------------------
542 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == -1)
543 GMmoins[TETRA4].push_back(it4);
545 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 0)
546 GMmoins[TETRA4].push_back(it4);
548 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 1)
554 V[4] = intersectionSegmentPlan(it4, 4);
559 // -------------------------------------------------------------------
561 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == -1)
562 { // Cas 2, arêtes 3 5
566 V[3] = intersectionSegmentPlan(it4, 3);
568 V[5] = intersectionSegmentPlan(it4, 5);
572 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 0)
577 V[3] = intersectionSegmentPlan(it4, 3);
583 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 1)
584 { // Cas 2, arêtes 3 4
588 V[3] = intersectionSegmentPlan(it4, 3);
589 V[4] = intersectionSegmentPlan(it4, 4);
594 // -------------------------------------------------------------------
596 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == -1)
597 GMmoins[TETRA4].push_back(it4);
599 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 0)
600 GMmoins[TETRA4].push_back(it4);
602 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 1)
609 V[5] = intersectionSegmentPlan(it4, 5);
613 // -------------------------------------------------------------------
615 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == -1)
616 GMmoins[TETRA4].push_back(it4);
618 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 0)
620 cout << "WARNING : TETRA4 numéro " << it4
621 << " entièrement inclus dans la zone de tolérance autour du plan de coupe" << endl;
622 cout << " --> affecté au groupe " << str_id_GMmoins << endl;
623 GMmoins[TETRA4].push_back(it4);
626 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 1)
627 GMplus[TETRA4].push_back(it4);
629 // -------------------------------------------------------------------
631 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == -1)
638 V[5] = intersectionSegmentPlan(it4, 5);
642 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 0)
643 GMplus[TETRA4].push_back(it4);
645 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 1)
646 GMplus[TETRA4].push_back(it4);
648 // -------------------------------------------------------------------
650 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == -1)
651 { // Cas 2, arêtes 3 4
655 V[3] = intersectionSegmentPlan(it4, 3);
656 V[4] = intersectionSegmentPlan(it4, 4);
661 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 0)
666 V[3] = intersectionSegmentPlan(it4, 3);
672 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 1)
673 { // Cas 2, arêtes 3 5
677 V[3] = intersectionSegmentPlan(it4, 3);
679 V[5] = intersectionSegmentPlan(it4, 5);
683 // -------------------------------------------------------------------
685 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == -1)
691 V[4] = intersectionSegmentPlan(it4, 4);
696 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 0)
697 GMplus[TETRA4].push_back(it4);
699 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 1)
700 GMplus[TETRA4].push_back(it4);
702 // -------------------------------------------------------------------
704 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == -1)
705 { // Cas 2, arêtes 4 5
710 V[4] = intersectionSegmentPlan(it4, 4);
711 V[5] = intersectionSegmentPlan(it4, 5);
715 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 0)
716 GMplus[TETRA4].push_back(it4);
718 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 1)
719 GMplus[TETRA4].push_back(it4);
721 // -------------------------------------------------------------------
723 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
724 { // Cas 3, arêtes 0 1 2
725 V[0] = intersectionSegmentPlan(it4, 0);
726 V[1] = intersectionSegmentPlan(it4, 1);
727 V[2] = intersectionSegmentPlan(it4, 2);
734 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
735 { // Cas 2, arêtes 0 1
736 V[0] = intersectionSegmentPlan(it4, 0);
737 V[1] = intersectionSegmentPlan(it4, 1);
745 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
746 { // Cas 4, arêtes 0 1 4 5
747 V[0] = intersectionSegmentPlan(it4, 0);
748 V[1] = intersectionSegmentPlan(it4, 1);
751 V[4] = intersectionSegmentPlan(it4, 4);
752 V[5] = intersectionSegmentPlan(it4, 5);
756 // -------------------------------------------------------------------
758 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
759 { // Cas 2, arêtes 0 2
760 V[0] = intersectionSegmentPlan(it4, 0);
762 V[2] = intersectionSegmentPlan(it4, 2);
769 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
771 V[0] = intersectionSegmentPlan(it4, 0);
780 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
781 { // Cas 2, arêtes 0 4
782 V[0] = intersectionSegmentPlan(it4, 0);
786 V[4] = intersectionSegmentPlan(it4, 4);
791 // -------------------------------------------------------------------
793 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
794 { // Cas 4, arêtes 0 2 3 5
795 V[0] = intersectionSegmentPlan(it4, 0);
797 V[2] = intersectionSegmentPlan(it4, 2);
798 V[3] = intersectionSegmentPlan(it4, 3);
800 V[5] = intersectionSegmentPlan(it4, 5);
804 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
805 { // Cas 2, arêtes 0 3
806 V[0] = intersectionSegmentPlan(it4, 0);
809 V[3] = intersectionSegmentPlan(it4, 3);
815 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
816 { // Cas 3, arêtes 0 3 4
817 V[0] = intersectionSegmentPlan(it4, 0);
820 V[3] = intersectionSegmentPlan(it4, 3);
821 V[4] = intersectionSegmentPlan(it4, 4);
826 // -------------------------------------------------------------------
828 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
829 { // Cas 2, arêtes 1 2
831 V[1] = intersectionSegmentPlan(it4, 1);
832 V[2] = intersectionSegmentPlan(it4, 2);
839 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
842 V[1] = intersectionSegmentPlan(it4, 1);
850 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
851 { // Cas 2, arêtes 1 5
853 V[1] = intersectionSegmentPlan(it4, 1);
857 V[5] = intersectionSegmentPlan(it4, 5);
861 // -------------------------------------------------------------------
863 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
867 V[2] = intersectionSegmentPlan(it4, 2);
874 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
875 GMplus[TETRA4].push_back(it4);
877 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
878 GMplus[TETRA4].push_back(it4);
880 // -------------------------------------------------------------------
882 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
883 { // Cas 2, arêtes 2 5
886 V[2] = intersectionSegmentPlan(it4, 2);
889 V[5] = intersectionSegmentPlan(it4, 5);
893 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
894 GMplus[TETRA4].push_back(it4);
896 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
897 GMplus[TETRA4].push_back(it4);
899 // -------------------------------------------------------------------
901 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
902 { // Cas 4, arêtes 1 2 3 4
904 V[1] = intersectionSegmentPlan(it4, 1);
905 V[2] = intersectionSegmentPlan(it4, 2);
906 V[3] = intersectionSegmentPlan(it4, 3);
907 V[4] = intersectionSegmentPlan(it4, 4);
912 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
913 { // Cas 2, arêtes 1 3
915 V[1] = intersectionSegmentPlan(it4, 1);
917 V[3] = intersectionSegmentPlan(it4, 3);
923 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
924 { // Cas 3, arêtes 1 3 5
926 V[1] = intersectionSegmentPlan(it4, 1);
928 V[3] = intersectionSegmentPlan(it4, 3);
930 V[5] = intersectionSegmentPlan(it4, 5);
934 // -------------------------------------------------------------------
936 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
937 { // Cas 2, arêtes 2 4
940 V[2] = intersectionSegmentPlan(it4, 2);
942 V[4] = intersectionSegmentPlan(it4, 4);
947 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
948 GMplus[TETRA4].push_back(it4);
950 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
951 GMplus[TETRA4].push_back(it4);
953 // -------------------------------------------------------------------
955 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
956 { // Cas 3, arêtes 2 4 5
959 V[2] = intersectionSegmentPlan(it4, 2);
961 V[4] = intersectionSegmentPlan(it4, 4);
962 V[5] = intersectionSegmentPlan(it4, 5);
966 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
967 GMplus[TETRA4].push_back(it4);
969 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
970 GMplus[TETRA4].push_back(it4);
973 ERREUR("Cas non prévu");
976 cout << chrono() << " - Fin de la boucle sur les TETRA4" << endl;
978 // cout << "indexNouveauxNoeuds = " << indexNouveauxNoeuds << endl;
979 newXX.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
980 newYY.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
981 newZZ.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
983 if (cptNouvellesMailles[TETRA4])
984 newCNX[TETRA4].resize(4 * cptNouvellesMailles[TETRA4]);
985 if (cptNouvellesMailles[PYRAM5])
986 newCNX[PYRAM5].resize(5 * cptNouvellesMailles[PYRAM5]);
987 if (cptNouvellesMailles[PENTA6])
988 newCNX[PENTA6].resize(6 * cptNouvellesMailles[PENTA6]);
990 // =========================================================================================
991 // 2. Constitution du maillage final
992 // =========================================================================================
994 cout << chrono() << " - Constitution du maillage final" << endl;
996 MAILLAGE2 = new Maillage(str_id_maillagenew);
997 MAILLAGE2->dimensionMaillage = MAILLAGE1->dimensionMaillage;
998 MAILLAGE2->dimensionEspace = MAILLAGE1->dimensionEspace;
999 strcpy(MAILLAGE2->axisname, MAILLAGE1->axisname);
1000 strcpy(MAILLAGE2->unitname, MAILLAGE1->unitname);
1001 MAILLAGE2->nombreNoeudsMaillage = indexNouveauxNoeuds;
1002 MAILLAGE2->nombreMaillesMaillage = MAILLAGE1->nombreMaillesMaillage + cptNouvellesMailles[TETRA4]
1003 + cptNouvellesMailles[PYRAM5] + cptNouvellesMailles[PENTA6];
1005 // ---------- Coordonnées
1006 // Optimisation de la mémoire au détriment du temps
1008 // Héritage des coordonnées MAILLAGE1
1009 MAILLAGE2->XX = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1010 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1011 *(MAILLAGE2->XX + i) = *(MAILLAGE1->XX + i);
1012 free(MAILLAGE1->XX);
1013 MAILLAGE2->YY = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1014 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1015 *(MAILLAGE2->YY + i) = *(MAILLAGE1->YY + i);
1016 free(MAILLAGE1->YY);
1017 MAILLAGE2->ZZ = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1018 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1019 *(MAILLAGE2->ZZ + i) = *(MAILLAGE1->ZZ + i);
1020 free(MAILLAGE1->ZZ);
1022 // Coordonnées des noeuds créés
1023 for (int i = 0; i < MAILLAGE2->nombreNoeudsMaillage - MAILLAGE1->nombreNoeudsMaillage; i++)
1025 *(MAILLAGE2->XX + MAILLAGE1->nombreNoeudsMaillage + i) = newXX[i];
1026 *(MAILLAGE2->YY + MAILLAGE1->nombreNoeudsMaillage + i) = newYY[i];
1027 *(MAILLAGE2->ZZ + MAILLAGE1->nombreNoeudsMaillage + i) = newZZ[i];
1028 // cout << "Nouveaux noeuds, indice " << i << " : " << newXX[i] << " " << newYY[i] << " " << newZZ[i] << " " << endl;
1031 // Legacy mailles maillage 1 +
1032 for (int itm = (int) POI1; itm <= (int) HEXA20; itm++)
1034 TYPE_MAILLE tm = (TYPE_MAILLE) itm;
1035 if (tm != TETRA4 && tm != PYRAM5 && tm != PENTA6)
1037 // Pour les types autres que TETRA4 PYRAM5 PENTA6 on fait seulement pointer CNX2 vers CNX1
1038 if (MAILLAGE1->EFFECTIFS_TYPES[tm])
1039 MAILLAGE2->CNX[tm] = MAILLAGE1->CNX[tm];
1040 MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm];
1044 // Pour les types TETRA4 PYRAM5 PENTA6 on recopie CNX1 et on ajoute à la suite les newCNX
1045 // cout << "Legacy " << tm << " effectif " << MAILLAGE1->EFFECTIFS_TYPES[tm] << endl;
1046 int tailleType = Nnoeuds(tm);
1048 MAILLAGE2->CNX[tm] = (int*) malloc(sizeof(int) * tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm]
1049 + cptNouvellesMailles[tm]));
1050 for (int i = 0; i < MAILLAGE1->EFFECTIFS_TYPES[tm]; i++)
1051 for (int j = 0; j < tailleType; j++)
1052 *(MAILLAGE2->CNX[tm] + tailleType * i + j) = *(MAILLAGE1->CNX[tm] + tailleType * i + j);
1054 for (int i = 0; i < cptNouvellesMailles[tm]; i++)
1055 for (int j = 0; j < tailleType; j++)
1056 *(MAILLAGE2->CNX[tm] + tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm] + i) + j) = newCNX[tm][i * tailleType
1059 MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm] + cptNouvellesMailles[tm];
1065 // cout << "Maillage 2 - CNX TETRA4 : " << endl;
1067 // for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[TETRA4]; i++)
1069 // cout << "Maille " << i << " : ";
1070 // for (int j = 0; j < 4; j++)
1071 // cout << MAILLAGE2->CNX[TETRA4][i * 4 + j] << " ";
1075 // cout << "Maillage 2 - CNX PENTA6 : " << endl;
1077 // for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[PENTA6]; i++)
1079 // cout << "Maille " << i << " : ";
1080 // for (int j = 0; j < 6; j++)
1081 // cout << MAILLAGE2->CNX[PENTA6][i * 6 + j] << " ";
1086 // Groupes de mailles
1087 MAILLAGE2->GM = MAILLAGE1->GM;
1088 MAILLAGE2->GM[str_id_GMplus] = GMplus;
1089 MAILLAGE2->GM[str_id_GMmoins] = GMmoins;
1091 MAILLAGE2->GN = MAILLAGE1->GN;
1093 cout << chrono() << " - Ecriture du fichier MED" << endl;
1095 MAILLAGE2->outputMED(ficMEDout);
1096 cout << chrono() << " - Terminé" << endl << endl;