1 // Copyright (C) 2006-2013 EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // Classes et fonctions XMeshLab
22 #include "MeshCut_Utils.hxx"
23 #include "MeshCut_Maillage.hxx"
25 #include "MeshCut_Carre.hxx"
26 #include "MeshCut_Cube.hxx"
28 #include "MeshCut_Fonctions.hxx"
29 #include "MeshCut_Cas.hxx"
31 #include "MeshCut_Globals.hxx"
38 using namespace MESHCUT;
41 // ================================== DECLARATION DES VARIABLES GLOBALES ==================================================
43 std::map<std::string, int> MESHCUT::intersections;
45 int MESHCUT::indexNouvellesMailles, MESHCUT::indexNouveauxNoeuds, MESHCUT::offsetMailles;
46 std::string MESHCUT::str_id_GMplus, MESHCUT::str_id_GMmoins;
47 Maillage *MESHCUT::MAILLAGE1, *MESHCUT::MAILLAGE2;
49 std::vector<float> MESHCUT::newXX, MESHCUT::newYY, MESHCUT::newZZ;
50 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::newCNX;
51 std::map<TYPE_MAILLE, int> MESHCUT::cptNouvellesMailles;
52 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::GMplus, MESHCUT::GMmoins;
53 std::vector<int> MESHCUT::cutTetras;
58 std::string MESHCUT::str_id_maillagenew;
60 float MESHCUT::normale[3], MESHCUT::pointPlan[3];
62 float MESHCUT::epsilon;
67 // ================================== PROGRAMME PRINCIPAL ==================================================
69 int main(int argc, char *argv[])
85 throw std::exception();
86 char *ficMEDin0 = argv[1];
87 ficMEDin = (string) ficMEDin0;
88 char *ficMEDout0 = argv[2];
89 ficMEDout = (string) ficMEDout0;
90 char *id_maillagenew = argv[3];
91 str_id_maillagenew = (string) id_maillagenew;
94 char *id_GMplus = argv[4];
95 str_id_GMplus = (string) id_GMplus;
96 char *id_GMmoins = argv[5];
97 str_id_GMmoins = (string) id_GMmoins;
99 // Vecteur normal au plan de coupe
100 char *charxn = argv[6];
101 xNormal = char2float(charxn);
102 char *charyn = argv[7];
103 yNormal = char2float(charyn);
104 char *charzn = argv[8];
105 zNormal = char2float(charzn);
107 // Point du plan de coupe
108 char *charxm = argv[9];
109 xm = char2float(charxm);
110 char *charym = argv[10];
111 ym = char2float(charym);
112 char *charzm = argv[11];
113 zm = char2float(charzm);
115 // Tolérance : epsilon = tolérance * longueur arête moyenne - où epsilon est la tolérance absolue (distance)
116 char *chtolerance = argv[12];
117 tolerance = char2float(chtolerance);
122 cout << " Cut a tetrahedron mesh by a plane" << endl;
123 cout << " ---------------------------------" << endl;
124 cout << "Syntax:" << endl << endl;
125 cout << argv[0] << " input.med output.med resuMeshName aboveGroup belowGroup nx ny nz px py pz T " << endl;
126 cout << endl << "where:" << endl;
127 cout << " input.med = name of the original mesh file in med format" << endl;
128 cout << " output.med = name of the result mesh file in med format" << endl;
129 cout << " resuMeshName = name of the result mesh" << endl;
130 cout << " aboveGroup = name of the group of volumes above the cut plane" << endl;
131 cout << " belowGroups = name of the group of volumes below the cut plane" << endl;
132 cout << " nx ny nz = vector normal to the cut plane" << endl;
133 cout << " px py pz = a point of the cut plane" << endl;
134 cout << " T = 0 < T < 1 : vertices of a tetrahedron are considered as belonging" << endl;
135 cout << " the cut plane if their distance to the plane is inferior to L*T" << endl;
136 cout << " where L is the mean edge size of the tetrahedron" << endl;
137 ERREUR("--> check arguments!");
140 cout << "Cut by a plane :" << endl;
141 cout << " source mesh: " << ficMEDin << endl;
142 cout << " result mesh: " << ficMEDout << endl;
143 cout << " mesh name: " << str_id_maillagenew << endl;
144 cout << " group above plane: " << str_id_GMplus << endl;
145 cout << " group below plane: " << str_id_GMmoins << endl;
146 cout << " vector normal to the cut plane: xn=" << xNormal << " yn=" << yNormal << " zn=" << zNormal << endl;
147 cout << " point in the cut plane: xm=" << xm << " ym=" << ym << " zm=" << zm << endl;
148 cout << " tolerance: " << tolerance << endl;
151 if (tolerance <= 0.0)
152 ERREUR("Tolerance must not be negative or null");
154 // Il faut normer la normale
155 float normeNormal = sqrt(xNormal * xNormal + yNormal * yNormal + zNormal * zNormal);
156 if (normeNormal == 0.0)
157 ERREUR("null normal vector");
158 normale[0] = xNormal / normeNormal;
159 normale[1] = yNormal / normeNormal;
160 normale[2] = zNormal / normeNormal;
166 // Calcul du coefficient d de l'équation du plan xn x + yn y + zn n + d = 0
167 d = -normale[0] * xm - normale[1] * ym - normale[2] * zm;
169 intersections.clear();
171 // Initialisation des compteurs de nouvelles mailles
172 for (int itm = (int) POI1; itm <= (int) HEXA20; itm++)
174 TYPE_MAILLE tm = (TYPE_MAILLE) itm;
175 cptNouvellesMailles[tm] = 0;
179 int S[4]; // Signature du T4 courant
180 //int NG[4]; // Num. globaux des sommets
182 // Acquisition maillage initial
183 //cout << chrono() << " - Acquisition du maillage initial" << endl;
184 MAILLAGE1 = new Maillage((string) "TEMP");
185 MAILLAGE1->inputMED(ficMEDin);
186 cout << chrono() << " - End of mesh read" << endl;
187 indexNouveauxNoeuds = MAILLAGE1->nombreNoeudsMaillage;
189 // Le maillage ne contient aucun TETRA4 : on rend le maillage initial sans modification
190 if (!MAILLAGE1->EFFECTIFS_TYPES[TETRA4])
192 cout << "WARNING: mesh does not contain tetra4 elements, it will not be modified" << endl;
193 MAILLAGE1->ID = str_id_maillagenew;
194 MAILLAGE1->outputMED(ficMEDout);
195 cout << chrono() << " - Finished!" << endl << endl;
198 // A partir de cet instant le maillage contient forcément des TETRA4
201 // Chargement des distances noeud-plan DNP
202 DNP = (float*) malloc(sizeof(float) * MAILLAGE1->nombreNoeudsMaillage);
203 for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
204 DNP[k] = distanceNoeudPlan(k + 1);
205 cout << chrono() << " - End of computation of distances between nodes and plane" << endl;
207 // Longueur d'arête moyenne des T4 intersectant le plan de coupe
208 float LONGUEURS = 0.0;
209 int cptLONGUEURS = 0;
210 for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
214 int *offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
215 for (int is = 0; is < 4; is++)
217 int ng = *(offset + is);
218 if (DNP[ng - 1] > 0.0)
220 else if (DNP[ng - 1] < 0.0)
225 // Ce tetra est à cheval sur le plan de coupe: on calcule ses longueurs d'arêtes
226 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 1));
228 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 2));
230 LONGUEURS += longueurSegment(*(offset + 0), *(offset + 3));
232 LONGUEURS += longueurSegment(*(offset + 1), *(offset + 2));
234 LONGUEURS += longueurSegment(*(offset + 1), *(offset + 3));
236 LONGUEURS += longueurSegment(*(offset + 2), *(offset + 3));
241 // Aucun TETRA4 intercepté par le plan de coupe : on rend MAILLAGE1
242 if (cptLONGUEURS == 0)
245 << "WARNING: the cut plane does not cut any tetra4 element, initial mesh will not be modified"
247 MAILLAGE1->ID = str_id_maillagenew;
248 MAILLAGE1->outputMED(ficMEDout);
249 cout << chrono() << " - Finished!" << endl << endl;
252 // A partir de cet instant le maillage contient forcément des TETRA4 intersectant le plan de coupe
255 float longueurMoyenne = LONGUEURS / cptLONGUEURS;
256 epsilon = tolerance * longueurMoyenne;
258 int nT4coupe = cptLONGUEURS / 6;
259 cout << chrono() << " - End of computation of mean length of tetra4 edges near the cut plane" << endl;
261 cout << "Number of tetra4 to be cut = " << nT4coupe << endl;
262 cout << "Mean length = " << longueurMoyenne << endl;
263 cout << "Tolerance = " << tolerance << endl;
264 cout << "Epsilon = " << epsilon << endl;
266 // Détermination des positions de noeuds par rapport au plan de coupe - POSN
267 POSN = (int*) malloc(sizeof(int) * MAILLAGE1->nombreNoeudsMaillage);
268 for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
270 if (DNP[k] > epsilon)
272 else if (DNP[k] < -epsilon)
277 cout << chrono() << " - End of nodes qualification above or below the cut plane" << endl;
278 cout << "Start of iteration on tetra4" << endl;
280 for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
283 for (int is = 0; is < 4; is++)
285 int ng = *(MAILLAGE1->CNX[TETRA4] + 4 * it4 + is);
287 S[is] = *(POSN + ng - 1);
290 // -------------------------------------------------------------------
292 if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
293 GMmoins[TETRA4].push_back(it4);
295 else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
296 GMmoins[TETRA4].push_back(it4);
298 else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
299 { // Cas 3 - Arêtes 2 4 5
302 V[2] = intersectionSegmentPlan(it4, 2);
304 V[4] = intersectionSegmentPlan(it4, 4);
305 V[5] = intersectionSegmentPlan(it4, 5);
309 // -------------------------------------------------------------------
311 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
312 GMmoins[TETRA4].push_back(it4);
314 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
315 GMmoins[TETRA4].push_back(it4);
317 else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
318 { // Cas 2, arêtes 2 4
321 V[2] = intersectionSegmentPlan(it4, 2);
323 V[4] = intersectionSegmentPlan(it4, 4);
328 // -------------------------------------------------------------------
330 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
331 { // Cas 3, arêtes 1 3 5
333 V[1] = intersectionSegmentPlan(it4, 1);
335 V[3] = intersectionSegmentPlan(it4, 3);
337 V[5] = intersectionSegmentPlan(it4, 5);
341 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
342 { // Cas 2, arêtes 1 3
344 V[1] = intersectionSegmentPlan(it4, 1);
346 V[3] = intersectionSegmentPlan(it4, 3);
352 else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
353 { // Cas 4, arêtes 1 2 3 4
355 V[1] = intersectionSegmentPlan(it4, 1);
356 V[2] = intersectionSegmentPlan(it4, 2);
357 V[3] = intersectionSegmentPlan(it4, 3);
358 V[4] = intersectionSegmentPlan(it4, 4);
363 // -------------------------------------------------------------------
365 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
366 GMmoins[TETRA4].push_back(it4);
368 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
369 GMmoins[TETRA4].push_back(it4);
371 else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
372 { // Cas 2, arêtes 2 5
375 V[2] = intersectionSegmentPlan(it4, 2);
378 V[5] = intersectionSegmentPlan(it4, 5);
382 // -------------------------------------------------------------------
384 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
385 GMmoins[TETRA4].push_back(it4);
387 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
388 GMmoins[TETRA4].push_back(it4);
390 else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
394 V[2] = intersectionSegmentPlan(it4, 2);
401 // -------------------------------------------------------------------
403 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
404 { // Cas 2, arêtes 1 5
406 V[1] = intersectionSegmentPlan(it4, 1);
410 V[5] = intersectionSegmentPlan(it4, 5);
414 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
417 V[1] = intersectionSegmentPlan(it4, 1);
425 else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
426 { // Cas 2, arêtes 1 2
428 V[1] = intersectionSegmentPlan(it4, 1);
429 V[2] = intersectionSegmentPlan(it4, 2);
436 // -------------------------------------------------------------------
438 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
439 { // Cas 3, arêtes 0 3 4
440 V[0] = intersectionSegmentPlan(it4, 0);
443 V[3] = intersectionSegmentPlan(it4, 3);
444 V[4] = intersectionSegmentPlan(it4, 4);
449 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
450 { // Cas 2, arêtes 0 3
451 V[0] = intersectionSegmentPlan(it4, 0);
454 V[3] = intersectionSegmentPlan(it4, 3);
460 else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
461 { // Cas 4, arêtes 0 2 3 5
462 V[0] = intersectionSegmentPlan(it4, 0);
464 V[2] = intersectionSegmentPlan(it4, 2);
465 V[3] = intersectionSegmentPlan(it4, 3);
467 V[5] = intersectionSegmentPlan(it4, 5);
471 // -------------------------------------------------------------------
473 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
474 { // Cas 2, arêtes 0 4
475 V[0] = intersectionSegmentPlan(it4, 0);
479 V[4] = intersectionSegmentPlan(it4, 4);
484 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
486 V[0] = intersectionSegmentPlan(it4, 0);
495 else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
496 { // Cas 2, arêtes 0 2
497 V[0] = intersectionSegmentPlan(it4, 0);
499 V[2] = intersectionSegmentPlan(it4, 2);
506 // -------------------------------------------------------------------
508 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
509 { // Cas 4, arêtes 0 1 4 5
510 V[0] = intersectionSegmentPlan(it4, 0);
511 V[1] = intersectionSegmentPlan(it4, 1);
514 V[4] = intersectionSegmentPlan(it4, 4);
515 V[5] = intersectionSegmentPlan(it4, 5);
519 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
520 { // Cas 2, arêtes 0 1
521 V[0] = intersectionSegmentPlan(it4, 0);
522 V[1] = intersectionSegmentPlan(it4, 1);
530 else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
531 { // Cas 3, arêtes 0 1 2
532 V[0] = intersectionSegmentPlan(it4, 0);
533 V[1] = intersectionSegmentPlan(it4, 1);
534 V[2] = intersectionSegmentPlan(it4, 2);
541 // -------------------------------------------------------------------
543 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == -1)
544 GMmoins[TETRA4].push_back(it4);
546 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 0)
547 GMmoins[TETRA4].push_back(it4);
549 else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 1)
550 { // Cas 2, arêtes 4 5
555 V[4] = intersectionSegmentPlan(it4, 4);
556 V[5] = intersectionSegmentPlan(it4, 5);
560 // -------------------------------------------------------------------
562 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == -1)
563 GMmoins[TETRA4].push_back(it4);
565 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 0)
566 GMmoins[TETRA4].push_back(it4);
568 else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 1)
574 V[4] = intersectionSegmentPlan(it4, 4);
579 // -------------------------------------------------------------------
581 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == -1)
582 { // Cas 2, arêtes 3 5
586 V[3] = intersectionSegmentPlan(it4, 3);
588 V[5] = intersectionSegmentPlan(it4, 5);
592 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 0)
597 V[3] = intersectionSegmentPlan(it4, 3);
603 else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 1)
604 { // Cas 2, arêtes 3 4
608 V[3] = intersectionSegmentPlan(it4, 3);
609 V[4] = intersectionSegmentPlan(it4, 4);
614 // -------------------------------------------------------------------
616 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == -1)
617 GMmoins[TETRA4].push_back(it4);
619 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 0)
620 GMmoins[TETRA4].push_back(it4);
622 else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 1)
629 V[5] = intersectionSegmentPlan(it4, 5);
633 // -------------------------------------------------------------------
635 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == -1)
636 GMmoins[TETRA4].push_back(it4);
638 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 0)
640 cout << "WARNING: TETRA4 number " << it4
641 << " entirely in the tolerance zone near the cut plane" << endl;
642 cout << " --> affected to group " << str_id_GMmoins << endl;
643 GMmoins[TETRA4].push_back(it4);
646 else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 1)
647 GMplus[TETRA4].push_back(it4);
649 // -------------------------------------------------------------------
651 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == -1)
658 V[5] = intersectionSegmentPlan(it4, 5);
662 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 0)
663 GMplus[TETRA4].push_back(it4);
665 else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 1)
666 GMplus[TETRA4].push_back(it4);
668 // -------------------------------------------------------------------
670 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == -1)
671 { // Cas 2, arêtes 3 4
675 V[3] = intersectionSegmentPlan(it4, 3);
676 V[4] = intersectionSegmentPlan(it4, 4);
681 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 0)
686 V[3] = intersectionSegmentPlan(it4, 3);
692 else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 1)
693 { // Cas 2, arêtes 3 5
697 V[3] = intersectionSegmentPlan(it4, 3);
699 V[5] = intersectionSegmentPlan(it4, 5);
703 // -------------------------------------------------------------------
705 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == -1)
711 V[4] = intersectionSegmentPlan(it4, 4);
716 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 0)
717 GMplus[TETRA4].push_back(it4);
719 else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 1)
720 GMplus[TETRA4].push_back(it4);
722 // -------------------------------------------------------------------
724 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == -1)
725 { // Cas 2, arêtes 4 5
730 V[4] = intersectionSegmentPlan(it4, 4);
731 V[5] = intersectionSegmentPlan(it4, 5);
735 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 0)
736 GMplus[TETRA4].push_back(it4);
738 else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 1)
739 GMplus[TETRA4].push_back(it4);
741 // -------------------------------------------------------------------
743 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
744 { // Cas 3, arêtes 0 1 2
745 V[0] = intersectionSegmentPlan(it4, 0);
746 V[1] = intersectionSegmentPlan(it4, 1);
747 V[2] = intersectionSegmentPlan(it4, 2);
754 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
755 { // Cas 2, arêtes 0 1
756 V[0] = intersectionSegmentPlan(it4, 0);
757 V[1] = intersectionSegmentPlan(it4, 1);
765 else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
766 { // Cas 4, arêtes 0 1 4 5
767 V[0] = intersectionSegmentPlan(it4, 0);
768 V[1] = intersectionSegmentPlan(it4, 1);
771 V[4] = intersectionSegmentPlan(it4, 4);
772 V[5] = intersectionSegmentPlan(it4, 5);
776 // -------------------------------------------------------------------
778 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
779 { // Cas 2, arêtes 0 2
780 V[0] = intersectionSegmentPlan(it4, 0);
782 V[2] = intersectionSegmentPlan(it4, 2);
789 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
791 V[0] = intersectionSegmentPlan(it4, 0);
800 else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
801 { // Cas 2, arêtes 0 4
802 V[0] = intersectionSegmentPlan(it4, 0);
806 V[4] = intersectionSegmentPlan(it4, 4);
811 // -------------------------------------------------------------------
813 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
814 { // Cas 4, arêtes 0 2 3 5
815 V[0] = intersectionSegmentPlan(it4, 0);
817 V[2] = intersectionSegmentPlan(it4, 2);
818 V[3] = intersectionSegmentPlan(it4, 3);
820 V[5] = intersectionSegmentPlan(it4, 5);
824 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
825 { // Cas 2, arêtes 0 3
826 V[0] = intersectionSegmentPlan(it4, 0);
829 V[3] = intersectionSegmentPlan(it4, 3);
835 else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
836 { // Cas 3, arêtes 0 3 4
837 V[0] = intersectionSegmentPlan(it4, 0);
840 V[3] = intersectionSegmentPlan(it4, 3);
841 V[4] = intersectionSegmentPlan(it4, 4);
846 // -------------------------------------------------------------------
848 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
849 { // Cas 2, arêtes 1 2
851 V[1] = intersectionSegmentPlan(it4, 1);
852 V[2] = intersectionSegmentPlan(it4, 2);
859 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
862 V[1] = intersectionSegmentPlan(it4, 1);
870 else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
871 { // Cas 2, arêtes 1 5
873 V[1] = intersectionSegmentPlan(it4, 1);
877 V[5] = intersectionSegmentPlan(it4, 5);
881 // -------------------------------------------------------------------
883 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
887 V[2] = intersectionSegmentPlan(it4, 2);
894 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
895 GMplus[TETRA4].push_back(it4);
897 else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
898 GMplus[TETRA4].push_back(it4);
900 // -------------------------------------------------------------------
902 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
903 { // Cas 2, arêtes 2 5
906 V[2] = intersectionSegmentPlan(it4, 2);
909 V[5] = intersectionSegmentPlan(it4, 5);
913 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
914 GMplus[TETRA4].push_back(it4);
916 else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
917 GMplus[TETRA4].push_back(it4);
919 // -------------------------------------------------------------------
921 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
922 { // Cas 4, arêtes 1 2 3 4
924 V[1] = intersectionSegmentPlan(it4, 1);
925 V[2] = intersectionSegmentPlan(it4, 2);
926 V[3] = intersectionSegmentPlan(it4, 3);
927 V[4] = intersectionSegmentPlan(it4, 4);
932 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
933 { // Cas 2, arêtes 1 3
935 V[1] = intersectionSegmentPlan(it4, 1);
937 V[3] = intersectionSegmentPlan(it4, 3);
943 else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
944 { // Cas 3, arêtes 1 3 5
946 V[1] = intersectionSegmentPlan(it4, 1);
948 V[3] = intersectionSegmentPlan(it4, 3);
950 V[5] = intersectionSegmentPlan(it4, 5);
954 // -------------------------------------------------------------------
956 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
957 { // Cas 2, arêtes 2 4
960 V[2] = intersectionSegmentPlan(it4, 2);
962 V[4] = intersectionSegmentPlan(it4, 4);
967 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
968 GMplus[TETRA4].push_back(it4);
970 else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
971 GMplus[TETRA4].push_back(it4);
973 // -------------------------------------------------------------------
975 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
976 { // Cas 3, arêtes 2 4 5
979 V[2] = intersectionSegmentPlan(it4, 2);
981 V[4] = intersectionSegmentPlan(it4, 4);
982 V[5] = intersectionSegmentPlan(it4, 5);
986 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
987 GMplus[TETRA4].push_back(it4);
989 else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
990 GMplus[TETRA4].push_back(it4);
993 ERREUR("Case not taken into account");
996 cout << chrono() << " - End of iteration on tetra4" << endl;
998 // cout << "indexNouveauxNoeuds = " << indexNouveauxNoeuds << endl;
999 newXX.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1000 newYY.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1001 newZZ.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1003 if (cptNouvellesMailles[TETRA4])
1004 newCNX[TETRA4].resize(4 * cptNouvellesMailles[TETRA4]);
1005 if (cptNouvellesMailles[PYRAM5])
1006 newCNX[PYRAM5].resize(5 * cptNouvellesMailles[PYRAM5]);
1007 if (cptNouvellesMailles[PENTA6])
1008 newCNX[PENTA6].resize(6 * cptNouvellesMailles[PENTA6]);
1010 // =========================================================================================
1011 // 2. Constitution du maillage final
1012 // =========================================================================================
1014 cout << chrono() << " - Constitution of final mesh" << endl;
1016 MAILLAGE2 = new Maillage(str_id_maillagenew);
1017 MAILLAGE2->dimensionMaillage = MAILLAGE1->dimensionMaillage;
1018 MAILLAGE2->dimensionEspace = MAILLAGE1->dimensionEspace;
1019 strcpy(MAILLAGE2->axisname, MAILLAGE1->axisname);
1020 strcpy(MAILLAGE2->unitname, MAILLAGE1->unitname);
1021 MAILLAGE2->nombreNoeudsMaillage = indexNouveauxNoeuds;
1022 MAILLAGE2->nombreMaillesMaillage = MAILLAGE1->nombreMaillesMaillage + cptNouvellesMailles[TETRA4]
1023 + cptNouvellesMailles[PYRAM5] + cptNouvellesMailles[PENTA6];
1025 // ---------- Coordonnées
1026 // Optimisation de la mémoire au détriment du temps
1028 // Héritage des coordonnées MAILLAGE1
1029 MAILLAGE2->XX = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1030 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1031 *(MAILLAGE2->XX + i) = *(MAILLAGE1->XX + i);
1032 free(MAILLAGE1->XX);
1033 MAILLAGE2->YY = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1034 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1035 *(MAILLAGE2->YY + i) = *(MAILLAGE1->YY + i);
1036 free(MAILLAGE1->YY);
1037 MAILLAGE2->ZZ = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1038 for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1039 *(MAILLAGE2->ZZ + i) = *(MAILLAGE1->ZZ + i);
1040 free(MAILLAGE1->ZZ);
1042 // Coordonnées des noeuds créés
1043 for (int i = 0; i < MAILLAGE2->nombreNoeudsMaillage - MAILLAGE1->nombreNoeudsMaillage; i++)
1045 *(MAILLAGE2->XX + MAILLAGE1->nombreNoeudsMaillage + i) = newXX[i];
1046 *(MAILLAGE2->YY + MAILLAGE1->nombreNoeudsMaillage + i) = newYY[i];
1047 *(MAILLAGE2->ZZ + MAILLAGE1->nombreNoeudsMaillage + i) = newZZ[i];
1048 // cout << "Nouveaux noeuds, indice " << i << " : " << newXX[i] << " " << newYY[i] << " " << newZZ[i] << " " << endl;
1051 // Legacy mailles maillage 1 (volumes seulement)
1052 for (int itm = (int) TETRA4; itm <= (int) HEXA20; itm++)
1054 TYPE_MAILLE tm = (TYPE_MAILLE) itm;
1055 if (tm != TETRA4 && tm != PYRAM5 && tm != PENTA6)
1057 // Pour les types autres que TETRA4 PYRAM5 PENTA6 on fait seulement pointer CNX2 vers CNX1
1058 if (MAILLAGE1->EFFECTIFS_TYPES[tm])
1059 MAILLAGE2->CNX[tm] = MAILLAGE1->CNX[tm];
1060 MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm];
1064 // Pour les types TETRA4 PYRAM5 PENTA6 on recopie CNX1 et on ajoute à la suite les newCNX
1065 // cout << "Legacy " << tm << " effectif " << MAILLAGE1->EFFECTIFS_TYPES[tm] << endl;
1066 int tailleType = Nnoeuds(tm);
1068 MAILLAGE2->CNX[tm] = (int*) malloc(sizeof(int) * tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm]
1069 + cptNouvellesMailles[tm]));
1070 for (int i = 0; i < MAILLAGE1->EFFECTIFS_TYPES[tm]; i++)
1071 for (int j = 0; j < tailleType; j++)
1072 *(MAILLAGE2->CNX[tm] + tailleType * i + j) = *(MAILLAGE1->CNX[tm] + tailleType * i + j);
1074 for (int i = 0; i < cptNouvellesMailles[tm]; i++)
1075 for (int j = 0; j < tailleType; j++)
1076 *(MAILLAGE2->CNX[tm] + tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm] + i) + j) = newCNX[tm][i * tailleType
1079 MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm] + cptNouvellesMailles[tm];
1085 // cout << "Maillage 2 - CNX TETRA4 : " << endl;
1087 // for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[TETRA4]; i++)
1089 // cout << "Maille " << i << " : ";
1090 // for (int j = 0; j < 4; j++)
1091 // cout << MAILLAGE2->CNX[TETRA4][i * 4 + j] << " ";
1095 // cout << "Maillage 2 - CNX PENTA6 : " << endl;
1097 // for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[PENTA6]; i++)
1099 // cout << "Maille " << i << " : ";
1100 // for (int j = 0; j < 6; j++)
1101 // cout << MAILLAGE2->CNX[PENTA6][i * 6 + j] << " ";
1106 // Groupes de mailles
1107 // MAILLAGE2->GM = MAILLAGE1->GM;
1108 MAILLAGE2->GN.clear();
1109 MAILLAGE2->GM.clear();
1110 MAILLAGE2->GM[str_id_GMplus] = GMplus;
1111 MAILLAGE2->GM[str_id_GMmoins] = GMmoins;
1113 // MAILLAGE2->GN = MAILLAGE1->GN;
1115 MAILLAGE2->eliminationMailles(TETRA4, cutTetras);
1117 cout << chrono() << " - MED file writing" << endl;
1119 MAILLAGE2->outputMED(ficMEDout);
1120 cout << chrono() << " - Finished!" << endl << endl;