1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #if defined(PETSC_HAVE_PRAGMATIC) 3 #include <pragmatic/cpragmatic.h> 4 #endif 5 #if defined(PETSC_HAVE_MMG) 6 #include <mmg/libmmg.h> 7 #endif 8 #if defined(PETSC_HAVE_PARMMG) 9 #include <parmmg/libparmmg.h> 10 #endif 11 12 static PetscErrorCode DMPlexLabelToVolumeConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscReal refRatio, PetscReal maxVolumes[]) 13 { 14 PetscInt dim, c; 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 19 refRatio = refRatio == PETSC_DEFAULT ? (PetscReal) ((PetscInt) 1 << dim) : refRatio; 20 for (c = cStart; c < cEnd; c++) { 21 PetscReal vol; 22 PetscInt closureSize = 0, cl; 23 PetscInt *closure = NULL; 24 PetscBool anyRefine = PETSC_FALSE; 25 PetscBool anyCoarsen = PETSC_FALSE; 26 PetscBool anyKeep = PETSC_FALSE; 27 28 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 29 maxVolumes[c - cStart] = vol; 30 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 31 for (cl = 0; cl < closureSize*2; cl += 2) { 32 const PetscInt point = closure[cl]; 33 PetscInt refFlag; 34 35 ierr = DMLabelGetValue(adaptLabel, point, &refFlag);CHKERRQ(ierr); 36 switch (refFlag) { 37 case DM_ADAPT_REFINE: 38 anyRefine = PETSC_TRUE;break; 39 case DM_ADAPT_COARSEN: 40 anyCoarsen = PETSC_TRUE;break; 41 case DM_ADAPT_KEEP: 42 anyKeep = PETSC_TRUE;break; 43 case DM_ADAPT_DETERMINE: 44 break; 45 default: 46 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "DMPlex does not support refinement flag %D\n", refFlag); 47 } 48 if (anyRefine) break; 49 } 50 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 51 if (anyRefine) { 52 maxVolumes[c - cStart] = vol / refRatio; 53 } else if (anyKeep) { 54 maxVolumes[c - cStart] = vol; 55 } else if (anyCoarsen) { 56 maxVolumes[c - cStart] = vol * refRatio; 57 } 58 } 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode DMPlexLabelToMetricConstraint(DM dm, DMLabel adaptLabel, PetscInt cStart, PetscInt cEnd, PetscInt vStart, PetscInt vEnd, PetscReal refRatio, Vec *metricVec) 63 { 64 DM udm, coordDM; 65 PetscSection coordSection; 66 Vec coordinates, mb, mx; 67 Mat A; 68 PetscScalar *metric, *eqns; 69 const PetscReal coarseRatio = refRatio == PETSC_DEFAULT ? PetscSqr(0.5) : 1/refRatio; 70 PetscInt dim, Nv, Neq, c, v; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 75 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 76 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 77 ierr = DMGetLocalSection(coordDM, &coordSection);CHKERRQ(ierr); 78 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 79 Nv = vEnd - vStart; 80 ierr = VecCreateSeq(PETSC_COMM_SELF, Nv*PetscSqr(dim), metricVec);CHKERRQ(ierr); 81 ierr = VecGetArray(*metricVec, &metric);CHKERRQ(ierr); 82 Neq = (dim*(dim+1))/2; 83 ierr = PetscMalloc1(PetscSqr(Neq), &eqns);CHKERRQ(ierr); 84 ierr = MatCreateSeqDense(PETSC_COMM_SELF, Neq, Neq, eqns, &A);CHKERRQ(ierr); 85 ierr = MatCreateVecs(A, &mx, &mb);CHKERRQ(ierr); 86 ierr = VecSet(mb, 1.0);CHKERRQ(ierr); 87 for (c = cStart; c < cEnd; ++c) { 88 const PetscScalar *sol; 89 PetscScalar *cellCoords = NULL; 90 PetscReal e[3], vol; 91 const PetscInt *cone; 92 PetscInt coneSize, cl, i, j, d, r; 93 94 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 95 /* Only works for simplices */ 96 for (i = 0, r = 0; i < dim+1; ++i) { 97 for (j = 0; j < i; ++j, ++r) { 98 for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i*dim+d] - cellCoords[j*dim+d]); 99 /* FORTRAN ORDERING */ 100 switch (dim) { 101 case 2: 102 eqns[0*Neq+r] = PetscSqr(e[0]); 103 eqns[1*Neq+r] = 2.0*e[0]*e[1]; 104 eqns[2*Neq+r] = PetscSqr(e[1]); 105 break; 106 case 3: 107 eqns[0*Neq+r] = PetscSqr(e[0]); 108 eqns[1*Neq+r] = 2.0*e[0]*e[1]; 109 eqns[2*Neq+r] = 2.0*e[0]*e[2]; 110 eqns[3*Neq+r] = PetscSqr(e[1]); 111 eqns[4*Neq+r] = 2.0*e[1]*e[2]; 112 eqns[5*Neq+r] = PetscSqr(e[2]); 113 break; 114 } 115 } 116 } 117 ierr = MatSetUnfactored(A);CHKERRQ(ierr); 118 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 119 ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr); 120 ierr = MatSolve(A, mb, mx);CHKERRQ(ierr); 121 ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr); 122 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 123 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 124 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 125 for (cl = 0; cl < coneSize; ++cl) { 126 const PetscInt v = cone[cl] - vStart; 127 128 if (dim == 2) { 129 metric[v*4+0] += vol*coarseRatio*sol[0]; 130 metric[v*4+1] += vol*coarseRatio*sol[1]; 131 metric[v*4+2] += vol*coarseRatio*sol[1]; 132 metric[v*4+3] += vol*coarseRatio*sol[2]; 133 } else { 134 metric[v*9+0] += vol*coarseRatio*sol[0]; 135 metric[v*9+1] += vol*coarseRatio*sol[1]; 136 metric[v*9+3] += vol*coarseRatio*sol[1]; 137 metric[v*9+2] += vol*coarseRatio*sol[2]; 138 metric[v*9+6] += vol*coarseRatio*sol[2]; 139 metric[v*9+4] += vol*coarseRatio*sol[3]; 140 metric[v*9+5] += vol*coarseRatio*sol[4]; 141 metric[v*9+7] += vol*coarseRatio*sol[4]; 142 metric[v*9+8] += vol*coarseRatio*sol[5]; 143 } 144 } 145 ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr); 146 } 147 for (v = 0; v < Nv; ++v) { 148 const PetscInt *support; 149 PetscInt supportSize, s; 150 PetscReal vol, totVol = 0.0; 151 152 ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr); 153 ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr); 154 for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;} 155 for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol; 156 } 157 ierr = PetscFree(eqns);CHKERRQ(ierr); 158 ierr = VecRestoreArray(*metricVec, &metric);CHKERRQ(ierr); 159 ierr = VecDestroy(&mx);CHKERRQ(ierr); 160 ierr = VecDestroy(&mb);CHKERRQ(ierr); 161 ierr = MatDestroy(&A);CHKERRQ(ierr); 162 ierr = DMDestroy(&udm);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 /* 167 Contains the list of registered DMPlexGenerators routines 168 */ 169 extern PlexGeneratorFunctionList DMPlexGenerateList; 170 171 PetscErrorCode DMPlexRefine_Internal(DM dm, DMLabel adaptLabel, DM *dmRefined) 172 { 173 PlexGeneratorFunctionList fl; 174 PetscErrorCode (*refine)(DM,PetscReal*,DM*); 175 PetscErrorCode (*adapt)(DM,DMLabel,DM*); 176 PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); 177 char genname[PETSC_MAX_PATH_LEN], *name = NULL; 178 PetscReal refinementLimit; 179 PetscReal *maxVolumes; 180 PetscInt dim, cStart, cEnd, c; 181 PetscBool flg, flg2, localized; 182 PetscErrorCode ierr; 183 184 PetscFunctionBegin; 185 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 186 ierr = DMPlexGetRefinementLimit(dm, &refinementLimit);CHKERRQ(ierr); 187 ierr = DMPlexGetRefinementFunction(dm, &refinementFunc);CHKERRQ(ierr); 188 if (refinementLimit == 0.0 && !refinementFunc && !adaptLabel) PetscFunctionReturn(0); 189 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 190 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 191 ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_adaptor", genname, sizeof(genname), &flg);CHKERRQ(ierr); 192 if (flg) name = genname; 193 else { 194 ierr = PetscOptionsGetString(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_generator", genname, sizeof(genname), &flg2);CHKERRQ(ierr); 195 if (flg2) name = genname; 196 } 197 198 fl = DMPlexGenerateList; 199 if (name) { 200 while (fl) { 201 ierr = PetscStrcmp(fl->name,name,&flg);CHKERRQ(ierr); 202 if (flg) { 203 refine = fl->refine; 204 adapt = fl->adaptlabel; 205 goto gotit; 206 } 207 fl = fl->next; 208 } 209 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Grid refiner %s not registered",name); 210 } else { 211 while (fl) { 212 if (fl->dim < 0 || dim-1 == fl->dim) { 213 refine = fl->refine; 214 adapt = fl->adaptlabel; 215 goto gotit; 216 } 217 fl = fl->next; 218 } 219 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No grid refiner of dimension %D registered",dim); 220 } 221 222 gotit: 223 switch (dim) { 224 case 1: 225 case 2: 226 case 3: 227 if (adapt) { 228 ierr = (*adapt)(dm, adaptLabel, dmRefined);CHKERRQ(ierr); 229 } else { 230 ierr = PetscMalloc1(cEnd - cStart, &maxVolumes);CHKERRQ(ierr); 231 if (adaptLabel) { 232 ierr = DMPlexLabelToVolumeConstraint(dm, adaptLabel, cStart, cEnd, PETSC_DEFAULT, maxVolumes);CHKERRQ(ierr); 233 } else if (refinementFunc) { 234 for (c = cStart; c < cEnd; ++c) { 235 PetscReal vol, centroid[3]; 236 237 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, centroid, NULL);CHKERRQ(ierr); 238 ierr = (*refinementFunc)(centroid, &maxVolumes[c-cStart]);CHKERRQ(ierr); 239 } 240 } else { 241 for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = refinementLimit; 242 } 243 ierr = (*refine)(dm, maxVolumes, dmRefined);CHKERRQ(ierr); 244 ierr = PetscFree(maxVolumes);CHKERRQ(ierr); 245 } 246 break; 247 default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh refinement in dimension %D is not supported.", dim); 248 } 249 ((DM_Plex *) (*dmRefined)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation; 250 ierr = DMCopyDisc(dm, *dmRefined);CHKERRQ(ierr); 251 if (localized) {ierr = DMLocalizeCoordinates(*dmRefined);CHKERRQ(ierr);} 252 PetscFunctionReturn(0); 253 } 254 255 PetscErrorCode DMPlexCoarsen_Internal(DM dm, DMLabel adaptLabel, DM *dmCoarsened) 256 { 257 Vec metricVec; 258 PetscInt cStart, cEnd, vStart, vEnd; 259 DMLabel bdLabel = NULL; 260 char bdLabelName[PETSC_MAX_PATH_LEN]; 261 PetscBool localized, flg; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 266 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 267 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 268 ierr = DMPlexLabelToMetricConstraint(dm, adaptLabel, cStart, cEnd, vStart, vEnd, PETSC_DEFAULT, &metricVec);CHKERRQ(ierr); 269 ierr = PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, sizeof(bdLabelName), &flg);CHKERRQ(ierr); 270 if (flg) {ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr);} 271 ierr = DMAdaptMetric_Plex(dm, metricVec, bdLabel, dmCoarsened);CHKERRQ(ierr); 272 ierr = VecDestroy(&metricVec);CHKERRQ(ierr); 273 ((DM_Plex *) (*dmCoarsened)->data)->useHashLocation = ((DM_Plex *) dm->data)->useHashLocation; 274 ierr = DMCopyDisc(dm, *dmCoarsened);CHKERRQ(ierr); 275 if (localized) {ierr = DMLocalizeCoordinates(*dmCoarsened);CHKERRQ(ierr);} 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode DMAdaptLabel_Plex(DM dm, DMLabel adaptLabel, DM *dmAdapted) 280 { 281 IS flagIS; 282 const PetscInt *flags; 283 PetscInt defFlag, minFlag, maxFlag, numFlags, f; 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 ierr = DMLabelGetDefaultValue(adaptLabel, &defFlag);CHKERRQ(ierr); 288 minFlag = defFlag; 289 maxFlag = defFlag; 290 ierr = DMLabelGetValueIS(adaptLabel, &flagIS);CHKERRQ(ierr); 291 ierr = ISGetLocalSize(flagIS, &numFlags);CHKERRQ(ierr); 292 ierr = ISGetIndices(flagIS, &flags);CHKERRQ(ierr); 293 for (f = 0; f < numFlags; ++f) { 294 const PetscInt flag = flags[f]; 295 296 minFlag = PetscMin(minFlag, flag); 297 maxFlag = PetscMax(maxFlag, flag); 298 } 299 ierr = ISRestoreIndices(flagIS, &flags);CHKERRQ(ierr); 300 ierr = ISDestroy(&flagIS);CHKERRQ(ierr); 301 { 302 PetscInt minMaxFlag[2], minMaxFlagGlobal[2]; 303 304 minMaxFlag[0] = minFlag; 305 minMaxFlag[1] = -maxFlag; 306 ierr = MPI_Allreduce(minMaxFlag, minMaxFlagGlobal, 2, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 307 minFlag = minMaxFlagGlobal[0]; 308 maxFlag = -minMaxFlagGlobal[1]; 309 } 310 if (minFlag == maxFlag) { 311 switch (minFlag) { 312 case DM_ADAPT_DETERMINE: 313 *dmAdapted = NULL;break; 314 case DM_ADAPT_REFINE: 315 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr); 316 ierr = DMRefine(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break; 317 case DM_ADAPT_COARSEN: 318 ierr = DMCoarsen(dm, MPI_COMM_NULL, dmAdapted);CHKERRQ(ierr);break; 319 default: 320 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP,"DMPlex does not support refinement flag %D\n", minFlag); 321 } 322 } else { 323 ierr = DMPlexSetRefinementUniform(dm, PETSC_FALSE);CHKERRQ(ierr); 324 ierr = DMPlexRefine_Internal(dm, adaptLabel, dmAdapted);CHKERRQ(ierr); 325 } 326 PetscFunctionReturn(0); 327 } 328 329 330 PetscErrorCode DMAdaptMetricPragmatic_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew) 331 { 332 #if defined(PETSC_HAVE_PRAGMATIC) 333 MPI_Comm comm; 334 const char *bdName = "_boundary_"; 335 #if 0 336 DM odm = dm; 337 #endif 338 DM udm, cdm; 339 DMLabel bdLabelFull; 340 const char *bdLabelName; 341 IS bdIS, globalVertexNum; 342 PetscSection coordSection; 343 Vec coordinates; 344 const PetscScalar *coords, *met; 345 const PetscInt *bdFacesFull, *gV; 346 PetscInt *bdFaces, *bdFaceIds, *l2gv; 347 PetscReal *x, *y, *z, *metric; 348 PetscInt *cells; 349 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v; 350 PetscInt off, maxConeSize, numBdFaces, f, bdSize; 351 PetscBool flg; 352 DMLabel bdLabelNew; 353 PetscReal *coordsNew; 354 PetscInt *bdTags; 355 PetscReal *xNew[3] = {NULL, NULL, NULL}; 356 PetscInt *cellsNew; 357 PetscInt d, numCellsNew, numVerticesNew; 358 PetscInt numCornersNew, fStart, fEnd; 359 PetscMPIInt numProcs; 360 PetscErrorCode ierr; 361 362 PetscFunctionBegin; 363 /* Check for FEM adjacency flags */ 364 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 365 ierr = MPI_Comm_size(comm, &numProcs);CHKERRMPI(ierr); 366 if (bdLabel) { 367 ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr); 368 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 369 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 370 } 371 /* Add overlap for Pragmatic */ 372 #if 0 373 /* Check for overlap by looking for cell in the SF */ 374 if (!overlapped) { 375 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr); 376 if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);} 377 } 378 #endif 379 /* Get mesh information */ 380 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 381 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 382 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 383 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 384 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 385 numCells = cEnd - cStart; 386 if (numCells == 0) { 387 PetscMPIInt rank; 388 389 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 390 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank); 391 } 392 numVertices = vEnd - vStart; 393 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 394 for (c = 0, coff = 0; c < numCells; ++c) { 395 const PetscInt *cone; 396 PetscInt coneSize, cl; 397 398 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 399 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 400 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 401 } 402 ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr); 403 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 404 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 405 for (v = 0, numLocVertices = 0; v < numVertices; ++v) { 406 if (gV[v] >= 0) ++numLocVertices; 407 l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v]; 408 } 409 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 410 ierr = DMDestroy(&udm);CHKERRQ(ierr); 411 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 412 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 413 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 414 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 415 for (v = vStart; v < vEnd; ++v) { 416 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 417 x[v-vStart] = PetscRealPart(coords[off+0]); 418 if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]); 419 if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]); 420 } 421 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 422 /* Get boundary mesh */ 423 ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr); 424 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr); 425 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 426 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 427 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 428 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 429 PetscInt *closure = NULL; 430 PetscInt closureSize, cl; 431 432 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 433 for (cl = 0; cl < closureSize*2; cl += 2) { 434 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 435 } 436 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 437 } 438 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 439 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 440 PetscInt *closure = NULL; 441 PetscInt closureSize, cl; 442 443 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 444 for (cl = 0; cl < closureSize*2; cl += 2) { 445 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 446 } 447 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 448 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 449 else {bdFaceIds[f] = 1;} 450 } 451 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 452 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 453 /* Get metric */ 454 ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr); 455 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 456 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]); 457 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 458 #if 0 459 /* Destroy overlap mesh */ 460 ierr = DMDestroy(&dm);CHKERRQ(ierr); 461 #endif 462 /* Create new mesh */ 463 switch (dim) { 464 case 2: 465 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break; 466 case 3: 467 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break; 468 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim); 469 } 470 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 471 pragmatic_set_metric(metric); 472 pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0); 473 ierr = PetscFree(l2gv);CHKERRQ(ierr); 474 /* Read out mesh */ 475 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 476 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr); 477 switch (dim) { 478 case 2: 479 numCornersNew = 3; 480 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr); 481 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]); 482 break; 483 case 3: 484 numCornersNew = 4; 485 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr); 486 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]); 487 break; 488 default: 489 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim); 490 } 491 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];} 492 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr); 493 pragmatic_get_elements(cellsNew); 494 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr); 495 /* Read out boundary label */ 496 pragmatic_get_boundaryTags(&bdTags); 497 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 498 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 499 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 500 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 501 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 502 for (c = cStart; c < cEnd; ++c) { 503 /* Only for simplicial meshes */ 504 coff = (c-cStart)*(dim+1); 505 /* d is the local cell number of the vertex opposite to the face we are marking */ 506 for (d = 0; d < dim+1; ++d) { 507 if (bdTags[coff+d]) { 508 const PetscInt perm[4][4] = {{-1, -1, -1, -1}, {-1, -1, -1, -1}, {1, 2, 0, -1}, {3, 2, 1, 0}}; /* perm[d] = face opposite */ 509 const PetscInt *cone; 510 511 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 512 ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr); 513 ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr); 514 } 515 } 516 } 517 /* Cleanup */ 518 switch (dim) { 519 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break; 520 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break; 521 } 522 ierr = PetscFree(cellsNew);CHKERRQ(ierr); 523 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr); 524 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 525 ierr = PetscFree(coordsNew);CHKERRQ(ierr); 526 pragmatic_finalize(); 527 PetscFunctionReturn(0); 528 #else 529 PetscFunctionBegin; 530 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic."); 531 #endif 532 } 533 534 int my_increasing_comparison_function(const void *left, const void *right, void *ctx) { 535 PetscInt l = *(PetscInt*) left, r = *(PetscInt *) right; 536 return (l < r) ? -1 : (l > r); 537 } 538 539 PetscErrorCode DMAdaptMetricMMG_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew) 540 { 541 #if defined(PETSC_HAVE_PRAGMATIC) 542 MPI_Comm comm; 543 const char *bdName = "_boundary_"; 544 #if 0 545 DM odm = dm; 546 #endif 547 DM udm, cdm; 548 DMLabel bdLabelFull; 549 const char *bdLabelName; 550 IS bdIS; // IS : index set: ensemble d'indices ordonnés BDIS pour les conditions aux limites , global vertex num pour le // 551 PetscSection coordSection; 552 Vec coordinates; 553 const PetscScalar *coords, *met; // double 554 const PetscInt *bdFacesFull; 555 PetscInt *bdFaces, *bdFaceIds; 556 PetscReal *vertices, *metric, *verticesNew; 557 PetscInt *cells; 558 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, v; 559 PetscInt off, maxConeSize, numBdFaces, f, bdSize, fStart, fEnd,j,k; 560 PetscBool flg; 561 DMLabel bdLabelNew; 562 PetscInt *cellsNew; 563 PetscInt numCellsNew, numVerticesNew; 564 PetscInt numCornersNew; 565 PetscMPIInt numProcs,me; 566 PetscErrorCode ierr; 567 // nos variables 568 PetscInt * tab_cl_verticies, * tab_cl_triangles, *tab_cl_verticies_new, *tab_cl_cells_new; 569 PetscInt * tab_areCorners, * tab_areRequiredCells, *tab_areRequiredVerticies; 570 PetscInt * faces, * tab_cl_faces, * tab_areRidges, * tab_areRequiredFaces; 571 PetscInt numFacesNew; 572 MMG5_pMesh mmg_mesh = NULL; 573 MMG5_pSol mmg_metric = NULL; 574 PetscInt i; 575 // Pour le parallèle 576 PMMG_pParMesh parmesh = NULL; 577 PetscSF starforest; 578 const PetscInt *gV; 579 PetscInt numleaves=0, numroots=0, num_communicators=0,ct=0,ctt=0,p; 580 PetscInt *flags_proc,*communicators_local,*communicators_global,*communicators_local_new; 581 PetscInt *num_per_proc, *offset; 582 IS globalVertexNum; 583 PetscInt *gv_new, *ranks_own,numVerticesNewNew=0; 584 PetscReal *VerticesNewNew; 585 586 PetscFunctionBegin; 587 588 // 0. Début du programme 589 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 590 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr); 591 ierr = MPI_Comm_rank(comm, &me);CHKERRQ(ierr); 592 if (bdLabel) { 593 ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr); 594 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 595 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 596 } 597 598 // 0. Dans le cas parallèlle il faut récupérer les sommets aux interfaces 599 // et le starforest 600 if(numProcs>1){ 601 ierr = DMPlexDistribute(dm, 0, NULL, &udm);dm=udm; 602 } 603 604 // 1. Chercher les informations dans le maillage 605 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 606 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 607 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); /*Récupération des cellulles*/ 608 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); /*Récupération des sommets*/ 609 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); /*On regarde uniquement les cellulles et les sommets*/ 610 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); /*regarde la taille maximum du cône ie combien de sommets au max a une cellule donc si ce sont des triangles, des quadrilatères*/ 611 numCells = cEnd - cStart; /*indice du début moins indice de fin des tranches*/ 612 numVertices = vEnd - vStart; 613 614 printf("DEBUG %d : nombre de tétra %d \n",me, numCells); 615 616 // 2. Récupération des cellules 617 ierr = PetscCalloc1(numCells*maxConeSize, &cells);CHKERRQ(ierr); 618 for (c = 0, coff = 0; c < numCells; ++c) { /*boucle sur les cellules*/ 619 const PetscInt *cone; 620 PetscInt coneSize, cl; 621 622 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); /*récupération de la taille du cone*/ 623 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); /*récupération du cône*/ 624 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart+1; /*translation du tableau*/ 625 } 626 627 // 3. Récupération de tous les sommets 628 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 629 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 630 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 631 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 632 if (dim==2) {ierr=PetscCalloc2(numVertices*3, &metric,2*numVertices, &vertices);CHKERRQ(ierr);} 633 if (dim==3) {ierr=PetscCalloc2(numVertices*6, &metric,3*numVertices, &vertices);CHKERRQ(ierr);} 634 635 for (v = 0; v < vEnd-vStart; ++v) { 636 ierr = PetscSectionGetOffset(coordSection, v+vStart, &off);CHKERRQ(ierr); 637 if (dim==2) { 638 vertices[2*v]=PetscRealPart(coords[off+0]); 639 vertices[2*v+1]=PetscRealPart(coords[off+1]); 640 } 641 else if (dim==3) { 642 vertices[3*v]=PetscRealPart(coords[off+0]); 643 vertices[3*v+1]=PetscRealPart(coords[off+1]); 644 vertices[3*v+2]=PetscRealPart(coords[off+2]); 645 } 646 } 647 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 648 649 // 3.5 récupération des interfaces 650 if (numProcs>1) { 651 652 PetscInt niranks,nranks; 653 const PetscMPIInt *iranks,*ranks; 654 const PetscInt *ioffset, *irootloc,*roffset,*rmine,*rremote; 655 PetscSection coordSection; 656 DM cdm; 657 PetscInt off; 658 Vec coordinates; 659 const PetscScalar *coords; 660 PetscScalar x, y, z; 661 662 ierr = DMGetPointSF(dm,&starforest);CHKERRQ(ierr); 663 ierr = PetscSFSetUp(starforest);CHKERRQ(ierr); 664 ierr = PetscSFGetLeafRanks(starforest, &niranks, &iranks, &ioffset, &irootloc); CHKERRQ(ierr); 665 ierr = PetscSFGetRootRanks(starforest, &nranks, &ranks, &roffset, &rmine, &rremote); CHKERRQ(ierr); 666 667 ierr = PetscCalloc1(numVertices*numProcs,&flags_proc);CHKERRQ(ierr); 668 ierr = PetscCalloc1(numProcs,&num_per_proc);CHKERRQ(ierr); 669 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 670 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 671 672 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 673 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 674 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 675 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 676 /// FIN TESTS 677 678 // Recherche des feuilles 679 for(p=0;p<nranks;p++) 680 for(i=roffset[p];i<roffset[p+1];i++) 681 if (rmine[i] >= vStart && rmine[i] < vEnd && flags_proc[(rmine[i]-vStart)*numProcs+ranks[p]]==0) { 682 numleaves++; 683 flags_proc[(rmine[i]-vStart)*numProcs+ranks[p]]=1; 684 num_per_proc[ranks[p]]++; 685 } 686 687 // recherche des racines 688 for(p=0;p<niranks;p++) 689 for(i=ioffset[p];i<ioffset[p+1];i++) 690 if (irootloc[i] >= vStart && irootloc[i] < vEnd && flags_proc[(irootloc[i]-vStart)*numProcs+iranks[p]]==0) { 691 numroots++; 692 flags_proc[(irootloc[i]-vStart)*numProcs+iranks[p]]=1; 693 num_per_proc[iranks[p]]++; 694 } 695 696 printf("numleaves + numroots %d sur %d total : %d\n",numleaves+numroots,me,numVertices); 697 698 // nombre de comm 699 for (p=0;p<numProcs;p++) if (p!=me && num_per_proc[p] !=0) num_communicators++; 700 701 // construction des tableaux non triés 702 ierr = PetscCalloc1(num_communicators+1,&offset);CHKERRQ(ierr); offset[0]=0; 703 ierr = PetscCalloc3(numroots+numleaves,&communicators_local,numroots+numleaves,&communicators_global,numroots+numleaves,&communicators_local_new); CHKERRQ(ierr); 704 for(p=0;p<numProcs;p++) if (p!=me && num_per_proc[p] !=0) { 705 for(i=0;i<numVertices;i++) if (flags_proc[i*numProcs+p]==1) { 706 communicators_local[ct]=i+1; 707 communicators_global[ct]=gV[i] < 0 ? -(gV[i]+1) : gV[i]; 708 communicators_global[ct]++; 709 ct++; 710 } 711 offset[++ctt]=ct; 712 } 713 // tri du tableau global 714 for(p=0;p<num_communicators;p++) { 715 ierr = PetscTimSort(offset[p+1]-offset[p],&communicators_global[offset[p]],sizeof(PetscInt),my_increasing_comparison_function,NULL); CHKERRQ(ierr); 716 } 717 // reconstruction du tableau local 718 for(p=0;p<numroots+numleaves;p++) { 719 for(i=0;i<numroots+numleaves;i++) { 720 PetscInt tempa=communicators_local[i]-1,tempb; 721 tempb=gV[tempa] < 0 ? -(gV[tempa]+1) : gV[tempa]; 722 tempb++; 723 if (tempb==communicators_global[p]) communicators_local_new[p]=communicators_local[i]; 724 } 725 } 726 727 for (p=0;p<numProcs;p++) { 728 if (p==me) for(v=0;v<num_communicators+1;v++) {printf("%d",offset[v]);printf(" fin offst\n");} 729 if (p==me) for(i=0;i<numroots+numleaves;i++) { 730 ierr = PetscSectionGetOffset(coordSection, vStart+communicators_local_new[i]-1, &off);CHKERRQ(ierr); 731 x = PetscRealPart(coords[off+0]); 732 y = PetscRealPart(coords[off+1]); 733 z = PetscRealPart(coords[off+2]); 734 //printf("%d %d %d %d diff %d coo x %1.2f coo y %1.2f coo z %1.2f \n",me,i,communicators_local_new[i],communicators_global[i],gV[communicators_local_new[i]-1],x,y,z); 735 } 736 737 MPI_Barrier(PETSC_COMM_WORLD); 738 } 739 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 740 } 741 ierr = DMDestroy(&udm);CHKERRQ(ierr); 742 743 744 // 4. Récupératin des conditions aux bords 745 ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr); 746 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr); 747 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 748 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 749 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 750 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 751 PetscInt *closure = NULL; 752 PetscInt closureSize, cl; 753 754 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 755 for (cl = 0; cl < closureSize*2; cl += 2) { 756 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 757 } 758 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 759 } 760 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 761 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 762 PetscInt *closure = NULL; 763 PetscInt closureSize, cl; 764 765 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 766 for (cl = 0; cl < closureSize*2; cl += 2) { 767 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart+1; 768 } 769 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 770 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 771 else {bdFaceIds[f] = 1;} 772 } 773 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 774 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 775 776 // 5. Récupération de la metric 777 ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr); 778 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 779 780 if (dim==2) { 781 for (v = 0; v < (vEnd-vStart); ++v) { // *PetscSqr(dim) 782 metric[3*v] = PetscRealPart(met[4*v]); 783 metric[3*v+1] = PetscRealPart(met[4*v+1]); 784 metric[3*v+2] = PetscRealPart(met[4*v+3]); 785 } 786 } 787 else if (dim==3) { 788 for (v = 0; v < (vEnd-vStart); ++v) { // *PetscSqr(dim) 789 metric[6*v] = PetscRealPart(met[9*v]); 790 metric[6*v+1] = PetscRealPart(met[9*v+1]); 791 metric[6*v+2] = PetscRealPart(met[9*v+2]); 792 metric[6*v+3] = PetscRealPart(met[9*v+4]); 793 metric[6*v+4] = PetscRealPart(met[9*v+5]); 794 metric[6*v+5] = PetscRealPart(met[9*v+8]); 795 } 796 } 797 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 798 799 800 /* PARTIE 2: Transformation du maillage avec mmg*/ 801 ierr = PetscCalloc2(numVertices,&tab_cl_verticies,numCells,&tab_cl_triangles);CHKERRQ(ierr); 802 803 if (numProcs==1) { 804 switch(dim) 805 { 806 case 2: 807 ierr = MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 808 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose,10); // quantité d'information à l'écran 10=toutes les informations 809 ierr = MMG2D_Set_meshSize(mmg_mesh,numVertices,numCells,0,numBdFaces); 810 // Passage des informations sur le maillage 811 812 // géométrie 813 ierr = MMG2D_Set_vertices(mmg_mesh,vertices,tab_cl_verticies); 814 ierr = MMG2D_Set_triangles(mmg_mesh,cells,tab_cl_triangles); 815 ierr = MMG2D_Set_edges(mmg_mesh,bdFaces,bdFaceIds); 816 817 // métrique 818 ierr = MMG2D_Set_solSize(mmg_mesh,mmg_metric,MMG5_Vertex,numVertices,MMG5_Tensor); 819 for (i=0;i<numVertices;i++) MMG2D_Set_tensorSol(mmg_metric,metric[3*i],metric[3*i+1],metric[3*i+2],i+1); 820 821 // Remaillage 822 ierr = MMG2D_saveMshMesh(mmg_mesh,NULL,"maillage_avant.msh"); 823 ierr = MMG2D_mmg2dlib(mmg_mesh,mmg_metric); printf("DEBUG remaillage 2D: %d \n", ierr); 824 ierr = MMG2D_saveMshMesh(mmg_mesh,NULL,"maillage_apres.msh"); 825 break; 826 827 case 3: 828 ierr = MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 829 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose,10); 830 ierr = MMG3D_Set_meshSize(mmg_mesh,numVertices,numCells,0,numBdFaces,0,0); 831 832 ierr = MMG3D_Set_vertices(mmg_mesh,vertices,tab_cl_verticies); 833 for (i=0;i<numCells;i++) ierr = MMG3D_Set_tetrahedron(mmg_mesh,cells[4*i+0],cells[4*i+1],cells[4*i+2],cells[4*i+3],0,i+1); 834 ierr = MMG3D_Set_triangles(mmg_mesh,bdFaces,bdFaceIds); 835 836 // Métrique 837 ierr = MMG3D_Set_solSize(mmg_mesh,mmg_metric,MMG5_Vertex,numVertices,MMG5_Tensor); 838 ierr = MMG3D_Set_tensorSols(mmg_metric,metric); 839 840 // Remaillage 841 ierr = MMG3D_saveMshMesh(mmg_mesh,NULL,"maillage_avant_3D.msh"); 842 ierr = MMG3D_mmg3dlib(mmg_mesh,mmg_metric); printf("DEBUG remaillage 3D: %d \n", ierr); 843 ierr = MMG3D_saveMshMesh(mmg_mesh,NULL,"maillage_apres_3D.msh"); 844 break; 845 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 846 } 847 } 848 else 849 { 850 // Initialisation 851 ierr = PMMG_Init_parMesh(PMMG_ARG_start,PMMG_ARG_ppParMesh,&parmesh,PMMG_ARG_pMesh,PMMG_ARG_pMet,PMMG_ARG_dim,3,PMMG_ARG_MPIComm,comm,PMMG_ARG_end); 852 ierr = PMMG_Set_meshSize(parmesh,numVertices,numCells,0,numBdFaces,0,0); 853 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_APImode,PMMG_APIDISTRIB_nodes); 854 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_verbose,10); 855 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_globalNum,1); 856 857 // Maillage et métrique 858 ierr = PMMG_Set_vertices(parmesh,vertices,tab_cl_verticies); 859 for (i=0;i<numCells;i++) ierr = PMMG_Set_tetrahedron(parmesh,cells[4*i+0],cells[4*i+1],cells[4*i+2],cells[4*i+3],0,i+1); 860 ierr = PMMG_Set_triangles(parmesh,bdFaces,bdFaceIds); 861 ierr = PMMG_Set_metSize(parmesh,MMG5_Vertex,numVertices,MMG5_Tensor); 862 for (i=0;i<numVertices;i++) PMMG_Set_tensorMet(parmesh,metric[6*i],metric[6*i+1],metric[6*i+2],metric[6*i+3],metric[6*i+4],metric[6*i+5],i+1); 863 864 // Interfaces de communication 865 ierr = PMMG_Set_numberOfNodeCommunicators(parmesh,num_communicators); 866 printf("nombre de communicateurs sur le proc %d : %d \n",me,num_communicators); 867 for(ct=0,p=0;p<numProcs;p++) 868 if (num_per_proc[p] !=0 && p!=me) { 869 printf("DEBUG %d, taille: %d ct: %d p:%d \n",me,offset[ct+1]-offset[ct],ct,p); 870 ierr = PMMG_Set_ithNodeCommunicatorSize(parmesh,ct,p,offset[ct+1]-offset[ct]); printf("DEBUG communicator size: %d \n", ierr); 871 ierr = PMMG_Set_ithNodeCommunicator_nodes(parmesh,ct,&communicators_local_new[offset[ct]],&communicators_global[offset[ct]],1); printf("DEBUG communicator: %d \n", ierr); 872 873 // printf("DEBUG (%d) ct: %d\n", me, ct); 874 // printf("DEBUG(%d) communicators_local_new: ", me); 875 // for (int iv=offset[ct]; iv<offset[ct+1]; ++iv) { 876 // printf(" %d,", communicators_local_new[iv]); 877 // } 878 // printf("\n"); 879 // printf("DEBUG(%d) communicators_global: ", me); 880 // for (int iv=offset[ct]; iv<offset[ct+1]; ++iv) { 881 // printf(" %d,", communicators_global[iv]); 882 // } 883 // printf("\n"); 884 885 ct++; 886 } 887 888 // Remaillage 889 MPI_Barrier(comm); 890 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_globalNum,1); 891 //ierr = PMMG_saveMesh_distributed(parmesh,"mesh_depart"); 892 ierr = PMMG_parmmglib_distributed(parmesh);printf("DEBUG remaillage //: %d \n", ierr); 893 //ierr = PMMG_saveMesh_distributed(parmesh,"mesh_apres"); 894 } 895 896 /*3. Passer du nouveau maillage mmg à un maillage lisible par petsc*/ 897 if (numProcs==1) { 898 switch(dim) 899 { 900 case(2): 901 numCornersNew = 3; 902 ierr = MMG2D_Get_meshSize(mmg_mesh,&numVerticesNew,&numCellsNew,0,&numFacesNew); 903 904 ierr = PetscCalloc4(2*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 905 ierr = PetscCalloc3(3*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 906 ierr = PetscCalloc4(2*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 907 ierr = MMG2D_Get_vertices(mmg_mesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 908 ierr = MMG2D_Get_triangles(mmg_mesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 909 ierr = MMG2D_Get_edges(mmg_mesh,faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces); 910 911 for (i=0;i<3*numCellsNew;i++) cellsNew[i]-=1; 912 for (i=0;i<2*numFacesNew;i++) faces[i]-=1; 913 914 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, dmNew);CHKERRQ(ierr); 915 ierr = MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); printf("DEBUG libération mémoire : %d \n", ierr); 916 break; 917 case(3): 918 numCornersNew = 4; 919 ierr = MMG3D_Get_meshSize(mmg_mesh,&numVerticesNew,&numCellsNew,0,&numFacesNew,0,0); 920 ierr = PetscCalloc4(3*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 921 ierr = PetscCalloc3(4*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 922 ierr = PetscCalloc4(3*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 923 ierr = MMG3D_Get_vertices(mmg_mesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 924 ierr = MMG3D_Get_tetrahedra(mmg_mesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 925 ierr = MMG3D_Get_triangles(mmg_mesh,faces,tab_cl_faces,tab_areRequiredFaces); 926 927 for (i=0;i<4*numCellsNew;i++) cellsNew[i]-=1; 928 for (i=0;i<3*numFacesNew;i++) faces[i]-=1; 929 930 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, dmNew);CHKERRQ(ierr); 931 ierr = MMG3D_Free_all(MMG5_ARG_start,MMG5_ARG_ppMesh,&mmg_mesh,MMG5_ARG_ppMet,&mmg_metric,MMG5_ARG_end); printf("DEBUG libération mémoire : %d \n", ierr); 932 break; 933 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 934 } 935 } 936 else { 937 numCornersNew = 4; 938 ierr = PMMG_Get_meshSize(parmesh,&numVerticesNew,&numCellsNew,0,&numFacesNew,0,0); 939 ierr = PetscCalloc4(3*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 940 ierr = PetscCalloc3(4*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 941 ierr = PetscCalloc4(3*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 942 ierr = PMMG_Get_vertices(parmesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 943 ierr = PMMG_Get_tetrahedra(parmesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 944 ierr = PMMG_Get_triangles(parmesh,faces,tab_cl_faces,tab_areRequiredFaces); 945 946 ierr = PetscCalloc2(numVerticesNew,&gv_new,numVerticesNew,&ranks_own); 947 ierr = PMMG_Get_verticesGloNum(parmesh,gv_new,ranks_own); 948 949 // Décallage et changement de numérotation 950 for (i=0;i<3*numFacesNew;i++) faces[i]-=1; 951 for (i=0;i<4*numCellsNew;i++) cellsNew[i]=gv_new[cellsNew[i]-1]-1; 952 953 // Calcul de la liste des sommets sur chacun des procs 954 for (i=0;i<numVerticesNew;i++) if (ranks_own[i]==me) numVerticesNewNew++; 955 ierr = PetscCalloc1(numVerticesNewNew,&VerticesNewNew); CHKERRQ(ierr); 956 for (ct=0,i=0;i<numVerticesNew;i++) if (ranks_own[i]==me) { 957 VerticesNewNew[ct++]=verticesNew[i]; 958 } 959 960 // tests 961 for (p=0;p<numProcs;p++) { 962 if (p==me) printf("me : %d , numVerticesNew : %d \n",me,numVerticesNew); 963 if (p==me) { for (i=0;i<numVerticesNew;i++) printf("%d ",gv_new[i]); printf("\n");} 964 if (p==me) { for (i=0;i<numVerticesNew;i++) printf("%d ",ranks_own[i]); printf("\n");} 965 MPI_Barrier(comm); 966 } 967 for(i=0;i<4*numCellsNew;i++) if (cellsNew[i]<0) printf("DEBUG -1 %d %d\n",i,cellsNew[i]); 968 969 for (p=0;p<numProcs;p++) { 970 if (p==me) { 971 for (i=0;i<numVerticesNew;i++) { 972 printf("DBEUG(%d) VerticesNewNew[%d]: %1.2f %1.2f %1.2f\n", me, i, 973 VerticesNewNew[3*i], VerticesNewNew[3*i+1], VerticesNewNew[3*i+2]); 974 } 975 } 976 MPI_Barrier(comm); 977 } 978 979 // printf("DEBUG %d nombre de sommets %d, nombre de tétra %d num total: %d \n",me,numVerticesNew,numCellsNew,num_total); 980 // if (me==1) for (i=0;i<numVerticesNew;i++) printf("%d sommets %d %lf %lf %lf\n",me,i,verticesNew[3*i],verticesNew[3*i+1],verticesNew[3*i+2]); 981 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, VerticesNewNew, NULL, dmNew);CHKERRQ(ierr); 982 ierr = PMMG_Free_all(PMMG_ARG_start,PMMG_ARG_ppParMesh,&parmesh,PMMG_ARG_end); 983 984 } 985 986 /*Reconstruction des conditions aux limites*/ 987 //ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 988 //ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 989 //for(i=0;i<numFacesNew;i++) ierr = DMLabelSetValue(bdLabelNew,faces[i],tab_cl_faces[i]);CHKERRQ(ierr); 990 991 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 992 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 993 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 994 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 995 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 996 997 switch(dim) 998 { 999 case(2): 1000 for (i=0;i<numFacesNew;i++){ 1001 const PetscInt *supportA, *supportB; 1002 PetscInt SA,SB,inter=-1; 1003 ierr = DMPlexGetSupportSize(*dmNew,faces[2*i]+vStart,&SA); 1004 ierr = DMPlexGetSupportSize(*dmNew,faces[2*i+1]+vStart,&SB); 1005 ierr = DMPlexGetSupport(*dmNew,faces[2*i]+vStart,&supportA); 1006 ierr = DMPlexGetSupport(*dmNew,faces[2*i+1]+vStart,&supportB); 1007 /*printf("Numero des points : %d %d \n",faces[2*i]+1,faces[2*i+1]+1); 1008 printf("Taille des supports : %d %d \n",SA,SB); 1009 printf("Support de A ");for(d=0;d<SA;d++) printf("%d ",supportA[d]); printf("\n"); 1010 printf("Support de B ");for(d=0;d<SB;d++) printf("%d ",supportB[d]); printf("\n");*/ 1011 1012 // Calcul de l'intersection: 1013 for(j=0;j<SA;j++){ 1014 for(k=0;k<SB;k++){ 1015 if(supportA[j]==supportB[k]) inter=supportA[j]; 1016 } 1017 } 1018 ierr = DMLabelSetValue(bdLabelNew,inter,tab_cl_faces[i]);CHKERRQ(ierr); 1019 } 1020 break; 1021 case(3): 1022 for (i=0;i<numFacesNew;i++){ 1023 const PetscInt *supportA, *supportB, *supportC; 1024 const PetscInt *supportiA, *supportiB, *supportiC; 1025 1026 PetscInt SA,SB,SC,inter=-1; 1027 PetscInt SiA,SiB,SiC; 1028 PetscInt ia, ib, ic; 1029 PetscInt ja, jb, jc; 1030 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i]+vStart,&SA); 1031 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i+1]+vStart,&SB); 1032 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i+2]+vStart,&SC); 1033 ierr = DMPlexGetSupport(*dmNew,faces[3*i]+vStart,&supportA); 1034 ierr = DMPlexGetSupport(*dmNew,faces[3*i+1]+vStart,&supportB); 1035 ierr = DMPlexGetSupport(*dmNew,faces[3*i+2]+vStart,&supportC); 1036 1037 //printf("%d %d %d\n",SA,SB,SC); 1038 inter=-1; 1039 for (ia=0;ia<SA;ia++) { 1040 ierr = DMPlexGetSupportSize(*dmNew,supportA[ia],&SiA); 1041 ierr = DMPlexGetSupport(*dmNew,supportA[ia],&supportiA); 1042 for (ib=0;ib<SB;ib++) { 1043 ierr = DMPlexGetSupportSize(*dmNew,supportB[ib],&SiB); 1044 ierr = DMPlexGetSupport(*dmNew,supportB[ib],&supportiB); 1045 for (ic=0;ic<SC;ic++) { 1046 ierr = DMPlexGetSupportSize(*dmNew,supportC[ic],&SiC); 1047 ierr = DMPlexGetSupport(*dmNew,supportC[ic],&supportiC); 1048 for(ja=0;ja<SiA;ja++) { 1049 for(jb=0;jb<SiB;jb++) { 1050 for(jc=0;jc<SiC;jc++) { 1051 if(supportiA[ja]==supportiB[jb] && supportiA[ja]==supportiC[jc]) inter=supportiA[ja]; 1052 } 1053 } 1054 } 1055 } 1056 } 1057 } 1058 ierr = DMLabelSetValue(bdLabelNew,inter,tab_cl_faces[i]);CHKERRQ(ierr); 1059 } 1060 break; 1061 default:SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 1062 } 1063 1064 /*libération des tableaux en mémoire TO DO*/ 1065 ierr = PetscFree3(cells,metric,vertices);CHKERRQ(ierr); printf("debug 1 \n"); 1066 ierr = PetscFree4(bdFaces,bdFaceIds,tab_cl_verticies,tab_cl_triangles);CHKERRQ(ierr);printf("debug 2 \n"); 1067 ierr = PetscFree4(verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies);CHKERRQ(ierr);printf("debug 3 \n"); 1068 ierr = PetscFree4(faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces);CHKERRQ(ierr);printf("debug 4 \n"); 1069 ierr = PetscFree4(faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces);CHKERRQ(ierr);printf("debug 5 \n"); 1070 1071 if (numProcs>1) ierr = PetscFree4(communicators_local,communicators_global,num_per_proc,flags_proc); CHKERRQ(ierr); 1072 if (numProcs>1) ierr = PetscFree3(VerticesNewNew,ranks_own,gv_new); CHKERRQ(ierr); 1073 1074 PetscFunctionReturn(0); 1075 #else 1076 PetscFunctionBegin; 1077 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic."); 1078 PetscFunctionReturn(0); 1079 #endif 1080 } 1081 1082 PetscErrorCode DMAdaptMetricParMMG_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew) 1083 { 1084 #if defined(PETSC_HAVE_PRAGMATIC) 1085 MPI_Comm comm; 1086 const char *bdName = "_boundary_"; 1087 #if 0 1088 DM odm = dm; 1089 #endif 1090 DM udm, cdm; 1091 DMLabel bdLabelFull; 1092 const char *bdLabelName; 1093 IS bdIS; // IS : index set: ensemble d'indices ordonnés BDIS pour les conditions aux limites , global vertex num pour le // 1094 PetscSection coordSection; 1095 Vec coordinates; 1096 const PetscScalar *coords, *met; // double 1097 const PetscInt *bdFacesFull; 1098 PetscInt *bdFaces, *bdFaceIds; 1099 PetscReal *vertices, *metric, *verticesNew; 1100 PetscInt *cells; 1101 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, v; 1102 PetscInt off, maxConeSize, numBdFaces, f, bdSize, fStart, fEnd,j,k; 1103 PetscBool flg; 1104 DMLabel bdLabelNew; 1105 PetscInt *cellsNew; 1106 PetscInt numCellsNew, numVerticesNew; 1107 PetscInt numCornersNew; 1108 PetscMPIInt numProcs,me; 1109 PetscErrorCode ierr; 1110 // nos variables 1111 PetscInt * tab_cl_verticies, * tab_cl_triangles, *tab_cl_verticies_new, *tab_cl_cells_new; 1112 PetscInt * tab_areCorners, * tab_areRequiredCells, *tab_areRequiredVerticies; 1113 PetscInt * faces, * tab_cl_faces, * tab_areRidges, * tab_areRequiredFaces; 1114 PetscInt numFacesNew; 1115 MMG5_pMesh mmg_mesh = NULL; 1116 MMG5_pSol mmg_metric = NULL; 1117 PetscInt i; 1118 // Pour le parallèle 1119 PMMG_pParMesh parmesh = NULL; 1120 PetscSF starforest; 1121 const PetscInt *gV; 1122 PetscInt numleaves=0, numroots=0, num_communicators=0,ct=0,ctt=0,p; 1123 PetscInt *flags_proc,*communicators_local,*communicators_global,*communicators_local_new; 1124 PetscInt *num_per_proc, *offset; 1125 IS globalVertexNum; 1126 PetscInt *gv_new, *ranks_own,numVerticesNewNew=0; 1127 PetscReal *VerticesNewNew; 1128 1129 PetscFunctionBegin; 1130 1131 // 0. Début du programme 1132 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 1133 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr); 1134 ierr = MPI_Comm_rank(comm, &me);CHKERRQ(ierr); 1135 if (bdLabel) { 1136 ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr); 1137 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 1138 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 1139 } 1140 1141 // 0. Dans le cas parallèlle il faut récupérer les sommets aux interfaces 1142 // et le starforest 1143 if(numProcs>1){ 1144 ierr = DMPlexDistribute(dm, 0, NULL, &udm);dm=udm; 1145 } 1146 1147 // 1. Chercher les informations dans le maillage 1148 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1149 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1150 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); /*Récupération des cellulles*/ 1151 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); /*Récupération des sommets*/ 1152 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); /*On regarde uniquement les cellulles et les sommets*/ 1153 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); /*regarde la taille maximum du cône ie combien de sommets au max a une cellule donc si ce sont des triangles, des quadrilatères*/ 1154 numCells = cEnd - cStart; /*indice du début moins indice de fin des tranches*/ 1155 numVertices = vEnd - vStart; 1156 1157 printf("DEBUG %d : nombre de tétra %d \n",me, numCells); 1158 1159 // 2. Récupération des cellules 1160 ierr = PetscCalloc1(numCells*maxConeSize, &cells);CHKERRQ(ierr); 1161 for (c = 0, coff = 0; c < numCells; ++c) { /*boucle sur les cellules*/ 1162 const PetscInt *cone; 1163 PetscInt coneSize, cl; 1164 1165 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); /*récupération de la taille du cone*/ 1166 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); /*récupération du cône*/ 1167 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart+1; /*translation du tableau*/ 1168 } 1169 1170 // 3. Récupération de tous les sommets 1171 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1172 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 1173 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1174 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 1175 if (dim==2) {ierr=PetscCalloc2(numVertices*3, &metric,2*numVertices, &vertices);CHKERRQ(ierr);} 1176 if (dim==3) {ierr=PetscCalloc2(numVertices*6, &metric,3*numVertices, &vertices);CHKERRQ(ierr);} 1177 1178 for (v = 0; v < vEnd-vStart; ++v) { 1179 ierr = PetscSectionGetOffset(coordSection, v+vStart, &off);CHKERRQ(ierr); 1180 if (dim==2) { 1181 vertices[2*v]=PetscRealPart(coords[off+0]); 1182 vertices[2*v+1]=PetscRealPart(coords[off+1]); 1183 } 1184 else if (dim==3) { 1185 vertices[3*v]=PetscRealPart(coords[off+0]); 1186 vertices[3*v+1]=PetscRealPart(coords[off+1]); 1187 vertices[3*v+2]=PetscRealPart(coords[off+2]); 1188 } 1189 } 1190 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 1191 1192 // 3.5 récupération des interfaces 1193 if (numProcs>1) { 1194 1195 PetscInt niranks,nranks; 1196 const PetscMPIInt *iranks,*ranks; 1197 const PetscInt *ioffset, *irootloc,*roffset,*rmine,*rremote; 1198 PetscSection coordSection; 1199 DM cdm; 1200 PetscInt off; 1201 Vec coordinates; 1202 const PetscScalar *coords; 1203 PetscScalar x, y, z; 1204 1205 ierr = DMGetPointSF(dm,&starforest);CHKERRQ(ierr); 1206 ierr = PetscSFSetUp(starforest);CHKERRQ(ierr); 1207 ierr = PetscSFGetLeafRanks(starforest, &niranks, &iranks, &ioffset, &irootloc); CHKERRQ(ierr); 1208 ierr = PetscSFGetRootRanks(starforest, &nranks, &ranks, &roffset, &rmine, &rremote); CHKERRQ(ierr); 1209 1210 ierr = PetscCalloc1(numVertices*numProcs,&flags_proc);CHKERRQ(ierr); 1211 ierr = PetscCalloc1(numProcs,&num_per_proc);CHKERRQ(ierr); 1212 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 1213 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 1214 1215 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 1216 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 1217 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 1218 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 1219 /// FIN TESTS 1220 1221 // Recherche des feuilles 1222 for(p=0;p<nranks;p++) 1223 for(i=roffset[p];i<roffset[p+1];i++) 1224 if (rmine[i] >= vStart && rmine[i] < vEnd && flags_proc[(rmine[i]-vStart)*numProcs+ranks[p]]==0) { 1225 numleaves++; 1226 flags_proc[(rmine[i]-vStart)*numProcs+ranks[p]]=1; 1227 num_per_proc[ranks[p]]++; 1228 } 1229 1230 // recherche des racines 1231 for(p=0;p<niranks;p++) 1232 for(i=ioffset[p];i<ioffset[p+1];i++) 1233 if (irootloc[i] >= vStart && irootloc[i] < vEnd && flags_proc[(irootloc[i]-vStart)*numProcs+iranks[p]]==0) { 1234 numroots++; 1235 flags_proc[(irootloc[i]-vStart)*numProcs+iranks[p]]=1; 1236 num_per_proc[iranks[p]]++; 1237 } 1238 1239 printf("numleaves + numroots %d sur %d total : %d\n",numleaves+numroots,me,numVertices); 1240 1241 // nombre de comm 1242 for (p=0;p<numProcs;p++) if (p!=me && num_per_proc[p] !=0) num_communicators++; 1243 1244 // construction des tableaux non triés 1245 ierr = PetscCalloc1(num_communicators+1,&offset);CHKERRQ(ierr); offset[0]=0; 1246 ierr = PetscCalloc3(numroots+numleaves,&communicators_local,numroots+numleaves,&communicators_global,numroots+numleaves,&communicators_local_new); CHKERRQ(ierr); 1247 for(p=0;p<numProcs;p++) if (p!=me && num_per_proc[p] !=0) { 1248 for(i=0;i<numVertices;i++) if (flags_proc[i*numProcs+p]==1) { 1249 communicators_local[ct]=i+1; 1250 communicators_global[ct]=gV[i] < 0 ? -(gV[i]+1) : gV[i]; 1251 communicators_global[ct]++; 1252 ct++; 1253 } 1254 offset[++ctt]=ct; 1255 } 1256 // tri du tableau global 1257 for(p=0;p<num_communicators;p++) { 1258 ierr = PetscTimSort(offset[p+1]-offset[p],&communicators_global[offset[p]],sizeof(PetscInt),my_increasing_comparison_function,NULL); CHKERRQ(ierr); 1259 } 1260 // reconstruction du tableau local 1261 for(p=0;p<numroots+numleaves;p++) { 1262 for(i=0;i<numroots+numleaves;i++) { 1263 PetscInt tempa=communicators_local[i]-1,tempb; 1264 tempb=gV[tempa] < 0 ? -(gV[tempa]+1) : gV[tempa]; 1265 tempb++; 1266 if (tempb==communicators_global[p]) communicators_local_new[p]=communicators_local[i]; 1267 } 1268 } 1269 1270 for (p=0;p<numProcs;p++) { 1271 if (p==me) for(v=0;v<num_communicators+1;v++) {printf("%d",offset[v]);printf(" fin offst\n");} 1272 if (p==me) for(i=0;i<numroots+numleaves;i++) { 1273 ierr = PetscSectionGetOffset(coordSection, vStart+communicators_local_new[i]-1, &off);CHKERRQ(ierr); 1274 x = PetscRealPart(coords[off+0]); 1275 y = PetscRealPart(coords[off+1]); 1276 z = PetscRealPart(coords[off+2]); 1277 //printf("%d %d %d %d diff %d coo x %1.2f coo y %1.2f coo z %1.2f \n",me,i,communicators_local_new[i],communicators_global[i],gV[communicators_local_new[i]-1],x,y,z); 1278 } 1279 1280 MPI_Barrier(PETSC_COMM_WORLD); 1281 } 1282 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 1283 } 1284 ierr = DMDestroy(&udm);CHKERRQ(ierr); 1285 1286 1287 // 4. Récupératin des conditions aux bords 1288 ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr); 1289 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr); 1290 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 1291 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 1292 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 1293 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 1294 PetscInt *closure = NULL; 1295 PetscInt closureSize, cl; 1296 1297 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1298 for (cl = 0; cl < closureSize*2; cl += 2) { 1299 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 1300 } 1301 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1302 } 1303 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 1304 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 1305 PetscInt *closure = NULL; 1306 PetscInt closureSize, cl; 1307 1308 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1309 for (cl = 0; cl < closureSize*2; cl += 2) { 1310 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart+1; 1311 } 1312 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 1313 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 1314 else {bdFaceIds[f] = 1;} 1315 } 1316 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 1317 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 1318 1319 // 5. Récupération de la metric 1320 ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr); 1321 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 1322 1323 if (dim==2) { 1324 for (v = 0; v < (vEnd-vStart); ++v) { // *PetscSqr(dim) 1325 metric[3*v] = PetscRealPart(met[4*v]); 1326 metric[3*v+1] = PetscRealPart(met[4*v+1]); 1327 metric[3*v+2] = PetscRealPart(met[4*v+3]); 1328 } 1329 } 1330 else if (dim==3) { 1331 for (v = 0; v < (vEnd-vStart); ++v) { // *PetscSqr(dim) 1332 metric[6*v] = PetscRealPart(met[9*v]); 1333 metric[6*v+1] = PetscRealPart(met[9*v+1]); 1334 metric[6*v+2] = PetscRealPart(met[9*v+2]); 1335 metric[6*v+3] = PetscRealPart(met[9*v+4]); 1336 metric[6*v+4] = PetscRealPart(met[9*v+5]); 1337 metric[6*v+5] = PetscRealPart(met[9*v+8]); 1338 } 1339 } 1340 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 1341 1342 1343 /* PARTIE 2: Transformation du maillage avec mmg*/ 1344 ierr = PetscCalloc2(numVertices,&tab_cl_verticies,numCells,&tab_cl_triangles);CHKERRQ(ierr); 1345 1346 if (numProcs==1) { 1347 switch(dim) 1348 { 1349 case 2: 1350 ierr = MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 1351 ierr = MMG2D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose,10); // quantité d'information à l'écran 10=toutes les informations 1352 ierr = MMG2D_Set_meshSize(mmg_mesh,numVertices,numCells,0,numBdFaces); 1353 // Passage des informations sur le maillage 1354 1355 // géométrie 1356 ierr = MMG2D_Set_vertices(mmg_mesh,vertices,tab_cl_verticies); 1357 ierr = MMG2D_Set_triangles(mmg_mesh,cells,tab_cl_triangles); 1358 ierr = MMG2D_Set_edges(mmg_mesh,bdFaces,bdFaceIds); 1359 1360 // métrique 1361 ierr = MMG2D_Set_solSize(mmg_mesh,mmg_metric,MMG5_Vertex,numVertices,MMG5_Tensor); 1362 for (i=0;i<numVertices;i++) MMG2D_Set_tensorSol(mmg_metric,metric[3*i],metric[3*i+1],metric[3*i+2],i+1); 1363 1364 // Remaillage 1365 ierr = MMG2D_saveMshMesh(mmg_mesh,NULL,"maillage_avant.msh"); 1366 ierr = MMG2D_mmg2dlib(mmg_mesh,mmg_metric); printf("DEBUG remaillage 2D: %d \n", ierr); 1367 ierr = MMG2D_saveMshMesh(mmg_mesh,NULL,"maillage_apres.msh"); 1368 break; 1369 1370 case 3: 1371 ierr = MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); 1372 ierr = MMG3D_Set_iparameter(mmg_mesh, mmg_metric, MMG2D_IPARAM_verbose,10); 1373 ierr = MMG3D_Set_meshSize(mmg_mesh,numVertices,numCells,0,numBdFaces,0,0); 1374 1375 ierr = MMG3D_Set_vertices(mmg_mesh,vertices,tab_cl_verticies); 1376 for (i=0;i<numCells;i++) ierr = MMG3D_Set_tetrahedron(mmg_mesh,cells[4*i+0],cells[4*i+1],cells[4*i+2],cells[4*i+3],0,i+1); 1377 ierr = MMG3D_Set_triangles(mmg_mesh,bdFaces,bdFaceIds); 1378 1379 // Métrique 1380 ierr = MMG3D_Set_solSize(mmg_mesh,mmg_metric,MMG5_Vertex,numVertices,MMG5_Tensor); 1381 ierr = MMG3D_Set_tensorSols(mmg_metric,metric); 1382 1383 // Remaillage 1384 ierr = MMG3D_saveMshMesh(mmg_mesh,NULL,"maillage_avant_3D.msh"); 1385 ierr = MMG3D_mmg3dlib(mmg_mesh,mmg_metric); printf("DEBUG remaillage 3D: %d \n", ierr); 1386 ierr = MMG3D_saveMshMesh(mmg_mesh,NULL,"maillage_apres_3D.msh"); 1387 break; 1388 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 1389 } 1390 } 1391 else 1392 { 1393 // Initialisation 1394 ierr = PMMG_Init_parMesh(PMMG_ARG_start,PMMG_ARG_ppParMesh,&parmesh,PMMG_ARG_pMesh,PMMG_ARG_pMet,PMMG_ARG_dim,3,PMMG_ARG_MPIComm,comm,PMMG_ARG_end); 1395 ierr = PMMG_Set_meshSize(parmesh,numVertices,numCells,0,numBdFaces,0,0); 1396 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_APImode,PMMG_APIDISTRIB_nodes); 1397 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_verbose,10); 1398 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_globalNum,1); 1399 1400 // Maillage et métrique 1401 ierr = PMMG_Set_vertices(parmesh,vertices,tab_cl_verticies); 1402 for (i=0;i<numCells;i++) ierr = PMMG_Set_tetrahedron(parmesh,cells[4*i+0],cells[4*i+1],cells[4*i+2],cells[4*i+3],0,i+1); 1403 ierr = PMMG_Set_triangles(parmesh,bdFaces,bdFaceIds); 1404 ierr = PMMG_Set_metSize(parmesh,MMG5_Vertex,numVertices,MMG5_Tensor); 1405 for (i=0;i<numVertices;i++) PMMG_Set_tensorMet(parmesh,metric[6*i],metric[6*i+1],metric[6*i+2],metric[6*i+3],metric[6*i+4],metric[6*i+5],i+1); 1406 1407 // Interfaces de communication 1408 ierr = PMMG_Set_numberOfNodeCommunicators(parmesh,num_communicators); 1409 printf("nombre de communicateurs sur le proc %d : %d \n",me,num_communicators); 1410 for(ct=0,p=0;p<numProcs;p++) 1411 if (num_per_proc[p] !=0 && p!=me) { 1412 printf("DEBUG %d, taille: %d ct: %d p:%d \n",me,offset[ct+1]-offset[ct],ct,p); 1413 ierr = PMMG_Set_ithNodeCommunicatorSize(parmesh,ct,p,offset[ct+1]-offset[ct]); printf("DEBUG communicator size: %d \n", ierr); 1414 ierr = PMMG_Set_ithNodeCommunicator_nodes(parmesh,ct,&communicators_local_new[offset[ct]],&communicators_global[offset[ct]],1); printf("DEBUG communicator: %d \n", ierr); 1415 1416 // printf("DEBUG (%d) ct: %d\n", me, ct); 1417 // printf("DEBUG(%d) communicators_local_new: ", me); 1418 // for (int iv=offset[ct]; iv<offset[ct+1]; ++iv) { 1419 // printf(" %d,", communicators_local_new[iv]); 1420 // } 1421 // printf("\n"); 1422 // printf("DEBUG(%d) communicators_global: ", me); 1423 // for (int iv=offset[ct]; iv<offset[ct+1]; ++iv) { 1424 // printf(" %d,", communicators_global[iv]); 1425 // } 1426 // printf("\n"); 1427 1428 ct++; 1429 } 1430 1431 // Remaillage 1432 MPI_Barrier(comm); 1433 ierr = PMMG_Set_iparameter(parmesh,PMMG_IPARAM_globalNum,1); 1434 //ierr = PMMG_saveMesh_distributed(parmesh,"mesh_depart"); 1435 ierr = PMMG_parmmglib_distributed(parmesh);printf("DEBUG remaillage //: %d \n", ierr); 1436 //ierr = PMMG_saveMesh_distributed(parmesh,"mesh_apres"); 1437 } 1438 1439 /*3. Passer du nouveau maillage mmg à un maillage lisible par petsc*/ 1440 if (numProcs==1) { 1441 switch(dim) 1442 { 1443 case(2): 1444 numCornersNew = 3; 1445 ierr = MMG2D_Get_meshSize(mmg_mesh,&numVerticesNew,&numCellsNew,0,&numFacesNew); 1446 1447 ierr = PetscCalloc4(2*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 1448 ierr = PetscCalloc3(3*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 1449 ierr = PetscCalloc4(2*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 1450 ierr = MMG2D_Get_vertices(mmg_mesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 1451 ierr = MMG2D_Get_triangles(mmg_mesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 1452 ierr = MMG2D_Get_edges(mmg_mesh,faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces); 1453 1454 for (i=0;i<3*numCellsNew;i++) cellsNew[i]-=1; 1455 for (i=0;i<2*numFacesNew;i++) faces[i]-=1; 1456 1457 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, dmNew);CHKERRQ(ierr); 1458 ierr = MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg_mesh, MMG5_ARG_ppMet, &mmg_metric, MMG5_ARG_end); printf("DEBUG libération mémoire : %d \n", ierr); 1459 break; 1460 case(3): 1461 numCornersNew = 4; 1462 ierr = MMG3D_Get_meshSize(mmg_mesh,&numVerticesNew,&numCellsNew,0,&numFacesNew,0,0); 1463 ierr = PetscCalloc4(3*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 1464 ierr = PetscCalloc3(4*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 1465 ierr = PetscCalloc4(3*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 1466 ierr = MMG3D_Get_vertices(mmg_mesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 1467 ierr = MMG3D_Get_tetrahedra(mmg_mesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 1468 ierr = MMG3D_Get_triangles(mmg_mesh,faces,tab_cl_faces,tab_areRequiredFaces); 1469 1470 for (i=0;i<4*numCellsNew;i++) cellsNew[i]-=1; 1471 for (i=0;i<3*numFacesNew;i++) faces[i]-=1; 1472 1473 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNew, NULL, dmNew);CHKERRQ(ierr); 1474 ierr = MMG3D_Free_all(MMG5_ARG_start,MMG5_ARG_ppMesh,&mmg_mesh,MMG5_ARG_ppMet,&mmg_metric,MMG5_ARG_end); printf("DEBUG libération mémoire : %d \n", ierr); 1475 break; 1476 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 1477 } 1478 } 1479 else { 1480 numCornersNew = 4; 1481 ierr = PMMG_Get_meshSize(parmesh,&numVerticesNew,&numCellsNew,0,&numFacesNew,0,0); 1482 ierr = PetscCalloc4(3*numVerticesNew,&verticesNew,numVerticesNew,&tab_cl_verticies_new,numVerticesNew,&tab_areCorners,numVerticesNew,&tab_areRequiredVerticies);CHKERRQ(ierr); 1483 ierr = PetscCalloc3(4*numCellsNew,&cellsNew,numCellsNew,&tab_cl_cells_new,numCellsNew,&tab_areRequiredCells);CHKERRQ(ierr); 1484 ierr = PetscCalloc4(3*numFacesNew,&faces,numFacesNew,&tab_cl_faces,numFacesNew,&tab_areRidges,numFacesNew,&tab_areRequiredFaces);CHKERRQ(ierr); 1485 ierr = PMMG_Get_vertices(parmesh,verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies); 1486 ierr = PMMG_Get_tetrahedra(parmesh,cellsNew,tab_cl_cells_new,tab_areRequiredCells); 1487 ierr = PMMG_Get_triangles(parmesh,faces,tab_cl_faces,tab_areRequiredFaces); 1488 1489 ierr = PetscCalloc2(numVerticesNew,&gv_new,numVerticesNew,&ranks_own); 1490 ierr = PMMG_Get_verticesGloNum(parmesh,gv_new,ranks_own); 1491 1492 // Décallage et changement de numérotation 1493 for (i=0;i<3*numFacesNew;i++) faces[i]-=1; 1494 for (i=0;i<4*numCellsNew;i++) cellsNew[i]=gv_new[cellsNew[i]-1]-1; 1495 1496 // Calcul de la liste des sommets sur chacun des procs 1497 for (i=0;i<numVerticesNew;i++) if (ranks_own[i]==me) numVerticesNewNew++; 1498 ierr = PetscCalloc1(numVerticesNewNew,&VerticesNewNew); CHKERRQ(ierr); 1499 for (ct=0,i=0;i<numVerticesNew;i++) if (ranks_own[i]==me) { 1500 VerticesNewNew[ct++]=verticesNew[i]; 1501 } 1502 1503 // tests 1504 for (p=0;p<numProcs;p++) { 1505 if (p==me) printf("me : %d , numVerticesNew : %d \n",me,numVerticesNew); 1506 if (p==me) { for (i=0;i<numVerticesNew;i++) printf("%d ",gv_new[i]); printf("\n");} 1507 if (p==me) { for (i=0;i<numVerticesNew;i++) printf("%d ",ranks_own[i]); printf("\n");} 1508 MPI_Barrier(comm); 1509 } 1510 for(i=0;i<4*numCellsNew;i++) if (cellsNew[i]<0) printf("DEBUG -1 %d %d\n",i,cellsNew[i]); 1511 1512 for (p=0;p<numProcs;p++) { 1513 if (p==me) { 1514 for (i=0;i<numVerticesNew;i++) { 1515 printf("DBEUG(%d) VerticesNewNew[%d]: %1.2f %1.2f %1.2f\n", me, i, 1516 VerticesNewNew[3*i], VerticesNewNew[3*i+1], VerticesNewNew[3*i+2]); 1517 } 1518 } 1519 MPI_Barrier(comm); 1520 } 1521 1522 // printf("DEBUG %d nombre de sommets %d, nombre de tétra %d num total: %d \n",me,numVerticesNew,numCellsNew,num_total); 1523 // if (me==1) for (i=0;i<numVerticesNew;i++) printf("%d sommets %d %lf %lf %lf\n",me,i,verticesNew[3*i],verticesNew[3*i+1],verticesNew[3*i+2]); 1524 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, VerticesNewNew, NULL, dmNew);CHKERRQ(ierr); 1525 ierr = PMMG_Free_all(PMMG_ARG_start,PMMG_ARG_ppParMesh,&parmesh,PMMG_ARG_end); 1526 1527 } 1528 1529 /*Reconstruction des conditions aux limites*/ 1530 //ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 1531 //ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 1532 //for(i=0;i<numFacesNew;i++) ierr = DMLabelSetValue(bdLabelNew,faces[i],tab_cl_faces[i]);CHKERRQ(ierr); 1533 1534 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 1535 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 1536 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 1537 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 1538 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 1539 1540 switch(dim) 1541 { 1542 case(2): 1543 for (i=0;i<numFacesNew;i++){ 1544 const PetscInt *supportA, *supportB; 1545 PetscInt SA,SB,inter=-1; 1546 ierr = DMPlexGetSupportSize(*dmNew,faces[2*i]+vStart,&SA); 1547 ierr = DMPlexGetSupportSize(*dmNew,faces[2*i+1]+vStart,&SB); 1548 ierr = DMPlexGetSupport(*dmNew,faces[2*i]+vStart,&supportA); 1549 ierr = DMPlexGetSupport(*dmNew,faces[2*i+1]+vStart,&supportB); 1550 /*printf("Numero des points : %d %d \n",faces[2*i]+1,faces[2*i+1]+1); 1551 printf("Taille des supports : %d %d \n",SA,SB); 1552 printf("Support de A ");for(d=0;d<SA;d++) printf("%d ",supportA[d]); printf("\n"); 1553 printf("Support de B ");for(d=0;d<SB;d++) printf("%d ",supportB[d]); printf("\n");*/ 1554 1555 // Calcul de l'intersection: 1556 for(j=0;j<SA;j++){ 1557 for(k=0;k<SB;k++){ 1558 if(supportA[j]==supportB[k]) inter=supportA[j]; 1559 } 1560 } 1561 ierr = DMLabelSetValue(bdLabelNew,inter,tab_cl_faces[i]);CHKERRQ(ierr); 1562 } 1563 break; 1564 case(3): 1565 for (i=0;i<numFacesNew;i++){ 1566 const PetscInt *supportA, *supportB, *supportC; 1567 const PetscInt *supportiA, *supportiB, *supportiC; 1568 1569 PetscInt SA,SB,SC,inter=-1; 1570 PetscInt SiA,SiB,SiC; 1571 PetscInt ia, ib, ic; 1572 PetscInt ja, jb, jc; 1573 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i]+vStart,&SA); 1574 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i+1]+vStart,&SB); 1575 ierr = DMPlexGetSupportSize(*dmNew,faces[3*i+2]+vStart,&SC); 1576 ierr = DMPlexGetSupport(*dmNew,faces[3*i]+vStart,&supportA); 1577 ierr = DMPlexGetSupport(*dmNew,faces[3*i+1]+vStart,&supportB); 1578 ierr = DMPlexGetSupport(*dmNew,faces[3*i+2]+vStart,&supportC); 1579 1580 //printf("%d %d %d\n",SA,SB,SC); 1581 inter=-1; 1582 for (ia=0;ia<SA;ia++) { 1583 ierr = DMPlexGetSupportSize(*dmNew,supportA[ia],&SiA); 1584 ierr = DMPlexGetSupport(*dmNew,supportA[ia],&supportiA); 1585 for (ib=0;ib<SB;ib++) { 1586 ierr = DMPlexGetSupportSize(*dmNew,supportB[ib],&SiB); 1587 ierr = DMPlexGetSupport(*dmNew,supportB[ib],&supportiB); 1588 for (ic=0;ic<SC;ic++) { 1589 ierr = DMPlexGetSupportSize(*dmNew,supportC[ic],&SiC); 1590 ierr = DMPlexGetSupport(*dmNew,supportC[ic],&supportiC); 1591 for(ja=0;ja<SiA;ja++) { 1592 for(jb=0;jb<SiB;jb++) { 1593 for(jc=0;jc<SiC;jc++) { 1594 if(supportiA[ja]==supportiB[jb] && supportiA[ja]==supportiC[jc]) inter=supportiA[ja]; 1595 } 1596 } 1597 } 1598 } 1599 } 1600 } 1601 ierr = DMLabelSetValue(bdLabelNew,inter,tab_cl_faces[i]);CHKERRQ(ierr); 1602 } 1603 break; 1604 default:SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No MMG adaptation defined for dimension %D", dim); 1605 } 1606 1607 /*libération des tableaux en mémoire TO DO*/ 1608 ierr = PetscFree3(cells,metric,vertices);CHKERRQ(ierr); printf("debug 1 \n"); 1609 ierr = PetscFree4(bdFaces,bdFaceIds,tab_cl_verticies,tab_cl_triangles);CHKERRQ(ierr);printf("debug 2 \n"); 1610 ierr = PetscFree4(verticesNew,tab_cl_verticies_new,tab_areCorners,tab_areRequiredVerticies);CHKERRQ(ierr);printf("debug 3 \n"); 1611 ierr = PetscFree4(faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces);CHKERRQ(ierr);printf("debug 4 \n"); 1612 ierr = PetscFree4(faces,tab_cl_faces,tab_areRidges,tab_areRequiredFaces);CHKERRQ(ierr);printf("debug 5 \n"); 1613 1614 if (numProcs>1) ierr = PetscFree4(communicators_local,communicators_global,num_per_proc,flags_proc); CHKERRQ(ierr); 1615 if (numProcs>1) ierr = PetscFree3(VerticesNewNew,ranks_own,gv_new); CHKERRQ(ierr); 1616 1617 PetscFunctionReturn(0); 1618 #else 1619 PetscFunctionBegin; 1620 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic."); 1621 PetscFunctionReturn(0); 1622 #endif 1623 } 1624 1625 1626 /* 1627 DMAdaptMetric_Plex - Generates a new mesh conforming to a metric field. 1628 1629 Input Parameters: 1630 + dm - The DM object 1631 . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector 1632 - bdLabel - Label for boundary tags which are preserved in dmNew, or NULL. Should not be named "_boundary_". 1633 1634 Output Parameter: 1635 . dmNew - the new DM 1636 1637 Level: advanced 1638 1639 .seealso: DMCoarsen(), DMRefine() 1640 */ 1641 PetscErrorCode DMAdaptMetric_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DM *dmNew) 1642 { 1643 PetscInt remesher = 0; 1644 switch (remesher) { 1645 case 0: 1646 DMAdaptMetricPragmatic_Plex(dm, vertexMetric, bdLabel, dmNew); 1647 break; 1648 case 1: 1649 DMAdaptMetricMMG_Plex(dm, vertexMetric, bdLabel, dmNew); 1650 break; 1651 case 2: 1652 DMAdaptMetricParMMG_Plex(dm, vertexMetric, bdLabel, dmNew); 1653 break; 1654 } 1655 }