/* Copyright 2012-2017 Inria ** This file is part of the PaMPA software package for parallel ** mesh partitioning and adaptation. ** ** PaMPA is free software: you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation, either version 3 of the License, or ** any later version. ** ** PaMPA is distributed in the hope that it will be useful, ** but WITHOUT ANY WARRANTY; without even the implied warranty of ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ** GNU General Public License for more details. ** ** In this respect, the user's attention is drawn to the risks associated ** with loading, using, modifying and/or developing or reproducing the ** software by the user in light of its specific status of free software, ** that may mean that it is complicated to manipulate, and that also ** therefore means that it is reserved for developers and experienced ** professionals having in-depth computer knowledge. Users are therefore ** encouraged to load and test the software's suitability as regards ** their requirements in conditions enabling the security of their ** systems and/or data to be ensured and, more generally, to use and ** operate it in the same conditions as regards security. ** ** The fact that you are presently reading this means that you have had ** knowledge of the GPLv3 license and that you accept its terms. */ /************************************************************/ //! //! \file dmesh_adapt_part.c //! //! \authors Cedric Lachat //! //! //! \date Version 1.0: from: 6 Dec 2012 //! to: 22 Sep 2017 //! /************************************************************/ #define DMESH #define DMESH_ADAPT #ifdef PAMPA_TIME_CHECK2 #define PAMPA_TIME_CHECK #endif /* PAMPA_TIME_CHECK2 */ //#define PAMPA_DEBUG_GRAPH //#define PAMPA_DEBUG_DMESH_SAVE #include "module.h" #include "common.h" #include "comm.h" #include "dvalues.h" #include "dmesh.h" #include "values.h" #include "mesh.h" #include "smesh.h" #include "pampa.h" #include "pampa.h" #include #include #include "dmesh_adapt.h" #include "dmesh_adapt_common.h" #include "dmesh_dgraph.h" #include "dmesh_gather_induce_multiple.h" #include "dmesh_rebuild.h" #include "dmesh_allreduce.h" #include "dmesh_zone_dist.h" DMESHALLREDUCEMINMAXOP (3, 3) typedef struct VertLoad_ { Gnum vertnum; Gnum loadval; Gnum edgenum; } VertLoad; /** XXX ** It returns: ** - 0 : on success. ** - !0 : on error. */ int dmeshAdaptPart ( Dmesh * const srcmeshptr, /**< mesh */ PAMPA_AdaptInfo * const infoptr, Dmesh * const dstmeshptr) { #ifdef PAMPA_DEBUG_ADAPT_SAVE3 MPI_Comm comm1p; #endif /* PAMPA_DEBUG_ADAPT_SAVE3 */ Dmesh indmeshdat; Dmesh * orgmeshptr; SCOTCH_Dgraph dgrfdat; SCOTCH_Dgraph bgrfdat; /* Band graph */ SCOTCH_Strat strtdat; Gnum * restrict rmshloctax; Gnum rmshlocnbr; Gnum rmshglbnbr; Gnum rmshlocnb2; Gnum rmshglbnb2; Gnum enodlocnbr; Gnum enodglbnbr; //Gnum rmshglbnb2; //Gnum rmshglbnb3; Gnum bvellocsum; Gnum bvelglbsum; Gnum baseval; Gnum vertlocmax; Gnum coarval; Gnum coarnbr; double coarrat; Gnum ballsi2; Gnum vertlocnum; Gnum vertlocnnd; Gnum dvrtlocnbr; Gnum emraval; Gnum switval; Gnum cvrtglbnbr; Gnum ibndval; int cheklocval; Gnum * restrict paroloctax; Gnum * restrict vmloloctab; Gnum * restrict cprcvrttab; Gnum * restrict srcprocvrttab; MPI_Comm proccomm; #ifdef PAMPA_DEBUG_DMESH_SAVE static Gnum wmshcnt = 0; #endif /* PAMPA_DEBUG_DMESH_SAVE */ //rmshglbnb2 = // rmshglbnb3 = -1; #ifdef PAMPA_INFO_ADAPT infoPrint ("*******************\n\nPart\n\n*******************"); #endif /* PAMPA_INFO_ADAPT */ orgmeshptr = srcmeshptr; cheklocval = 0; coarnbr = 1; coarrat = 1; switval = 1; for (coarval = 1, ballsi2 = 1; ballsi2 < infoptr->ballsiz; coarval ++, ballsi2 <<= 1); baseval = orgmeshptr->baseval; proccomm = orgmeshptr->proccomm; #ifdef PAMPA_DEBUG_ADAPT_SAVE3 MPI_Comm_split(proccomm, orgmeshptr->proclocnum, 0, &comm1p); #endif /* PAMPA_DEBUG_ADAPT_SAVE3 */ CHECK_FERR (dmeshInit (&indmeshdat, proccomm), proccomm); CHECK_FERR (SCOTCH_stratInit (&strtdat), proccomm); CHECK_FERR (dmeshDgraphInit (orgmeshptr, &dgrfdat), proccomm); CHECK_FERR (SCOTCH_dgraphInit (&bgrfdat, orgmeshptr->proccomm), proccomm); CHECK_VERR (cheklocval, proccomm); infoptr->iitenum = 1; if (memAllocGroup ((void **) (void *) &cprcvrttab, (size_t) ((orgmeshptr->procglbnbr + 1) * sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); cvrtglbnbr = 0; ibndval = infoptr->ibndval; while (1) { //while (infoptr->iitenum < 10) { Gnum dvrtlocbas; Gnum bvrtlocnbr; Gnum bvrtglbnbr; Gnum bvrtlocnum; Gnum bvrtlocnnd; Gnum bvrtgstnbr; Gnum bvrtgstmax; Gnum cvrtgstnbr; Gnum cvrtlocnbr; Gnum cvrtlocnum; Gnum cvrtlocnnd; Gnum edgelocnum; Gnum partlocidx; Gnum meshlocnbr; Gnum meshlocnum; Gnum fronlocidx; Gnum nodeent; int procngbnum; Gnum * restrict srcpartloctax; Gnum * restrict bvlbloctax; Gnum * restrict rmshloctx2; Gnum * restrict bvelloctax; Gnum * restrict dvelloctax; double * restrict dvelloctx2; Gnum * restrict vnumgsttax; Gnum * restrict bnumgsttax; Gnum * restrict bvrtloctax; Gnum * restrict bvndloctax; Gnum * restrict bedggsttax; Gnum * restrict cvelloctab; Gnum * restrict parttax; VertLoad * restrict veloglbtax; Mesh * imshloctab; Mesh * omshloctab; nodeent = orgmeshptr->enttglbnbr - (1 - baseval); // XXX temporaire, récupérer la valeur de data passé en argument CHECK_FERR (dmeshValueData (orgmeshptr, baseval, PAMPA_TAG_REMESH, NULL, NULL, (void **) &rmshloctax), proccomm); rmshloctax -= baseval; CHECK_FERR (PAMPA_dmeshHaloValue ((PAMPA_Dmesh *) orgmeshptr, baseval, PAMPA_TAG_REMESH), proccomm); // XXX début temporaire //#define PAMPA_DEBUG_DMESH_SAVE // XXX fin temporaire #ifdef PAMPA_DEBUG_DMESH_SAVE //// XXX début temporaire // if (infoptr->eitenum == 2) //// XXX fin temporaire { char s1[50]; if (memAllocGroup ((void **) (void *) &vnumgsttax, (size_t) (orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)), &parttax, (size_t) (sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); vnumgsttax -= baseval; parttax -= baseval; memSet (vnumgsttax + baseval, ~0, orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)); for (vertlocnum = baseval, vertlocnnd = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) if (rmshloctax[vertlocnum] == PAMPA_TAG_VERT_REMESH) vnumgsttax[vertlocnum] = 0; //for (procngbnum = 0; procngbnum < orgmeshptr->procglbnbr; procngbnum ++) parttax[baseval] = 0; imshloctab = NULL; meshlocnbr = -1; // FIXME l'allocation de imshloctab ne devrait pas être faite ici ? CHECK_FERR (dmeshGatherInduceMultiple (orgmeshptr, 0, 1, vnumgsttax, parttax + baseval, &meshlocnbr, &imshloctab), proccomm); if (meshlocnbr > 0) { sprintf (s1, "dmesh-adapt-rmsh1-%d-%d.mesh", infoptr->eitenum, infoptr->iitenum); CHECK_FERR (infoptr->EXT_meshSave((PAMPA_Mesh *) (imshloctab), 1, infoptr->dataptr, NULL, s1), proccomm); meshExit (imshloctab); } memFree (imshloctab); memFreeGroup (vnumgsttax + baseval); } #endif /* PAMPA_DEBUG_DMESH_SAVE */ // XXX début temporaire //#undef PAMPA_DEBUG_DMESH_SAVE // XXX fin temporaire vertlocmax = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval - 1; for (rmshlocnbr = 0, vertlocnum = baseval; vertlocnum <= vertlocmax; vertlocnum ++) { if (rmshloctax[vertlocnum] == PAMPA_TAG_VERT_REMESH) rmshlocnbr ++; } if (commAllreduce (&rmshlocnbr, &rmshglbnbr, 1, GNUM_MPI, MPI_SUM, proccomm) != MPI_SUCCESS) { errorPrint ("communication error"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); //#ifndef PAMPA_INFO_ADAPT if (orgmeshptr->proclocnum == 0) { //#endif /* PAMPA_INFO_ADAPT */ infoPrint ("-- Iteration number: %d", infoptr->iitenum); infoPrint ("Number of elements to be remeshed before growing band: %d", rmshglbnbr); //#ifndef PAMPA_INFO_ADAPT } //#endif /* PAMPA_INFO_ADAPT */ //rmshglbnb3 = rmshglbnb2; //rmshglbnb2 = rmshglbnbr; //if ((rmshglbnbr == 0) || (infoptr->iitenum > 2)) {//CHANGER NB ITER ICI --- CECILE if (rmshglbnbr == 0) { CHECK_FERR (dmeshDgraphFree (&dgrfdat), proccomm); break; } // XXX ne devrait-on pas faire cette alloc une bonne fois pour toutes les // itérations ? if (memAllocGroup ((void **) (void *) &rmshloctx2, (size_t) (orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)), &dvelloctax, (size_t) (orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); rmshloctx2 -= baseval; dvelloctax -= baseval; CHECK_FERR (dmeshValueData (orgmeshptr, baseval, PAMPA_TAG_WEIGHT, NULL, NULL, (void **) &dvelloctx2), proccomm); dvelloctx2 -= baseval; memCpy (rmshloctx2 + baseval, rmshloctax + baseval, orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)); //if (infoptr->iitenum > 1) CHECK_FERR (infoptr->EXT_dmeshBand ((PAMPA_Dmesh *) orgmeshptr, infoptr->dataptr, rmshloctx2 + baseval, ibndval), proccomm); for (rmshlocnb2 = 0, vertlocnum = baseval; vertlocnum <= vertlocmax; vertlocnum ++) //rmshlocnb2 += ceil (dvelloctx2[vertlocnum]); if (rmshloctx2[vertlocnum] == PAMPA_TAG_VERT_REMESH) rmshlocnb2 ++; if (commAllreduce (&rmshlocnb2, &rmshglbnb2, 1, GNUM_MPI, MPI_SUM, proccomm) != MPI_SUCCESS) { errorPrint ("communication error"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); ////printf("loc %d\n",enodlocnbr); //CHECK_VDBG (cheklocval, proccomm); //if (commAllreduce (&enodlocnbr, &enodglbnbr, 1, GNUM_MPI, MPI_SUM, proccomm) != MPI_SUCCESS) { // errorPrint ("communication error"); // cheklocval = 1; //} //CHECK_VERR (cheklocval, proccomm); //if (orgmeshptr->proclocnum == 0) //printf("npcible %d ntetra %d\n",enodglbnbr, rmshglbnbr); //cvrtglbnbr = ((rmshglbnbr / 6) > 10000) ? MAX((enodglbnbr - orgmeshptr->enttloctax[baseval]->vertglbnbr / 6) / 300000, rmshglbnbr / 6000000): 1; double cond, val1, val2; cond = 10000; //cond = 100; if (infoptr->iitenum > 2) if ((infoptr->enodglbnbr / (orgmeshptr->enttloctax[baseval]->vertglbnbr / 6)) > 2) infoptr->enodglbnbr = MAX (orgmeshptr->enttloctax[baseval]->vertglbnbr / 3, infoptr->enodglbnbr / 2); else if ((1.1 * infoptr->enodglbnbr / (orgmeshptr->enttloctax[baseval]->vertglbnbr / 6)) > 1) infoptr->enodglbnbr = ((1.1 * orgmeshptr->enttloctax[baseval]->vertglbnbr / 6) > (0.8 * infoptr->enodglbnbr)) ? 1.1 * orgmeshptr->enttloctax[baseval]->vertglbnbr / 6 : 0.8 * infoptr->enodglbnbr; val1 = (1.0 * infoptr->enodglbnbr - (1.0 * orgmeshptr->enttloctax[baseval]->vertglbnbr - rmshglbnb2) / 6) / 250000; //val1 = (1.0 * infoptr->enodglbnbr - (1.0 * orgmeshptr->enttloctax[baseval]->vertglbnbr - rmshglbnb2) / 6) / 2500; val2 = 1.0 * rmshglbnb2 / 1000000; cvrtglbnbr = ((1.0 * rmshglbnb2 / 6) > cond) ? MAX(ceil(val1), ceil(val2)) : 1; // XXX tmp //if (infoptr->iitenum == 1) //cvrtglbnbr = 10 / infoptr->iitenum; //cvrtglbnbr = 2; // XXX fin tmp //if ((1.0 * rmshglbnb2 / cvrtglbnbr) < 10000) // XXX test à remettre pour éviter qu'il y ait moins de 10000 éléments par zone // cvrtglbnbr /= ceil (1.0 * cvrtglbnbr * 10000 / rmshglbnb2); //} // XXX temporaire //if (infoptr->iitenum == 1) // cvrtglbnbr = 4; //if (infoptr->iitenum == 2) // cvrtglbnbr = coldglbnbr; //else if (infoptr->iitenum == 3) // cvrtglbnbr = 1; // XXX fin temporaire //if ((coldglbnbr == cvrtglbnbr) && (ibndval != 0) && (cvrtglbnbr != 1)) { // ibndval = MAX(0, ibndval - 1); // CHECK_FERR (dmeshDgraphFree (&dgrfdat), proccomm); // memFree (hashtab); // continue; //} //else if ((coldglbnbr == cvrtglbnbr) && (cvrtglbnbr == 1)) { // CHECK_FERR (dmeshDgraphFree (&dgrfdat), proccomm); // memFree (hashtab); // break; //} if (infoptr->zonenb2 != ~0) cvrtglbnbr = infoptr->zonenb2; if (infoptr->iitenum >= 4) { cvrtglbnbr /= 2; infoptr->zonenb2 = cvrtglbnbr; } //cvrtglbnbr = ceil (1.0 * rmshglbnb2 / 25000); if (orgmeshptr->proclocnum == 0) infoPrint ("vertglb: " GNUMSTRING ", rmshglb: " GNUMSTRING ", rmshglb2: " GNUMSTRING ", enodglbnbr: " GNUMSTRING ", cond: %lf, val1: %lf, val2: %lf, cvrt:" GNUMSTRING "\n", orgmeshptr->enttloctax[baseval]->vertglbnbr, rmshglbnbr, rmshglbnb2, infoptr->enodglbnbr, cond, val1, val2, cvrtglbnbr); #ifdef PAMPA_DEBUG_ADAPT if (0 != PAMPA_TAG_VERT_NOREMESH) { errorPrint ("PAMPA_TAG_VERT_NOREMESH is not 0"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); #endif /* PAMPA_DEBUG_ADAPT */ #ifdef PAMPA_DEBUG_DMESH_SAVE { char s1[50]; if (memAllocGroup ((void **) (void *) &vnumgsttax, (size_t) (orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)), &parttax, (size_t) (sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); vnumgsttax -= baseval; parttax -= baseval; memSet (vnumgsttax + baseval, ~0, orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)); for (vertlocnum = baseval, vertlocnnd = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) if (rmshloctx2[vertlocnum] == PAMPA_TAG_VERT_REMESH) vnumgsttax[vertlocnum] = 0; //for (procngbnum = 0; procngbnum < orgmeshptr->procglbnbr; procngbnum ++) parttax[baseval] = 0; imshloctab = NULL; meshlocnbr = -1; // FIXME l'allocation de imshloctab ne devrait pas être faite ici ? CHECK_FERR (dmeshGatherInduceMultiple (orgmeshptr, 0, 1, vnumgsttax, parttax + baseval, &meshlocnbr, &imshloctab), proccomm); if (meshlocnbr > 0) { sprintf (s1, "dmesh-adapt-rmsh2-%d-%d.mesh", infoptr->eitenum, infoptr->iitenum); CHECK_FERR (infoptr->EXT_meshSave((PAMPA_Mesh *) (imshloctab), 1, infoptr->dataptr, NULL, s1), proccomm); meshExit (imshloctab); } memFree (imshloctab); memFreeGroup (vnumgsttax + baseval); } #endif /* PAMPA_DEBUG_DMESH_SAVE */ // XXX début temporaire //#undef PAMPA_DEBUG_DMESH_SAVE // XXX fin temporaire for (vertlocnum = baseval; vertlocnum <= vertlocmax; vertlocnum ++) { dvelloctax[vertlocnum] = round(dvelloctx2[vertlocnum]); } // 1) Identification des zones à remailler // a) conversion du maillage en graphe uniquement pour les sommets qui ont un drapeau // XXX FIXME bonne idée d'utiliser DMESH_ENTT_MAIN_PART, cela évite de // calculer ensuite le graphe induit, mais faut-il encore avoir // l'ancienne numérotation (vlblloctax à rajouter dans dmeshDgraphBuild // pour DMESH_ENTT_MAIN_PART). En attendant, on continue comme avant... //CHECK_FERR (dmeshDgraphBuild (orgmeshptr, DMESH_ENTT_MAIN_PART, rmshlocnbr, dvrtgsttax, dvelloctax, &dgrfdat, &dvrtlocbas), proccomm); //CHECK_FERR (dmeshDgraphBuild (orgmeshptr, DMESH_ENTT_MAIN, ~0, NULL, dvelloctax, NULL, &dgrfdat, &dvrtlocbas), proccomm); // XXX pas de edlo ?? // XXX temporaire pour éviter l'effet météorite dû aux poids sur les éléments CHECK_FERR (dmeshDgraphBuild (orgmeshptr, DMESH_ENTT_MAIN, ~0, NULL, NULL, NULL, &dgrfdat, &dvrtlocbas), proccomm); // XXX pas de edlo ?? // XXX fin temporaire CHECK_FERR (SCOTCH_dgraphGhst (&dgrfdat), proccomm); #ifdef PAMPA_DEBUG_DGRAPH char s[30]; sprintf(s, "dgrfdat-adapt-%d", orgmeshptr->proclocnum); FILE *dgraph_file; dgraph_file = fopen(s, "w"); SCOTCH_dgraphSave (&dgrfdat, dgraph_file); fclose(dgraph_file); #endif /* PAMPA_DEBUG_DGRAPH */ //CHECK_VDBG (cheklocval, proccomm); //if (commAllreduce (&rmshlocnbr, &rmshglbnbr, 1, GNUM_MPI, MPI_SUM, proccomm) != MPI_SUCCESS) { // errorPrint ("communication error"); // cheklocval = 1; //} //CHECK_VERR (cheklocval, proccomm); // a bis) élargissement des zones avec un dgraphBand ok //if (ibndval > 0) { // CHECK_FERR (SCOTCH_dgraphBand (&dgrfdat, rmshlocnbr, fronloctax + baseval, ibndval, &bgrfdat), proccomm); //} //else { // b) sous-graphe induit des zones ok memCpy (rmshloctx2 + baseval, rmshloctax + baseval, orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)); CHECK_FERR (SCOTCH_dgraphInducePart (&dgrfdat, rmshloctx2 + baseval, PAMPA_TAG_VERT_REMESH, rmshlocnbr, &bgrfdat), proccomm); //} #ifdef PAMPA_DEBUG_ADAPT2 CHECK_FERR (SCOTCH_dgraphCheck (&bgrfdat), proccomm); #endif /* PAMPA_DEBUG_ADAPT2 */ CHECK_FERR (SCOTCH_dgraphGhst (&bgrfdat), proccomm); //SCOTCH_dgraphData (&bgrfdat, NULL, NULL, &bvrtlocnbr, NULL, &bvrtgstnbr, (SCOTCH_Num **) &bvrtloctax, (SCOTCH_Num **) &bvndloctax, NULL, (SCOTCH_Num **) NULL, NULL, NULL, NULL, NULL, (SCOTCH_Num **) &bedggsttax, NULL, NULL); SCOTCH_dgraphData (&bgrfdat, NULL, &bvrtglbnbr, &bvrtlocnbr, NULL, &bvrtgstnbr, (SCOTCH_Num **) &bvrtloctax, (SCOTCH_Num **) &bvndloctax, (SCOTCH_Num **) &bvelloctax, (SCOTCH_Num **) &bvlbloctax, NULL, NULL, NULL, NULL, (SCOTCH_Num **) &bedggsttax, NULL, NULL); bvlbloctax -= baseval; bvelloctax -= baseval; bvrtloctax -= baseval; bvndloctax -= baseval; bedggsttax -= baseval; for (bvellocsum = bvrtlocnum = baseval, bvrtlocnnd = bvrtlocnbr + baseval; bvrtlocnum < bvrtlocnnd; bvrtlocnum ++) bvellocsum += bvelloctax[bvrtlocnum]; CHECK_VDBG (cheklocval, proccomm); if (commAllreduce (&bvellocsum, &bvelglbsum, 1, GNUM_MPI, MPI_SUM, proccomm) != MPI_SUCCESS) { errorPrint ("communication error"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); //if (switval == 1) // cvrtglbnbr = (bvelglbsum > ballsiz) ? bvelglbsum / ballsiz : 1; //else // cvrtglbnbr = (bvelglbsum > (2 * ballsiz)) ? bvelglbsum / (2 * ballsiz) : 1; //cvrtglbnbr = (bvelglbsum > ballsiz) ? bvelglbsum / ballsiz : 1; switval = 1 - switval; #ifdef PAMPA_INFO_ADAPT //infoPrint ("dvelglbsum : %d, rmshglbnbr : %d, ballsiz : %d, cvrtglbnbr : %d", bvelglbsum, rmshglbnbr, ballsiz, cvrtglbnbr); #else /* PAMPA_INFO_ADAPT */ if (orgmeshptr->proclocnum == 0) { infoPrint ("estimated weight: %d, number of elements to be remeshed: %d, number of zones: %d", bvelglbsum, rmshglbnbr, cvrtglbnbr); } #endif /* PAMPA_INFO_ADAPT */ // c) contraction du graphe des zones ok if (memAllocGroup ((void **) (void *) &bnumgsttax, (size_t) (bvrtgstnbr * sizeof (Gnum)), // XXX supprimer ci-dessous « * 2 » &cvelloctab, (size_t) (cvrtglbnbr * sizeof (Gnum)), //&bvlbloctax, (size_t) (bvrtlocnbr * sizeof (Gnum)), &vnumgsttax, (size_t) (orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)), &parttax, (size_t) ((cvrtglbnbr + orgmeshptr->procglbnbr) * sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); bnumgsttax -= baseval; vnumgsttax -= baseval; parttax -= baseval; //bvlbloctax -= baseval; memSet (cvelloctab, 0, cvrtglbnbr * sizeof (Gnum)); memSet (vnumgsttax + baseval, ~0, orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)); //if ((bvelglbsum > ballsiz) && (rmshglbnb2 != rmshglbnb3)) { if (cvrtglbnbr > 1) { //CHECK_FERR (dmeshAdapt2 (orgmeshptr, &bgrfdat, bnumgsttax, &cgrfdat, cprcvrttab, coarval - 1, coarnbr, coarrat), orgmeshptr->proccomm); #ifdef PAMPA_DEBUG_ADAPT2 CHECK_FERR (SCOTCH_dgraphCheck (&bgrfdat), proccomm); #endif /* PAMPA_DEBUG_ADAPT2 */ //#define PAMPA_DEBUG_DGRAPH #ifdef PAMPA_DEBUG_DGRAPH { char s[30]; FILE *dgraph_file; errorPrint ("cvrtglbnbr: %d\n", cvrtglbnbr); sprintf(s, "dgrf-adapt-%d-%d.grf", infoptr->iitenum, orgmeshptr->proclocnum); dgraph_file = fopen(s, "w"); SCOTCH_dgraphSave (&bgrfdat, dgraph_file); fclose(dgraph_file); Gnum tetentt, nodentt, i; PAMPA_Iterator it_nghb; double * coorloctax; tetentt = 0; nodentt = 1; CHECK_FERR (PAMPA_dmeshHaloValue ((PAMPA_Dmesh *) orgmeshptr, nodentt, PAMPA_TAG_GEOM), proccomm); dmeshValueData(orgmeshptr, nodentt, PAMPA_TAG_GEOM, NULL, NULL, (void **) &coorloctax); coorloctax -= 3 * baseval; sprintf(s, "dgrf-adapt-%d-%d.xyz", infoptr->iitenum, orgmeshptr->proclocnum); dgraph_file = fopen(s, "w"); if (orgmeshptr->proclocnum == 0) { fprintf (dgraph_file, "3\n%d\n", orgmeshptr->enttloctax[baseval]->vertglbnbr); } //if (commAllgather (&bvrtlocnbr, 1, GNUM_MPI, // vnbrloctab, 1, GNUM_MPI, orgmeshptr->proccomm) != MPI_SUCCESS) { // errorPrint ("communication error"); // return (1); //} //vertlocbas = 0; //for (procngbnum = 0; procngbnum < orgmeshptr->proclocnum ; procngbnum ++) // vertlocbas += vnbrloctab[procngbnum]; PAMPA_dmeshItInit((PAMPA_Dmesh *) orgmeshptr, tetentt, nodentt, &it_nghb); for (vertlocnum = baseval, vertlocnnd = bvrtlocnbr; vertlocnum < vertlocnnd; vertlocnum ++) { PAMPA_Num tetnum; double coorval[3]; tetnum = bvlbloctax[vertlocnum] - dvrtlocbas; coorval[0] = coorval[1] = coorval[2] = 0.0; PAMPA_itStart(&it_nghb, tetnum); while (PAMPA_itHasMore(&it_nghb)) { PAMPA_Num nodenum; nodenum = PAMPA_itCurEnttVertNum(&it_nghb); for (i = 0; i < 3; i++) coorval[i] += coorloctax[nodenum * 3 + i]; PAMPA_itNext(&it_nghb); } //fprintf(dgraph_file, "%d", vertlocbas + tetnum); fprintf(dgraph_file, "%d", vertlocnum); for (i = 0; i < 3; i++) fprintf(dgraph_file, " %lf", coorval[i] / 4); fprintf(dgraph_file, "\n"); } fclose (dgraph_file); } #endif /* PAMPA_DEBUG_DGRAPH */ //#undef PAMPA_DEBUG_DGRAPH //CHECK_FERR (SCOTCH_stratDgraphMap (&strtdat, "r{sep=q{strat=g},seq=r{sep=g}}"), orgmeshptr->proccomm); CHECK_FERR (SCOTCH_dgraphPart (&bgrfdat, cvrtglbnbr, &strtdat, bnumgsttax + baseval), orgmeshptr->proccomm); // XXX TODO vérifier que chaque zone ait le nombre minimum d'éléments // émettre un avertissement si moins que prévu #ifdef PAMPA_DEBUG_ADAPT_PART { if (orgmeshptr->proclocnum == 0) { SCOTCH_Graph * graflvl; SCOTCH_Graph * grafptr; int tmpval; Gnum vertnbr; Gnum vertnum; Gnum vertnnd; Gnum colrnum; Gnum colrflg; Gnum * restrict finematetax; Gnum * restrict coarmulttax; int * restrict drcvcnttab; int * restrict drcvdsptab; if (memAllocGroup ((void **) (void *) &drcvcnttab, (size_t) (orgmeshptr->procglbnbr * sizeof (int)), &drcvdsptab, (size_t) ((orgmeshptr->procglbnbr + 1) * sizeof (int)), &parttax, (size_t) (bvrtglbnbr * sizeof (Gnum)), &parttx2, (size_t) (MAX(cvrtglbnbr, bvrtglbnbr) * sizeof (Gnum)), &finematetax, (size_t) (bvrtglbnbr * sizeof (Gnum)), &coarmulttax, (size_t) (bvrtglbnbr * 2 * sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } parttax -= baseval; parttx2 -= baseval; finematetax -= baseval; coarmulttax -= 2 * baseval; if ((grafptr = (SCOTCH_Graph *) memAlloc (sizeof (SCOTCH_Graph))) == NULL) { errorPrint ("out of memory"); cheklocval = 1; } CHECK_VDBG2 (cheklocval); CHECK_FDBG2 (SCOTCH_graphInit (grafptr)); CHECK_FERR (SCOTCH_dgraphGather (&bgrfdat, grafptr), orgmeshptr->proccomm); tmpval = bvrtlocnbr; CHECK_FMPI (cheklocval, commGather (&tmpval, 1, MPI_INT, drcvcnttab, 1, MPI_INT, 0, proccomm), proccomm); drcvdsptab[0] = 0; for (procngbnum = 0; procngbnum < orgmeshptr->procglbnbr; procngbnum ++) { drcvdsptab[procngbnum + 1] = drcvdsptab[procngbnum] + drcvcnttab[procngbnum]; } CHECK_FMPI (cheklocval, commGatherv (bnumgsttax + baseval, bvrtlocnbr, GNUM_MPI, parttax + baseval, drcvcnttab, drcvdsptab, GNUM_MPI, 0, proccomm), proccomm); do { Gnum coarvertnbr; Gnum * restrict verttax; Gnum * restrict vendtax; Gnum * restrict edgetax; if ((graflvl = (SCOTCH_Graph *) memAlloc (sizeof (SCOTCH_Graph))) == NULL) { errorPrint ("out of memory"); cheklocval = 1; } CHECK_VDBG2 (cheklocval); CHECK_FDBG2 (SCOTCH_graphInit (graflvl)); SCOTCH_graphData (grafptr, NULL, &vertnbr, &verttax, &vendtax, NULL, NULL, NULL, &edgetax, NULL); verttax -= baseval; vendtax -= baseval; edgetax -= baseval; memSet (finematetax + baseval, ~0, vertnbr * sizeof (Gnum)); for (colrflg = 0, coarvertnbr = vertnbr, vertnum = baseval, vertnnd = vertnbr + baseval; vertnum < vertnnd; vertnum ++) { Gnum edgenum; Gnum edgennd; if (finematetax[vertnum] == ~0) for (edgenum = verttax[vertnum], edgennd = vendtax[vertnum]; edgenum < edgennd; edgenum ++) if ((finematetax[edgetax[edgenum]] == ~0) && (parttax[vertnum] == parttax[edgetax[edgenum]])) { if (colrflg == 0) colrflg = 1; finematetax[vertnum] = edgetax[edgenum]; finematetax[edgetax[edgenum]] = vertnum; coarvertnbr --; break; } } for (vertnum = baseval, vertnnd = vertnbr + baseval; vertnum < vertnbr; vertnum ++) /* If there are still alone vertices */ if (finematetax[vertnum] == ~0) finematetax[vertnum] = vertnum; CHECK_FDBG2 (SCOTCH_graphCoarsenBuild (grafptr, coarvertnbr, finematetax + baseval, graflvl, coarmulttax + 2 * baseval)); for (vertnum = baseval, vertnnd = coarvertnbr + baseval; vertnum < vertnnd; vertnum ++) parttx2[vertnum] = parttax[coarmulttax[vertnum << 1]]; memCpy (parttax + baseval, parttx2 + baseval, coarvertnbr * sizeof (Gnum)); SCOTCH_graphExit (grafptr); memFree (grafptr); grafptr = graflvl; } while (colrflg == 1); SCOTCH_graphSize (grafptr, &vertnbr, NULL); memSet (parttx2 + baseval, 0, cvrtglbnbr * sizeof (Gnum)); parttx2 += baseval; for (vertnum = baseval, vertnnd = baseval + vertnbr; vertnum < vertnnd; vertnum ++) parttx2[parttax[vertnum]] ++; for (colrflg = colrnum = 0; (colrnum < cvrtglbnbr) && (colrflg < 5); colrnum ++) // TRICK: start from 0 because it is a color if (parttx2[colrnum] != 1) { colrflg ++; errorPrint (GNUMSTRING "CC for color " GNUMSTRING, parttx2[colrnum], colrnum); } if (colrflg != 0) cheklocval = 1; SCOTCH_graphExit (grafptr); memFree (grafptr); memFreeGroup (drcvcnttab); } else { /* proclocnum != 0) */ int tmpval; CHECK_FERR (SCOTCH_dgraphGather (&bgrfdat, NULL), orgmeshptr->proccomm); tmpval = bvrtlocnbr; CHECK_FMPI (cheklocval, commGather (&tmpval, 1, MPI_INT, NULL, 0, MPI_INT, 0, proccomm), proccomm); CHECK_FMPI (cheklocval, commGatherv (bnumgsttax + baseval, bvrtlocnbr, GNUM_MPI, NULL, NULL, NULL, GNUM_MPI, 0, proccomm), proccomm); } } CHECK_VERR (cheklocval, proccomm); #endif /* PAMPA_DEBUG_ADAPT_PART */ cprcvrttab[0] = 0; for (procngbnum = 0; procngbnum < orgmeshptr->procglbnbr; procngbnum ++) cprcvrttab[procngbnum + 1] = cprcvrttab[procngbnum] + DATASIZE(cvrtglbnbr, orgmeshptr->procglbnbr, procngbnum); cvrtlocnbr = cprcvrttab[orgmeshptr->proclocnum + 1] - cprcvrttab[orgmeshptr->proclocnum]; //#define PAMPA_DEBUG_DMESH_SAVE2 #ifdef PAMPA_DEBUG_DMESH_SAVE2 { char s1[100]; Gnum * vnumloctax; Gnum * vflgloctax; Gnum * vflgglbtax; Gnum parttab[1]; parttab[0] = 0; if (memAllocGroup ((void **) (void *) &vnumloctax, (size_t) (orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)), &vflgloctax, (size_t) (bvrtlocnbr * sizeof (Gnum)), &vflgglbtax, (size_t) (bvrtglbnbr * sizeof (Gnum)), NULL) == NULL) { errorPrint ("Out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); vnumloctax -= baseval; vflgloctax -= baseval; vflgglbtax -= baseval; memSet (vnumloctax + baseval, ~0, orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)); for (vertlocnum = baseval, vertlocnnd = bvrtlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) if (rmshloctx2[bvlbloctax[vertlocnum] - dvrtlocbas] == PAMPA_TAG_VERT_REMESH) { vnumloctax[bvlbloctax[vertlocnum] - dvrtlocbas] = 0; vflgloctax[vertlocnum] = bnumgsttax[vertlocnum]; } imshloctab = NULL; meshlocnbr = -1; // FIXME l'allocation de imshloctab ne devrait pas être faite ici ? CHECK_FERR (dmeshGatherInduceMultiple (orgmeshptr, 0, 1, vnumloctax, parttab - baseval, &meshlocnbr, &imshloctab), proccomm); sprintf (s1, "dmesh-adapt-part-%d.mesh", infoptr->iitenum); CHECK_FERR (dmeshAdaptGatherSave3 (orgmeshptr, bvrtlocnbr, vflgloctax, vflgglbtax), proccomm); if (meshlocnbr > 0) { CHECK_FERR (infoptr->EXT_meshSave((PAMPA_Mesh *) (imshloctab), 1, infoptr->dataptr, vflgglbtax + baseval, s1), proccomm); meshExit (imshloctab); } memFreeGroup (vnumloctax + baseval); memFree (imshloctab); } #endif /* PAMPA_DEBUG_DMESH_SAVE2 */ //#undef PAMPA_DEBUG_DMESH_SAVE2 // XXX pour enregistrer la partition, il faut ne garder que la partie // du maillage correspondand au graphe bgrfdat (cf. rmshloctx2) //#ifdef PAMPA_DEBUG_DMESH_SAVE // char s1[50]; // sprintf (s1, "dmesh-adapt-part-%d.mesh", wmshcnt); // CHECK_FERR (dmeshAdaptGatherSave2 (orgmeshptr, bnumgsttax, dataptr, s1), proccomm); //#endif /* PAMPA_DEBUG_DMESH_SAVE */ SCOTCH_stratExit (&strtdat); SCOTCH_stratInit (&strtdat); CHECK_FERR (SCOTCH_dgraphHalo (&bgrfdat, bnumgsttax + baseval, GNUM_MPI), orgmeshptr->proccomm); //CHECK_FERR (SCOTCH_dgraphGhst (&cgrfdat), proccomm); //SCOTCH_dgraphData (&cgrfdat, NULL, &cvrtglbnbr, &cvrtlocnbr, NULL, &cvrtgstnbr, NULL, NULL, (SCOTCH_Num **) &cvelloctax, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); //cvelloctax -= baseval; //#define PAMPA_INFO_ADAPT #ifdef PAMPA_INFO_ADAPT infoPrint ("-- Iteration number: %d", infoptr->iitenum); infoPrint ("Number of vertices to be remeshed : %d", rmshglbnbr); infoPrint ("Number of colors : %d", cvrtglbnbr); infoPrint ("Number of local colors : %d", cvrtlocnbr); #endif /* PAMPA_INFO_ADAPT */ #undef PAMPA_INFO_ADAPT bvrtgstmax = bvrtlocnbr + baseval - 1; // d) partitionnement des gros sommets ko // * partitionnement de graphes // * les gros sommets sont indépendants donc absence de relations SCOTCH_dgraphSize (&dgrfdat, NULL, &dvrtlocnbr, NULL, NULL); //for (bvrtlocnum = dvrtlocnum = baseval, dvrtlocnnd = dvrtlocnbr + baseval; dvrtlocnum < dvrtlocnnd; dvrtlocnum ++) // if (rmshloctx2[dvrtlocnum] == PAMPA_TAG_VERT_REMESH) // bvlbloctax[bvrtlocnum ++] = dvrtlocnum + dvrtlocbas; memSet (vnumgsttax + baseval, ~0, orgmeshptr->enttloctax[baseval]->vertlocnbr * sizeof (Gnum)); for (bvrtlocnum = baseval, bvrtlocnnd = bvrtlocnbr + baseval; bvrtlocnum < bvrtlocnnd; bvrtlocnum ++) { //if (cvelgsttax[bnumgsttax[bvrtlocnum]] > (ballsiz / 10)) //FIXME peut-on faire un test sur les tailles de boule, le //problème est que bnumgsttax contient un sommet qui peut ne pas //être dans les ghst //if ((bnumgsttax[bvrtlocnum] < 0) || (bnumgsttax[bvrtlocnum] >= cvrtglbnbr)) { // errorPrint ("color " GNUMSTRING " invalid, max is " GNUMSTRING, bnumgsttax[bvrtlocnum], cvrtglbnbr); // return (1); //} vnumgsttax[bvlbloctax[bvrtlocnum] - dvrtlocbas] = bnumgsttax[bvrtlocnum]; /* Propagate colors on elements */ // XXX temporaire YYY if (bnumgsttax[bvrtlocnum] >= 0) // XXX fin temporaire YYY cvelloctab[bnumgsttax[bvrtlocnum]] += bvelloctax[bvrtlocnum]; } #ifdef PAMPA_DEBUG_ADAPT for (vertlocnum = baseval, vertlocnnd = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) { if ((vnumgsttax[vertlocnum] < ~0) || (vnumgsttax[vertlocnum] >= cvrtglbnbr)) { errorPrint ("Color " GNUMSTRING " invalid for vertex " GNUMSTRING, vnumgsttax[vertlocnum] , vertlocnum); return (1); } } #endif /* PAMPA_DEBUG_ADAPT */ CHECK_FERR (dmeshGrow ( orgmeshptr, baseval, vnumgsttax, cvrtglbnbr, nodeent, ibndval), orgmeshptr->proccomm); // XXX on fait des suppositions sur les entités, il faut passer par la bib externe /////* First pass to remove elements which have neighbors with different color */ //for (bvrtlocnum = baseval, bvrtlocnnd = bvrtlocnbr + baseval; bvrtlocnum < bvrtlocnnd; bvrtlocnum ++) { // Gnum cvrtlocnum; // Gnum bedglocnum; // Gnum bedglocnnd; // Gnum ucollocval; /* If undo coloring */ // cvrtlocnum = bnumgsttax[bvrtlocnum]; // for (ucollocval = 0, bedglocnum = bvrtloctax[bvrtlocnum], bedglocnnd = bvndloctax[bvrtlocnum]; bedglocnum < bedglocnnd; bedglocnum ++) /* Undo coloring on elements which are at the boundary */ // if (bnumgsttax[bedggsttax[bedglocnum]] != cvrtlocnum) { /* If vertex and neighbor have different colors */ // ucollocval = 1; // if (bedggsttax[bedglocnum] <= bvrtgstmax) // vnumgsttax[bvlbloctax[bedggsttax[bedglocnum]] - dvrtlocbas] = ~0; // } // if (ucollocval == 1) { /* There are at least one neighbor with different color */ // vnumgsttax[bvlbloctax[bvrtlocnum] - dvrtlocbas] = ~0; // // FIXME il faudrait diminuer le cvelloctax pour cette couleur // } //} memFreeGroup (rmshloctx2 + baseval); //CHECK_FERR (dmeshHaloSync (orgmeshptr, orgmeshptr->enttloctax[baseval], vnumgsttax + baseval, GNUM_MPI), orgmeshptr->proccomm); //// FIXME temporarly disabled ///* Second pass to remove elements which are isolated */ //for (vertlocnum = baseval, vertlocnnd = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) { // Gnum ucollocval; /* If undo coloring */ // Gnum vnumlocval; // Gnum edgelocnnd; // const DmeshEntity * enttlocptr = orgmeshptr->enttloctax[baseval]; // vnumlocval = vnumgsttax[vertlocnum]; // if (vnumlocval != ~0) { // for (ucollocval = 0, edgelocnum = enttlocptr->nghbloctax[baseval]->vindloctax[vertlocnum].vertlocidx, // edgelocnnd = enttlocptr->nghbloctax[baseval]->vindloctax[vertlocnum].vendlocidx; // ucollocval == 0 && edgelocnum < edgelocnnd; edgelocnum ++) { // Gnum nghblocnum; // nghblocnum = orgmeshptr->edgeloctax[edgelocnum]; // if (vnumgsttax[nghblocnum] == vnumlocval) /* If vertex and neighbor have same colors */ // ucollocval = 1; // } // if (ucollocval == 0) { /* If no neighbor with same color */ // vnumgsttax[vertlocnum] = ~0; // // FIXME il faudrait diminuer le cvelloctax pour cette couleur // } // } //} CHECK_FERR (dmeshZoneDist (orgmeshptr, infoptr, bvrtlocnbr, bnumgsttax, cvrtlocnbr, cvrtglbnbr, cprcvrttab, cvelloctab, parttax), orgmeshptr->proccomm); SCOTCH_dgraphFree (&bgrfdat); } else { memFreeGroup (rmshloctx2 + baseval); #ifdef PAMPA_INFO_ADAPT infoPrint ("-- Iteration number: %d", infoptr->iitenum); infoPrint ("Number of vertices to be remeshed : %d", rmshglbnbr); infoPrint ("Last iteration"); #endif /* PAMPA_INFO_ADAPT */ memSet (bnumgsttax + baseval, 0, bvrtgstnbr * sizeof (Gnum)); for (bvrtlocnum = baseval, bvrtlocnnd = bvrtlocnbr + baseval; bvrtlocnum < bvrtlocnnd; bvrtlocnum ++) vnumgsttax[bvlbloctax[bvrtlocnum] - dvrtlocbas] = bnumgsttax[bvrtlocnum]; /* Propagate colors on elements */ CHECK_FERR (dmeshGrow ( orgmeshptr, baseval, vnumgsttax, cvrtglbnbr, nodeent, ibndval), orgmeshptr->proccomm); // XXX on fait des suppositions sur les entités, il faut passer par la bib externe if ((parttax = (Gnum *) memAlloc (1 * sizeof (Gnum))) == NULL) { errorPrint ("out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); parttax -= baseval; parttax[baseval] = 0; // FIXME pourquoi 0 ??? cvrtglbnbr = 1; } // 2) Extraction des sous-maillages et envoie sur chaque proc ok meshlocnbr = -1; imshloctab = NULL; CHECK_FERR_STR (dmeshGatherInduceMultiple2 (orgmeshptr, 1, infoptr, cvrtglbnbr, vnumgsttax, parttax, &meshlocnbr, &imshloctab), "extraction", proccomm); //CHECK_FERR (dmeshValueLink (orgmeshptr, (void **) &srcprocvrttab, PAMPA_VALUE_PRIVATE, NULL, NULL, GNUM_MPI, PAMPA_ENTT_VIRT_PROC, TAG_INT_PROCVRT), proccomm); //memCpy (srcprocvrttab, orgmeshptr->procvrttab, (orgmeshptr->procglbnbr + 1) * sizeof (Gnum)); CHECK_FERR (dmeshValueData (orgmeshptr, PAMPA_ENTT_VIRT_VERT, PAMPA_TAG_STATUS, NULL, NULL, (void **) &srcpartloctax), proccomm); srcpartloctax -= baseval; #ifdef PAMPA_INFO_ADAPT infoPrint ("dmeshInducePart, vertex number before: %d", orgmeshptr->vertlocnbr); #endif /* PAMPA_INFO_ADAPT */ //CHECK_FERR_STR (dmeshInducePart (orgmeshptr, srcpartloctax, PAMPA_TAG_VERT_INTERNAL, &indmeshdat), "mesh induit", proccomm); #ifdef PAMPA_INFO_ADAPT infoPrint ("dmeshInducePart, vertex number after: %d", indmeshdat.vertlocnbr); #endif /* PAMPA_INFO_ADAPT */ //orgmeshptr = &indmeshdat; //// XXX debut temporaire //if (orgmeshptr->proclocnum == 0) { // PAMPA_Mesh m2; // // Initialisation of PaMPA mesh structure // CHECK_FERR(PAMPA_meshInit(&m2), proccomm); // CHECK_FERR(PAMPA_dmeshGather (orgmeshptr, &m2), proccomm); // CHECK_FERR(infoptr->EXT_meshSave(&m2, 1, dataptr, NULL, "apres_gather.mesh"), proccomm); // // Finalisation of PaMPA mesh structure // PAMPA_meshExit(&m2); //} //else { // // Gather PaMPA mesh // CHECK_FERR (PAMPA_dmeshGather (orgmeshptr, NULL), proccomm); //} //// XXX fin temporaire #ifdef PAMPA_INFO_ADAPT infoPrint ("Number of received meshes : %d", meshlocnbr); infoPrint ("Number of local vertices in source distributed mesh : %d", orgmeshptr->vertlocnbr); #endif /* PAMPA_INFO_ADAPT */ memFreeGroup (bnumgsttax + baseval); if (cvrtglbnbr == 1) memFree (parttax + baseval); // 3) Remaillage des zones // a) conversion pampa -> mmg ok // b) remaillage ok // c) conversion mmg -> pampa quasi ok if ((omshloctab = (Mesh *) memAlloc (meshlocnbr * sizeof (Mesh))) == NULL) { errorPrint ("out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); infoptr->cuntval = 0; //for (meshlocnum = 0; meshlocnum < meshlocnbr ; meshlocnum ++) { // XXX cette partie pourrait être parallélisé avec pthread for (meshlocnum = meshlocnbr - 1; meshlocnum >= 0 ; meshlocnum --) { // XXX cette partie pourrait être parallélisé avec pthread Gnum * restrict inbrtab; Gnum * restrict onbrtab; #ifdef PAMPA_DEBUG_ADAPT2 //CHECK_FERR (meshCheck (imshloctab + meshlocnum), proccomm); #endif /* PAMPA_DEBUG_ADAPT2 */ meshInit (omshloctab + meshlocnum); #ifdef PAMPA_DEBUG_ADAPT_SAVE3 { Dmesh tmshdat; char s[100]; int rank; int size; File f; f.mode = "w"; f.pntr = stdout; sprintf (s, "zone-%d-%d-%d-%d.bz2", infoptr->eitenum, infoptr->iitenum, orgmeshptr->proclocnum, meshlocnum); f.name = s; dmeshInit (&tmshdat, comm1p); MPI_Comm_rank (comm1p, &rank); MPI_Comm_size (comm1p, &size); infoPrint ("rank: %d, size: %d", rank, size); CHECK_FERR (dmeshScatter (&tmshdat, imshloctab + meshlocnum, infoptr->typenbr, infoptr->entttab, infoptr->tagtab, infoptr->typetab, (Gnum) 0), proccomm); CHECK_FERR (fileBlockOpen (&f, 1), proccomm); CHECK_FERR (dmeshSave (&tmshdat, f.pntr), proccomm); CHECK_FERR (fileBlockClose (&f, 1), proccomm); infoPrint ("file %s saved", s); dmeshExit (&tmshdat); } #endif /* PAMPA_DEBUG_ADAPT_SAVE3 */ infoPrint ("remeshing zone " GNUMSTRING, meshlocnum); CHECK_FERR_STR (infoptr->EXT_meshAdapt ((PAMPA_Mesh *) (imshloctab + meshlocnum), (PAMPA_Mesh *) (omshloctab + meshlocnum), infoptr, 0), "remesh zone", proccomm); // FIXME remplacer le 0 par un drapeau omshloctab[meshlocnum].idnum = imshloctab[meshlocnum].idnum; omshloctab[meshlocnum].procglbnbr = imshloctab[meshlocnum].procglbnbr; CHECK_FERR (meshValueData (imshloctab + meshlocnum, PAMPA_ENTT_VIRT_PROC, PAMPA_TAG_PART, NULL, NULL, (void **) &inbrtab), orgmeshptr->proccomm); if (meshValueData (omshloctab + meshlocnum, PAMPA_ENTT_VIRT_PROC, PAMPA_TAG_PART, NULL, NULL, (void **) &onbrtab) == 2) { CHECK_FERR (meshValueLink(omshloctab + meshlocnum, (void **) &onbrtab, PAMPA_VALUE_NONE, NULL, NULL, GNUM_MPI, PAMPA_ENTT_VIRT_PROC, PAMPA_TAG_PART), orgmeshptr->proccomm); memCpy (onbrtab, inbrtab, orgmeshptr->procglbnbr * sizeof (Gnum)); } // FIXME : Et s'il y a d'autres valeurs associées, ne devrait-on pas // les recopier, sachant que : // 1) elles sont récupérées lors du rebuild // 2) pour les entt réelles, seule la routine PAMPA_remesh peut le faire #ifdef PAMPA_DEBUG_ADAPT2 //CHECK_FERR (meshCheck (omshloctab + meshlocnum), proccomm); #endif /* PAMPA_DEBUG_ADAPT2 */ meshExit (imshloctab + meshlocnum); } memFree (imshloctab); CHECK_VDBG (cheklocval, proccomm); // 4) Reintegration des zones dans maillage distribué en cours CHECK_FERR_STR (dmeshRebuild (orgmeshptr, meshlocnbr, omshloctab), "reintegration", proccomm); if (infoptr->itdebug != NULL) { //if ((infoptr->itdebug != NULL) && (infoptr->iitenum == 5)) { //CHECK_FERR_STR (dmeshCheck (orgmeshptr), "check", proccomm); //if (orgmeshptr->proclocnum == 0) // infoPrint ("apres le check"); File f; char s[500]; //Dsmesh dmshdat; infoPrint ("je passe la"); sprintf (s, infoptr->itdebug, infoptr->eitenum, infoptr->iitenum); f.name = s; f.mode = "w"; f.pntr = stdout; CHECK_FERR (fileBlockOpen (&f, 1), proccomm); CHECK_FERR (PAMPA_dmeshSaveAll (orgmeshptr, f.pntr), proccomm); fileBlockClose (&f, 1); //MPI_Barrier (orgmeshptr->proccomm); //MPI_Finalize (); //exit (0); //MPI_Barrier (orgmeshptr->proccomm); //dsmeshInit (&dmshdat, orgmeshptr->proccomm); //CHECK_FERR (dmesh2dsmesh(orgmeshptr, &dmshdat), proccomm); //dsmeshExit (&dmshdat); } #ifdef PAMPA_DEBUG_ADAPT2 CHECK_FERR (dmeshCheck (orgmeshptr), proccomm); #endif /* PAMPA_DEBUG_ADAPT2 */ //// XXX debut temporaire //if (orgmeshptr->proclocnum == 0) { // PAMPA_Mesh m2; // // Initialisation of PaMPA mesh structure // CHECK_FERR(PAMPA_meshInit(&m2), proccomm); // CHECK_FERR(PAMPA_dmeshGather (orgmeshptr, &m2), proccomm); // CHECK_FERR(infoptr->EXT_meshSave(&m2, 1, dataptr, NULL, "apres_rebuild.mesh"), proccomm); // // Finalisation of PaMPA mesh structure // PAMPA_meshExit(&m2); //} //else { // // Gather PaMPA mesh // CHECK_FERR (PAMPA_dmeshGather (orgmeshptr, NULL), proccomm); //} //// XXX fin temporaire #ifdef PAMPA_INFO_ADAPT infoPrint ("Number of local vertices in destination distributed mesh : %d", orgmeshptr->vertlocnbr); #endif /* PAMPA_INFO_ADAPT */ #ifdef PAMPA_DEBUG_DMESH_SAVE2 { char s1[50]; sprintf (s1, "dmesh-adapt-apres-%d", wmshcnt ++); CHECK_FERR (dmeshAdaptGatherSave (orgmeshptr, NULL, infoptr->EXT_meshSave, infoptr->dataptr, s1), proccomm); } #endif /* PAMPA_DEBUG_DMESH_SAVE2 */ for (meshlocnum = 0; meshlocnum < meshlocnbr ; meshlocnum ++) meshFree (omshloctab + meshlocnum); memFree (omshloctab); //#ifdef PAMPA_DEBUG_DMESH_SAVE //#ifndef PAMPA_DEBUG_DMESH_SAVE2 // //char s1[50]; //#endif /* PAMPA_DEBUG_DMESH_SAVE2 */ // //sprintf (s1, "dmesh-adapt-%d-%d", orgmeshptr->proclocnum, wmshcnt ++); // //CHECK_FERR (dmeshAdaptGatherSave (orgmeshptr, NULL, dataptr, s1), proccomm); // if (memAllocGroup ((void **) (void *) // &vnumgsttax, (size_t) (orgmeshptr->enttloctax[baseval]->vgstlocnbr * sizeof (Gnum)), // &parttax, (size_t) (orgmeshptr->procglbnbr * sizeof (Gnum)), // NULL) == NULL) { // errorPrint ("Out of memory"); // cheklocval = 1; // } // CHECK_VERR (cheklocval, proccomm); // vnumgsttax -= baseval; // parttax -= baseval; // // for (vertlocnum = baseval, vertlocnnd = orgmeshptr->enttloctax[baseval]->vertlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) // vnumgsttax[vertlocnum] = orgmeshptr->proclocnum; // for (procngbnum = 0; procngbnum < orgmeshptr->procglbnbr; procngbnum ++) // parttax[procngbnum + baseval] = procngbnum; // // imshloctab = NULL; // meshlocnbr = -1; // FIXME l'allocation de imshloctab ne devrait pas être faite ici ? // CHECK_FERR (dmeshGatherInduceMultiple (orgmeshptr, 0, orgmeshptr->procglbnbr, vnumgsttax, parttax, &meshlocnbr, &imshloctab), proccomm); // sprintf (s1, "dmsh-adapt-%d-%d.mesh", orgmeshptr->proclocnum, wmshcnt ++); // CHECK_FERR (infoptr->EXT_meshSave((PAMPA_Mesh *) (imshloctab), 1, dataptr, vnumgsttax + baseval, s1), proccomm); // // meshExit (imshloctab); // memFree (imshloctab); // memFreeGroup (vnumgsttax + baseval); //#endif /* PAMPA_DEBUG_DMESH_SAVE */ CHECK_FERR (dmeshDgraphFree (&dgrfdat), proccomm); { Gnum vertmax; Gnum vertmin; if (commAllreduce (&orgmeshptr->vertlocnbr, &vertmax, 1, GNUM_MPI, MPI_MAX, proccomm) != MPI_SUCCESS) { errorPrint ("communication error"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); if (commAllreduce (&orgmeshptr->vertlocnbr, &vertmin, 1, GNUM_MPI, MPI_MIN, proccomm) != MPI_SUCCESS) { errorPrint ("communication error"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); if (orgmeshptr->proclocnum == 0) infoPrint ("vertmin: " GNUMSTRING ", vertmax: " GNUMSTRING, vertmin, vertmax); if (vertmax > (vertmin * 1.25)) { if (orgmeshptr->proclocnum == 0) infoPrint ("je rentre dedans"); Dmesh orgmeshdat; PAMPA_Strat strtdt2; Gnum * partloctx2; CHECK_FERR (PAMPA_stratInit (&strtdt2), proccomm); CHECK_FERR (dmeshInit (&orgmeshdat, proccomm), proccomm); if ((partloctx2 = (Gnum *) memAlloc (orgmeshptr->vertlocnbr * sizeof (Gnum))) == NULL) { errorPrint ("out of memory"); cheklocval = 1; } CHECK_VERR (cheklocval, proccomm); partloctx2 -= baseval; CHECK_FERR_STR(dmeshMapElem(orgmeshptr, orgmeshptr->procglbnbr, NULL, &strtdt2, partloctx2), "part again", proccomm); CHECK_FERR(dmeshRedist(orgmeshptr, partloctx2, NULL, -1, -1, -1, &orgmeshdat), proccomm); memFree (partloctx2 + baseval); PAMPA_stratExit(&strtdt2); PAMPA_dmeshExit(orgmeshptr); *orgmeshptr = orgmeshdat; #ifdef PAMPA_INFO_ADAPT infoPrint ("Number of local vertices in destination redistributed mesh : %d", orgmeshptr->vertlocnbr); #endif /* PAMPA_INFO_ADAPT */ } } infoptr->iitenum ++; } memFreeGroup (cprcvrttab); //CHECK_FERR_STR (dmeshDgraphBuild (orgmeshptr, DMESH_ENTT_ANY, ~0, NULL, NULL, &dgrfdat, NULL), "graphe final", proccomm); //SCOTCH_dgraphSize (&dgrfdat, NULL, &dvrtlocnbr, NULL, NULL); //if (memAllocGroup ((void **) (void *) // &paroloctax, (size_t) (dvrtlocnbr * sizeof (Gnum)), // &partloctax, (size_t) (dvrtlocnbr * sizeof (Gnum)), // NULL) == NULL) { // errorPrint ("Out of memory"); // cheklocval = 1; //} //CHECK_VERR (cheklocval, proccomm); //paroloctax -= baseval; //partloctax -= baseval; //for (vertlocnum = baseval, vertlocnnd = dvrtlocnbr + baseval; vertlocnum < vertlocnnd; vertlocnum ++) // paroloctax[vertlocnum] = orgmeshptr->proclocnum; //// 5) Repartitionnement //vmloloctab = NULL; //emraval = 1; //CHECK_FERR (SCOTCH_stratDgraphMap (&strtdat, "q{strat=m{asc=b{width=3,bnd=(d{pass=40,dif=1,rem=0}|)f{move=80,pass=-1,bal=0.01},org=f{move=80,pass=-1,bal=0.01}},low=r{job=t,bal=0.01,map=t,poli=S,sep=(m{asc=b{bnd=f{move=120,pass=-1,bal=0.01,type=b},org=f{move=120,pass=-1,bal=0.01,type=b},width=3},low=h{pass=10}f{move=120,pass=-1,bal=0.01,type=b},vert=120,rat=0.8}|m{asc=b{bnd=f{move=120,pass=-1,bal=0.01,type=b},org=f{move=120,pass=-1,bal=0.01,type=b},width=3},low=h{pass=10}f{move=120,pass=-1,bal=0.01,type=b},vert=120,rat=0.8})},vert=10000,rat=0.8,type=0}}"), proccomm); //CHECK_FERR (SCOTCH_dgraphRepart (&dgrfdat, orgmeshptr->procglbnbr, paroloctax + baseval, emraval, vmloloctab, &strtdat, partloctax + baseval), proccomm); dmeshDgraphExit (&dgrfdat); SCOTCH_dgraphExit (&bgrfdat); SCOTCH_stratExit (&strtdat); // 5 bis) Redistribution ok //CHECK_FERR (dmeshRedist (orgmeshptr, partloctax, NULL, -1, -1, -1, dstmeshptr), proccomm); // XXX peut-être utiliser un dsmeshRedist qui serait moins gourmand en mémoire à mon avis //memFreeGroup (paroloctax + baseval); #ifdef PAMPA_DEBUG_ADAPT_SAVE3 MPI_Comm_free (&comm1p); #endif /* PAMPA_DEBUG_ADAPT_SAVE3 */ CHECK_VDBG (cheklocval, proccomm); return (0); }