1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <../src/sys/utils/hash.h> 3 #include <petsc-private/isimpl.h> 4 #include <petscsf.h> 5 6 /* Logging support */ 7 PetscLogEvent DMPLEX_Interpolate, DMPLEX_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM; 8 9 PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer); 10 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer); 11 PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer); 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "DMPlexGetFieldType_Internal" 15 PetscErrorCode DMPlexGetFieldType_Internal(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft) 16 { 17 PetscInt dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, vdof = 0, cdof = 0; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 *ft = PETSC_VTK_POINT_FIELD; 22 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 23 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 24 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 25 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 26 if (field >= 0) { 27 if ((vStart >= pStart) && (vStart < pEnd)) {ierr = PetscSectionGetFieldDof(section, vStart, field, &vdof);CHKERRQ(ierr);} 28 if ((cStart >= pStart) && (cStart < pEnd)) {ierr = PetscSectionGetFieldDof(section, cStart, field, &cdof);CHKERRQ(ierr);} 29 } else { 30 if ((vStart >= pStart) && (vStart < pEnd)) {ierr = PetscSectionGetDof(section, vStart, &vdof);CHKERRQ(ierr);} 31 if ((cStart >= pStart) && (cStart < pEnd)) {ierr = PetscSectionGetDof(section, cStart, &cdof);CHKERRQ(ierr);} 32 } 33 if (vdof) { 34 *sStart = vStart; 35 *sEnd = vEnd; 36 if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD; 37 else *ft = PETSC_VTK_POINT_FIELD; 38 } else if (cdof) { 39 *sStart = cStart; 40 *sEnd = cEnd; 41 if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD; 42 else *ft = PETSC_VTK_CELL_FIELD; 43 } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK"); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "VecView_Plex_Local" 49 PetscErrorCode VecView_Plex_Local(Vec v, PetscViewer viewer) 50 { 51 DM dm; 52 PetscBool isvtk, ishdf5, isseq; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 57 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM"); 58 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr); 59 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 60 ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr); 61 if (isvtk || ishdf5) { 62 PetscInt numFields; 63 PetscBool fem = PETSC_FALSE; 64 65 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 66 if (numFields) { 67 PetscObject fe; 68 69 ierr = DMGetField(dm, 0, &fe);CHKERRQ(ierr); 70 if (fe->classid == PETSCFE_CLASSID) fem = PETSC_TRUE; 71 } 72 if (fem) {ierr = DMPlexInsertBoundaryValuesFEM(dm, v);CHKERRQ(ierr);} 73 } 74 if (isvtk) { 75 PetscSection section; 76 PetscViewerVTKFieldType ft; 77 PetscInt pStart, pEnd; 78 79 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 80 ierr = DMPlexGetFieldType_Internal(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);CHKERRQ(ierr); 81 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); /* viewer drops reference */ 82 ierr = PetscObjectReference((PetscObject) v);CHKERRQ(ierr); /* viewer drops reference */ 83 ierr = PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);CHKERRQ(ierr); 84 } else if (ishdf5) { 85 #if defined(PETSC_HAVE_HDF5) 86 ierr = VecView_Plex_Local_HDF5(v, viewer);CHKERRQ(ierr); 87 #else 88 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 89 #endif 90 } else { 91 if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);} 92 else {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);} 93 } 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "VecView_Plex" 99 PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer) 100 { 101 DM dm; 102 PetscBool isvtk, ishdf5, isseq; 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 107 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM"); 108 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr); 109 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 110 ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr); 111 if (isvtk) { 112 Vec locv; 113 const char *name; 114 115 ierr = DMGetLocalVector(dm, &locv);CHKERRQ(ierr); 116 ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr); 117 ierr = PetscObjectSetName((PetscObject) locv, name);CHKERRQ(ierr); 118 ierr = DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr); 119 ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr); 120 ierr = VecView_Plex_Local(locv, viewer);CHKERRQ(ierr); 121 ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr); 122 } else if (ishdf5) { 123 #if defined(PETSC_HAVE_HDF5) 124 ierr = VecView_Plex_HDF5(v, viewer);CHKERRQ(ierr); 125 #else 126 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 127 #endif 128 } else { 129 if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);} 130 else {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);} 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "VecLoad_Plex_Local" 137 PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer) 138 { 139 DM dm; 140 PetscBool ishdf5; 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 145 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM"); 146 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 147 if (ishdf5) { 148 DM dmBC; 149 Vec gv; 150 const char *name; 151 152 ierr = DMGetOutputDM(dm, &dmBC);CHKERRQ(ierr); 153 ierr = DMGetGlobalVector(dmBC, &gv);CHKERRQ(ierr); 154 ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr); 155 ierr = PetscObjectSetName((PetscObject) gv, name);CHKERRQ(ierr); 156 ierr = VecLoad_Default(gv, viewer);CHKERRQ(ierr); 157 ierr = DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);CHKERRQ(ierr); 158 ierr = DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);CHKERRQ(ierr); 159 ierr = DMRestoreGlobalVector(dmBC, &gv);CHKERRQ(ierr); 160 } else { 161 ierr = VecLoad_Default(v, viewer);CHKERRQ(ierr); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "VecLoad_Plex" 168 PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer) 169 { 170 DM dm; 171 PetscBool ishdf5; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = VecGetDM(v, &dm);CHKERRQ(ierr); 176 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM"); 177 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 178 if (ishdf5) { 179 #if defined(PETSC_HAVE_HDF5) 180 ierr = VecLoad_Plex_HDF5(v, viewer);CHKERRQ(ierr); 181 #else 182 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 183 #endif 184 } else { 185 ierr = VecLoad_Default(v, viewer);CHKERRQ(ierr); 186 } 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "DMPlexView_Ascii" 192 PetscErrorCode DMPlexView_Ascii(DM dm, PetscViewer viewer) 193 { 194 DM_Plex *mesh = (DM_Plex*) dm->data; 195 DM cdm; 196 DMLabel markers; 197 PetscSection coordSection; 198 Vec coordinates; 199 PetscViewerFormat format; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 204 ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr); 205 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 206 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 207 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 208 const char *name; 209 PetscInt maxConeSize, maxSupportSize; 210 PetscInt pStart, pEnd, p; 211 PetscMPIInt rank, size; 212 213 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 214 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRQ(ierr); 215 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 216 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 217 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 218 ierr = PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);CHKERRQ(ierr); 219 ierr = PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);CHKERRQ(ierr); 220 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);CHKERRQ(ierr); 221 ierr = PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);CHKERRQ(ierr); 222 ierr = PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);CHKERRQ(ierr); 223 for (p = pStart; p < pEnd; ++p) { 224 PetscInt dof, off, s; 225 226 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 227 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 228 for (s = off; s < off+dof; ++s) { 229 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D ----> %D\n", rank, p, mesh->supports[s]);CHKERRQ(ierr); 230 } 231 } 232 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);CHKERRQ(ierr); 234 for (p = pStart; p < pEnd; ++p) { 235 PetscInt dof, off, c; 236 237 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 238 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 239 for (c = off; c < off+dof; ++c) { 240 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);CHKERRQ(ierr); 241 } 242 } 243 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 244 ierr = PetscSectionGetChart(coordSection, &pStart, NULL);CHKERRQ(ierr); 245 if (pStart >= 0) {ierr = PetscSectionVecView(coordSection, coordinates, viewer);CHKERRQ(ierr);} 246 ierr = DMPlexGetLabel(dm, "marker", &markers);CHKERRQ(ierr); 247 ierr = DMLabelView(markers,viewer);CHKERRQ(ierr); 248 if (size > 1) { 249 PetscSF sf; 250 251 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 252 ierr = PetscSFView(sf, viewer);CHKERRQ(ierr); 253 } 254 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 255 } else if (format == PETSC_VIEWER_ASCII_LATEX) { 256 const char *name; 257 const char *colors[3] = {"red", "blue", "green"}; 258 const int numColors = 3; 259 PetscReal scale = 2.0; 260 PetscScalar *coords; 261 PetscInt depth, cStart, cEnd, c, vStart, vEnd, v, eStart = 0, eEnd = 0, e, p; 262 PetscMPIInt rank, size; 263 264 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 265 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 266 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRQ(ierr); 267 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 268 ierr = PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);CHKERRQ(ierr); 269 ierr = PetscViewerASCIIPrintf(viewer, "\ 270 \\documentclass[crop,multi=false]{standalone}\n\n\ 271 \\usepackage{tikz}\n\ 272 \\usepackage{pgflibraryshapes}\n\ 273 \\usetikzlibrary{backgrounds}\n\ 274 \\usetikzlibrary{arrows}\n\ 275 \\begin{document}\n\ 276 \\section{%s}\n\ 277 \\begin{center}\n", name, 8.0/scale);CHKERRQ(ierr); 278 ierr = PetscViewerASCIIPrintf(viewer, "Mesh for process ");CHKERRQ(ierr); 279 for (p = 0; p < size; ++p) { 280 if (p > 0 && p == size-1) { 281 ierr = PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);CHKERRQ(ierr); 282 } else if (p > 0) { 283 ierr = PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);CHKERRQ(ierr); 284 } 285 ierr = PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);CHKERRQ(ierr); 286 } 287 ierr = PetscViewerASCIIPrintf(viewer, ".\n\n\n\ 288 \\begin{tikzpicture}[scale = %g,font=\\fontsize{8}{8}\\selectfont]\n");CHKERRQ(ierr); 289 /* Plot vertices */ 290 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 291 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 292 ierr = PetscViewerASCIIPrintf(viewer, "\\path\n");CHKERRQ(ierr); 293 for (v = vStart; v < vEnd; ++v) { 294 PetscInt off, dof, d; 295 296 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 297 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 298 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "(");CHKERRQ(ierr); 299 for (d = 0; d < dof; ++d) { 300 if (d > 0) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, ",");CHKERRQ(ierr);} 301 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*PetscRealPart(coords[off+d])));CHKERRQ(ierr); 302 } 303 ierr = PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", v, rank, colors[rank%numColors], v);CHKERRQ(ierr); 304 } 305 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 306 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 307 ierr = PetscViewerASCIIPrintf(viewer, "(0,0);\n");CHKERRQ(ierr); 308 /* Plot edges */ 309 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 310 ierr = PetscViewerASCIIPrintf(viewer, "\\path\n");CHKERRQ(ierr); 311 if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);} 312 for (e = eStart; e < eEnd; ++e) { 313 const PetscInt *cone; 314 PetscInt coneSize, offA, offB, dof, d; 315 316 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 317 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Edge %d cone should have two vertices, not %d", e, coneSize); 318 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 319 ierr = PetscSectionGetDof(coordSection, cone[0], &dof);CHKERRQ(ierr); 320 ierr = PetscSectionGetOffset(coordSection, cone[0], &offA);CHKERRQ(ierr); 321 ierr = PetscSectionGetOffset(coordSection, cone[1], &offB);CHKERRQ(ierr); 322 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "(");CHKERRQ(ierr); 323 for (d = 0; d < dof; ++d) { 324 if (d > 0) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, ",");CHKERRQ(ierr);} 325 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)(scale*0.5*PetscRealPart(coords[offA+d]+coords[offB+d])));CHKERRQ(ierr); 326 } 327 ierr = PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", e, rank, colors[rank%numColors], e);CHKERRQ(ierr); 328 } 329 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 330 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPrintf(viewer, "(0,0);\n");CHKERRQ(ierr); 332 /* Plot cells */ 333 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 334 for (c = cStart; c < cEnd; ++c) { 335 PetscInt *closure = NULL; 336 PetscInt closureSize, firstPoint = -1; 337 338 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 339 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);CHKERRQ(ierr); 340 for (p = 0; p < closureSize*2; p += 2) { 341 const PetscInt point = closure[p]; 342 343 if ((point < vStart) || (point >= vEnd)) continue; 344 if (firstPoint >= 0) {ierr = PetscViewerASCIISynchronizedPrintf(viewer, " -- ");CHKERRQ(ierr);} 345 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%D)", point, rank);CHKERRQ(ierr); 346 if (firstPoint < 0) firstPoint = point; 347 } 348 /* Why doesn't this work? ierr = PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n");CHKERRQ(ierr); */ 349 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%D);\n", firstPoint, rank);CHKERRQ(ierr); 350 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 351 } 352 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 353 ierr = PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n\\end{center}\n");CHKERRQ(ierr); 354 ierr = PetscViewerASCIIPrintf(viewer, "\\end{document}\n", name);CHKERRQ(ierr); 355 } else { 356 MPI_Comm comm; 357 PetscInt *sizes, *hybsizes; 358 PetscInt locDepth, depth, dim, d, pMax[4]; 359 PetscInt pStart, pEnd, p; 360 PetscInt numLabels, l; 361 const char *name; 362 PetscMPIInt size; 363 364 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 365 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 366 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 367 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr); 368 if (name) {ierr = PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dim);CHKERRQ(ierr);} 369 else {ierr = PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);CHKERRQ(ierr);} 370 ierr = DMPlexGetDepth(dm, &locDepth);CHKERRQ(ierr); 371 ierr = MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 372 ierr = DMPlexGetHybridBounds(dm, &pMax[depth], depth > 0 ? &pMax[depth-1] : NULL, &pMax[1], &pMax[0]);CHKERRQ(ierr); 373 ierr = PetscMalloc2(size,&sizes,size,&hybsizes);CHKERRQ(ierr); 374 if (depth == 1) { 375 ierr = DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);CHKERRQ(ierr); 376 pEnd = pEnd - pStart; 377 ierr = MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 378 ierr = PetscViewerASCIIPrintf(viewer, " %D-cells:", 0);CHKERRQ(ierr); 379 for (p = 0; p < size; ++p) {ierr = PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);CHKERRQ(ierr);} 380 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 381 ierr = DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd);CHKERRQ(ierr); 382 pEnd = pEnd - pStart; 383 ierr = MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 384 ierr = PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);CHKERRQ(ierr); 385 for (p = 0; p < size; ++p) {ierr = PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);CHKERRQ(ierr);} 386 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 387 } else { 388 for (d = 0; d <= dim; d++) { 389 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 390 pEnd -= pStart; 391 pMax[d] -= pStart; 392 ierr = MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 393 ierr = MPI_Gather(&pMax[d], 1, MPIU_INT, hybsizes, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 394 ierr = PetscViewerASCIIPrintf(viewer, " %D-cells:", d);CHKERRQ(ierr); 395 for (p = 0; p < size; ++p) { 396 if (hybsizes[p] >= 0) {ierr = PetscViewerASCIIPrintf(viewer, " %D (%D)", sizes[p], sizes[p] - hybsizes[p]);CHKERRQ(ierr);} 397 else {ierr = PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);CHKERRQ(ierr);} 398 } 399 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 400 } 401 } 402 ierr = PetscFree2(sizes,hybsizes);CHKERRQ(ierr); 403 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 404 if (numLabels) {ierr = PetscViewerASCIIPrintf(viewer, "Labels:\n");CHKERRQ(ierr);} 405 for (l = 0; l < numLabels; ++l) { 406 DMLabel label; 407 const char *name; 408 IS valueIS; 409 const PetscInt *values; 410 PetscInt numValues, v; 411 412 ierr = DMPlexGetLabelName(dm, l, &name);CHKERRQ(ierr); 413 ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); 414 ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr); 415 ierr = PetscViewerASCIIPrintf(viewer, " %s: %d strata of sizes (", name, numValues);CHKERRQ(ierr); 416 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 417 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 418 for (v = 0; v < numValues; ++v) { 419 PetscInt size; 420 421 ierr = DMLabelGetStratumSize(label, values[v], &size);CHKERRQ(ierr); 422 if (v > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);} 423 ierr = PetscViewerASCIIPrintf(viewer, "%d", size);CHKERRQ(ierr); 424 } 425 ierr = PetscViewerASCIIPrintf(viewer, ")\n");CHKERRQ(ierr); 426 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 427 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 428 } 429 } 430 PetscFunctionReturn(0); 431 } 432 433 #undef __FUNCT__ 434 #define __FUNCT__ "DMView_Plex" 435 PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer) 436 { 437 PetscBool iascii, ishdf5; 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 442 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 443 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 444 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 445 if (iascii) { 446 ierr = DMPlexView_Ascii(dm, viewer);CHKERRQ(ierr); 447 } else if (ishdf5) { 448 #if defined(PETSC_HAVE_HDF5) 449 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);CHKERRQ(ierr); 450 ierr = DMPlexView_HDF5(dm, viewer);CHKERRQ(ierr); 451 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 452 #else 453 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 454 #endif 455 } 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "DMLoad_Plex" 461 PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer) 462 { 463 PetscBool isbinary, ishdf5; 464 PetscErrorCode ierr; 465 466 PetscFunctionBegin; 467 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 468 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 469 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr); 470 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 471 if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");} 472 else if (ishdf5) { 473 #if defined(PETSC_HAVE_HDF5) 474 ierr = DMPlexLoad_HDF5(dm, viewer);CHKERRQ(ierr); 475 #else 476 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 477 #endif 478 } 479 PetscFunctionReturn(0); 480 } 481 482 #undef __FUNCT__ 483 #define __FUNCT__ "BoundaryDestroy" 484 static PetscErrorCode BoundaryDestroy(DMBoundary *boundary) 485 { 486 DMBoundary b, next; 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 if (!boundary) PetscFunctionReturn(0); 491 b = *boundary; 492 *boundary = NULL; 493 for (; b; b = next) { 494 next = b->next; 495 ierr = PetscFree(b->ids);CHKERRQ(ierr); 496 ierr = PetscFree(b->name);CHKERRQ(ierr); 497 ierr = PetscFree(b->labelname);CHKERRQ(ierr); 498 ierr = PetscFree(b);CHKERRQ(ierr); 499 } 500 PetscFunctionReturn(0); 501 } 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "DMDestroy_Plex" 505 PetscErrorCode DMDestroy_Plex(DM dm) 506 { 507 DM_Plex *mesh = (DM_Plex*) dm->data; 508 DMLabel next = mesh->labels; 509 PetscErrorCode ierr; 510 511 PetscFunctionBegin; 512 if (--mesh->refct > 0) PetscFunctionReturn(0); 513 ierr = PetscSectionDestroy(&mesh->coneSection);CHKERRQ(ierr); 514 ierr = PetscFree(mesh->cones);CHKERRQ(ierr); 515 ierr = PetscFree(mesh->coneOrientations);CHKERRQ(ierr); 516 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr); 517 ierr = PetscFree(mesh->supports);CHKERRQ(ierr); 518 ierr = PetscFree(mesh->facesTmp);CHKERRQ(ierr); 519 ierr = PetscFree(mesh->tetgenOpts);CHKERRQ(ierr); 520 ierr = PetscFree(mesh->triangleOpts);CHKERRQ(ierr); 521 while (next) { 522 DMLabel tmp = next->next; 523 524 ierr = DMLabelDestroy(&next);CHKERRQ(ierr); 525 next = tmp; 526 } 527 ierr = DMDestroy(&mesh->coarseMesh);CHKERRQ(ierr); 528 ierr = DMLabelDestroy(&mesh->subpointMap);CHKERRQ(ierr); 529 ierr = ISDestroy(&mesh->globalVertexNumbers);CHKERRQ(ierr); 530 ierr = ISDestroy(&mesh->globalCellNumbers);CHKERRQ(ierr); 531 ierr = BoundaryDestroy(&mesh->boundary);CHKERRQ(ierr); 532 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */ 533 ierr = PetscFree(mesh);CHKERRQ(ierr); 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "DMCreateMatrix_Plex" 539 PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J) 540 { 541 PetscSection section, sectionGlobal; 542 PetscInt bs = -1; 543 PetscInt localSize; 544 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock; 545 PetscErrorCode ierr; 546 MatType mtype; 547 548 PetscFunctionBegin; 549 ierr = MatInitializePackage();CHKERRQ(ierr); 550 mtype = dm->mattype; 551 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 552 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 553 /* ierr = PetscSectionGetStorageSize(sectionGlobal, &localSize);CHKERRQ(ierr); */ 554 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);CHKERRQ(ierr); 555 ierr = MatCreate(PetscObjectComm((PetscObject)dm), J);CHKERRQ(ierr); 556 ierr = MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 557 ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 558 ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 559 ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr); 560 ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr); 561 ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr); 562 ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr); 563 ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr); 564 ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr); 565 ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr); 566 if (!isShell) { 567 PetscBool fillMatrix = (PetscBool) !dm->prealloc_only; 568 PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal, bsMax, bsMin; 569 570 if (bs < 0) { 571 if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) { 572 PetscInt pStart, pEnd, p, dof, cdof; 573 574 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 575 for (p = pStart; p < pEnd; ++p) { 576 ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr); 577 ierr = PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);CHKERRQ(ierr); 578 if (dof-cdof) { 579 if (bs < 0) { 580 bs = dof-cdof; 581 } else if (bs != dof-cdof) { 582 /* Layout does not admit a pointwise block size */ 583 bs = 1; 584 break; 585 } 586 } 587 } 588 /* Must have same blocksize on all procs (some might have no points) */ 589 bsLocal = bs; 590 ierr = MPI_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 591 bsLocal = bs < 0 ? bsMax : bs; 592 ierr = MPI_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 593 if (bsMin != bsMax) { 594 bs = 1; 595 } else { 596 bs = bsMax; 597 } 598 } else { 599 bs = 1; 600 } 601 } 602 ierr = PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);CHKERRQ(ierr); 603 ierr = DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);CHKERRQ(ierr); 604 ierr = PetscFree4(dnz, onz, dnzu, onzu);CHKERRQ(ierr); 605 } 606 PetscFunctionReturn(0); 607 } 608 609 #undef __FUNCT__ 610 #define __FUNCT__ "DMPlexGetChart" 611 /*@ 612 DMPlexGetChart - Return the interval for all mesh points [pStart, pEnd) 613 614 Not collective 615 616 Input Parameter: 617 . mesh - The DMPlex 618 619 Output Parameters: 620 + pStart - The first mesh point 621 - pEnd - The upper bound for mesh points 622 623 Level: beginner 624 625 .seealso: DMPlexCreate(), DMPlexSetChart() 626 @*/ 627 PetscErrorCode DMPlexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd) 628 { 629 DM_Plex *mesh = (DM_Plex*) dm->data; 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 634 ierr = PetscSectionGetChart(mesh->coneSection, pStart, pEnd);CHKERRQ(ierr); 635 PetscFunctionReturn(0); 636 } 637 638 #undef __FUNCT__ 639 #define __FUNCT__ "DMPlexSetChart" 640 /*@ 641 DMPlexSetChart - Set the interval for all mesh points [pStart, pEnd) 642 643 Not collective 644 645 Input Parameters: 646 + mesh - The DMPlex 647 . pStart - The first mesh point 648 - pEnd - The upper bound for mesh points 649 650 Output Parameters: 651 652 Level: beginner 653 654 .seealso: DMPlexCreate(), DMPlexGetChart() 655 @*/ 656 PetscErrorCode DMPlexSetChart(DM dm, PetscInt pStart, PetscInt pEnd) 657 { 658 DM_Plex *mesh = (DM_Plex*) dm->data; 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 663 ierr = PetscSectionSetChart(mesh->coneSection, pStart, pEnd);CHKERRQ(ierr); 664 ierr = PetscSectionSetChart(mesh->supportSection, pStart, pEnd);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 #undef __FUNCT__ 669 #define __FUNCT__ "DMPlexGetConeSize" 670 /*@ 671 DMPlexGetConeSize - Return the number of in-edges for this point in the Sieve DAG 672 673 Not collective 674 675 Input Parameters: 676 + mesh - The DMPlex 677 - p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 678 679 Output Parameter: 680 . size - The cone size for point p 681 682 Level: beginner 683 684 .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart() 685 @*/ 686 PetscErrorCode DMPlexGetConeSize(DM dm, PetscInt p, PetscInt *size) 687 { 688 DM_Plex *mesh = (DM_Plex*) dm->data; 689 PetscErrorCode ierr; 690 691 PetscFunctionBegin; 692 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 693 PetscValidPointer(size, 3); 694 ierr = PetscSectionGetDof(mesh->coneSection, p, size);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 #undef __FUNCT__ 699 #define __FUNCT__ "DMPlexSetConeSize" 700 /*@ 701 DMPlexSetConeSize - Set the number of in-edges for this point in the Sieve DAG 702 703 Not collective 704 705 Input Parameters: 706 + mesh - The DMPlex 707 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 708 - size - The cone size for point p 709 710 Output Parameter: 711 712 Note: 713 This should be called after DMPlexSetChart(). 714 715 Level: beginner 716 717 .seealso: DMPlexCreate(), DMPlexGetConeSize(), DMPlexSetChart() 718 @*/ 719 PetscErrorCode DMPlexSetConeSize(DM dm, PetscInt p, PetscInt size) 720 { 721 DM_Plex *mesh = (DM_Plex*) dm->data; 722 PetscErrorCode ierr; 723 724 PetscFunctionBegin; 725 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 726 ierr = PetscSectionSetDof(mesh->coneSection, p, size);CHKERRQ(ierr); 727 728 mesh->maxConeSize = PetscMax(mesh->maxConeSize, size); 729 PetscFunctionReturn(0); 730 } 731 732 #undef __FUNCT__ 733 #define __FUNCT__ "DMPlexAddConeSize" 734 /*@ 735 DMPlexAddConeSize - Add the given number of in-edges to this point in the Sieve DAG 736 737 Not collective 738 739 Input Parameters: 740 + mesh - The DMPlex 741 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 742 - size - The additional cone size for point p 743 744 Output Parameter: 745 746 Note: 747 This should be called after DMPlexSetChart(). 748 749 Level: beginner 750 751 .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexGetConeSize(), DMPlexSetChart() 752 @*/ 753 PetscErrorCode DMPlexAddConeSize(DM dm, PetscInt p, PetscInt size) 754 { 755 DM_Plex *mesh = (DM_Plex*) dm->data; 756 PetscInt csize; 757 PetscErrorCode ierr; 758 759 PetscFunctionBegin; 760 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 761 ierr = PetscSectionAddDof(mesh->coneSection, p, size);CHKERRQ(ierr); 762 ierr = PetscSectionGetDof(mesh->coneSection, p, &csize);CHKERRQ(ierr); 763 764 mesh->maxConeSize = PetscMax(mesh->maxConeSize, csize); 765 PetscFunctionReturn(0); 766 } 767 768 #undef __FUNCT__ 769 #define __FUNCT__ "DMPlexGetCone" 770 /*@C 771 DMPlexGetCone - Return the points on the in-edges for this point in the Sieve DAG 772 773 Not collective 774 775 Input Parameters: 776 + mesh - The DMPlex 777 - p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 778 779 Output Parameter: 780 . cone - An array of points which are on the in-edges for point p 781 782 Level: beginner 783 784 Fortran Notes: 785 Since it returns an array, this routine is only available in Fortran 90, and you must 786 include petsc.h90 in your code. 787 788 You must also call DMPlexRestoreCone() after you finish using the returned array. 789 790 .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart() 791 @*/ 792 PetscErrorCode DMPlexGetCone(DM dm, PetscInt p, const PetscInt *cone[]) 793 { 794 DM_Plex *mesh = (DM_Plex*) dm->data; 795 PetscInt off; 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 800 PetscValidPointer(cone, 3); 801 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 802 *cone = &mesh->cones[off]; 803 PetscFunctionReturn(0); 804 } 805 806 #undef __FUNCT__ 807 #define __FUNCT__ "DMPlexSetCone" 808 /*@ 809 DMPlexSetCone - Set the points on the in-edges for this point in the Sieve DAG 810 811 Not collective 812 813 Input Parameters: 814 + mesh - The DMPlex 815 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 816 - cone - An array of points which are on the in-edges for point p 817 818 Output Parameter: 819 820 Note: 821 This should be called after all calls to DMPlexSetConeSize() and DMSetUp(). 822 823 Level: beginner 824 825 .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp() 826 @*/ 827 PetscErrorCode DMPlexSetCone(DM dm, PetscInt p, const PetscInt cone[]) 828 { 829 DM_Plex *mesh = (DM_Plex*) dm->data; 830 PetscInt pStart, pEnd; 831 PetscInt dof, off, c; 832 PetscErrorCode ierr; 833 834 PetscFunctionBegin; 835 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 836 ierr = PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 837 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 838 if (dof) PetscValidPointer(cone, 3); 839 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 840 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 841 for (c = 0; c < dof; ++c) { 842 if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", cone[c], pStart, pEnd); 843 mesh->cones[off+c] = cone[c]; 844 } 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "DMPlexGetConeOrientation" 850 /*@C 851 DMPlexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG 852 853 Not collective 854 855 Input Parameters: 856 + mesh - The DMPlex 857 - p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 858 859 Output Parameter: 860 . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an 861 integer giving the prescription for cone traversal. If it is negative, the cone is 862 traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives 863 the index of the cone point on which to start. 864 865 Level: beginner 866 867 Fortran Notes: 868 Since it returns an array, this routine is only available in Fortran 90, and you must 869 include petsc.h90 in your code. 870 871 You must also call DMPlexRestoreConeOrientation() after you finish using the returned array. 872 873 .seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart() 874 @*/ 875 PetscErrorCode DMPlexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[]) 876 { 877 DM_Plex *mesh = (DM_Plex*) dm->data; 878 PetscInt off; 879 PetscErrorCode ierr; 880 881 PetscFunctionBegin; 882 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 883 #if defined(PETSC_USE_DEBUG) 884 { 885 PetscInt dof; 886 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 887 if (dof) PetscValidPointer(coneOrientation, 3); 888 } 889 #endif 890 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 891 892 *coneOrientation = &mesh->coneOrientations[off]; 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNCT__ 897 #define __FUNCT__ "DMPlexSetConeOrientation" 898 /*@ 899 DMPlexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG 900 901 Not collective 902 903 Input Parameters: 904 + mesh - The DMPlex 905 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 906 - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an 907 integer giving the prescription for cone traversal. If it is negative, the cone is 908 traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives 909 the index of the cone point on which to start. 910 911 Output Parameter: 912 913 Note: 914 This should be called after all calls to DMPlexSetConeSize() and DMSetUp(). 915 916 Level: beginner 917 918 .seealso: DMPlexCreate(), DMPlexGetConeOrientation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp() 919 @*/ 920 PetscErrorCode DMPlexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[]) 921 { 922 DM_Plex *mesh = (DM_Plex*) dm->data; 923 PetscInt pStart, pEnd; 924 PetscInt dof, off, c; 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 929 ierr = PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 930 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 931 if (dof) PetscValidPointer(coneOrientation, 3); 932 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 933 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 934 for (c = 0; c < dof; ++c) { 935 PetscInt cdof, o = coneOrientation[c]; 936 937 ierr = PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);CHKERRQ(ierr); 938 if (o && ((o < -(cdof+1)) || (o >= cdof))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone orientation %D is not in the valid range [%D. %D)", o, -(cdof+1), cdof); 939 mesh->coneOrientations[off+c] = o; 940 } 941 PetscFunctionReturn(0); 942 } 943 944 #undef __FUNCT__ 945 #define __FUNCT__ "DMPlexInsertCone" 946 PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint) 947 { 948 DM_Plex *mesh = (DM_Plex*) dm->data; 949 PetscInt pStart, pEnd; 950 PetscInt dof, off; 951 PetscErrorCode ierr; 952 953 PetscFunctionBegin; 954 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 955 ierr = PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 956 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 957 if ((conePoint < pStart) || (conePoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", conePoint, pStart, pEnd); 958 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 959 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 960 if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof); 961 mesh->cones[off+conePos] = conePoint; 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "DMPlexInsertConeOrientation" 967 PetscErrorCode DMPlexInsertConeOrientation(DM dm, PetscInt p, PetscInt conePos, PetscInt coneOrientation) 968 { 969 DM_Plex *mesh = (DM_Plex*) dm->data; 970 PetscInt pStart, pEnd; 971 PetscInt dof, off; 972 PetscErrorCode ierr; 973 974 PetscFunctionBegin; 975 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 976 ierr = PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 977 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 978 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 979 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 980 if ((conePos < 0) || (conePos >= dof)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D of point %D is not in the valid range [0, %D)", conePos, p, dof); 981 mesh->coneOrientations[off+conePos] = coneOrientation; 982 PetscFunctionReturn(0); 983 } 984 985 #undef __FUNCT__ 986 #define __FUNCT__ "DMPlexGetSupportSize" 987 /*@ 988 DMPlexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG 989 990 Not collective 991 992 Input Parameters: 993 + mesh - The DMPlex 994 - p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 995 996 Output Parameter: 997 . size - The support size for point p 998 999 Level: beginner 1000 1001 .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart(), DMPlexGetConeSize() 1002 @*/ 1003 PetscErrorCode DMPlexGetSupportSize(DM dm, PetscInt p, PetscInt *size) 1004 { 1005 DM_Plex *mesh = (DM_Plex*) dm->data; 1006 PetscErrorCode ierr; 1007 1008 PetscFunctionBegin; 1009 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1010 PetscValidPointer(size, 3); 1011 ierr = PetscSectionGetDof(mesh->supportSection, p, size);CHKERRQ(ierr); 1012 PetscFunctionReturn(0); 1013 } 1014 1015 #undef __FUNCT__ 1016 #define __FUNCT__ "DMPlexSetSupportSize" 1017 /*@ 1018 DMPlexSetSupportSize - Set the number of out-edges for this point in the Sieve DAG 1019 1020 Not collective 1021 1022 Input Parameters: 1023 + mesh - The DMPlex 1024 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1025 - size - The support size for point p 1026 1027 Output Parameter: 1028 1029 Note: 1030 This should be called after DMPlexSetChart(). 1031 1032 Level: beginner 1033 1034 .seealso: DMPlexCreate(), DMPlexGetSupportSize(), DMPlexSetChart() 1035 @*/ 1036 PetscErrorCode DMPlexSetSupportSize(DM dm, PetscInt p, PetscInt size) 1037 { 1038 DM_Plex *mesh = (DM_Plex*) dm->data; 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1043 ierr = PetscSectionSetDof(mesh->supportSection, p, size);CHKERRQ(ierr); 1044 1045 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, size); 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "DMPlexGetSupport" 1051 /*@C 1052 DMPlexGetSupport - Return the points on the out-edges for this point in the Sieve DAG 1053 1054 Not collective 1055 1056 Input Parameters: 1057 + mesh - The DMPlex 1058 - p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1059 1060 Output Parameter: 1061 . support - An array of points which are on the out-edges for point p 1062 1063 Level: beginner 1064 1065 Fortran Notes: 1066 Since it returns an array, this routine is only available in Fortran 90, and you must 1067 include petsc.h90 in your code. 1068 1069 You must also call DMPlexRestoreSupport() after you finish using the returned array. 1070 1071 .seealso: DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone() 1072 @*/ 1073 PetscErrorCode DMPlexGetSupport(DM dm, PetscInt p, const PetscInt *support[]) 1074 { 1075 DM_Plex *mesh = (DM_Plex*) dm->data; 1076 PetscInt off; 1077 PetscErrorCode ierr; 1078 1079 PetscFunctionBegin; 1080 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1081 PetscValidPointer(support, 3); 1082 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 1083 *support = &mesh->supports[off]; 1084 PetscFunctionReturn(0); 1085 } 1086 1087 #undef __FUNCT__ 1088 #define __FUNCT__ "DMPlexSetSupport" 1089 /*@ 1090 DMPlexSetSupport - Set the points on the out-edges for this point in the Sieve DAG 1091 1092 Not collective 1093 1094 Input Parameters: 1095 + mesh - The DMPlex 1096 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1097 - support - An array of points which are on the in-edges for point p 1098 1099 Output Parameter: 1100 1101 Note: 1102 This should be called after all calls to DMPlexSetSupportSize() and DMSetUp(). 1103 1104 Level: beginner 1105 1106 .seealso: DMPlexCreate(), DMPlexGetSupport(), DMPlexSetChart(), DMPlexSetSupportSize(), DMSetUp() 1107 @*/ 1108 PetscErrorCode DMPlexSetSupport(DM dm, PetscInt p, const PetscInt support[]) 1109 { 1110 DM_Plex *mesh = (DM_Plex*) dm->data; 1111 PetscInt pStart, pEnd; 1112 PetscInt dof, off, c; 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1117 ierr = PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);CHKERRQ(ierr); 1118 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 1119 if (dof) PetscValidPointer(support, 3); 1120 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 1121 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 1122 for (c = 0; c < dof; ++c) { 1123 if ((support[c] < pStart) || (support[c] >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", support[c], pStart, pEnd); 1124 mesh->supports[off+c] = support[c]; 1125 } 1126 PetscFunctionReturn(0); 1127 } 1128 1129 #undef __FUNCT__ 1130 #define __FUNCT__ "DMPlexInsertSupport" 1131 PetscErrorCode DMPlexInsertSupport(DM dm, PetscInt p, PetscInt supportPos, PetscInt supportPoint) 1132 { 1133 DM_Plex *mesh = (DM_Plex*) dm->data; 1134 PetscInt pStart, pEnd; 1135 PetscInt dof, off; 1136 PetscErrorCode ierr; 1137 1138 PetscFunctionBegin; 1139 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1140 ierr = PetscSectionGetChart(mesh->supportSection, &pStart, &pEnd);CHKERRQ(ierr); 1141 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 1142 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 1143 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd); 1144 if ((supportPoint < pStart) || (supportPoint >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support point %D is not in the valid range [%D, %D)", supportPoint, pStart, pEnd); 1145 if (supportPos >= dof) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Support position %D of point %D is not in the valid range [0, %D)", supportPos, p, dof); 1146 mesh->supports[off+supportPos] = supportPoint; 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "DMPlexGetTransitiveClosure" 1152 /*@C 1153 DMPlexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG 1154 1155 Not collective 1156 1157 Input Parameters: 1158 + mesh - The DMPlex 1159 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1160 . useCone - PETSC_TRUE for in-edges, otherwise use out-edges 1161 - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used 1162 1163 Output Parameters: 1164 + numPoints - The number of points in the closure, so points[] is of size 2*numPoints 1165 - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...] 1166 1167 Note: 1168 If using internal storage (points is NULL on input), each call overwrites the last output. 1169 1170 Fortran Notes: 1171 Since it returns an array, this routine is only available in Fortran 90, and you must 1172 include petsc.h90 in your code. 1173 1174 The numPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1175 1176 Level: beginner 1177 1178 .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone() 1179 @*/ 1180 PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[]) 1181 { 1182 DM_Plex *mesh = (DM_Plex*) dm->data; 1183 PetscInt *closure, *fifo; 1184 const PetscInt *tmp = NULL, *tmpO = NULL; 1185 PetscInt tmpSize, t; 1186 PetscInt depth = 0, maxSize; 1187 PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0; 1188 PetscErrorCode ierr; 1189 1190 PetscFunctionBegin; 1191 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1192 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1193 /* This is only 1-level */ 1194 if (useCone) { 1195 ierr = DMPlexGetConeSize(dm, p, &tmpSize);CHKERRQ(ierr); 1196 ierr = DMPlexGetCone(dm, p, &tmp);CHKERRQ(ierr); 1197 ierr = DMPlexGetConeOrientation(dm, p, &tmpO);CHKERRQ(ierr); 1198 } else { 1199 ierr = DMPlexGetSupportSize(dm, p, &tmpSize);CHKERRQ(ierr); 1200 ierr = DMPlexGetSupport(dm, p, &tmp);CHKERRQ(ierr); 1201 } 1202 if (depth == 1) { 1203 if (*points) { 1204 closure = *points; 1205 } else { 1206 maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1); 1207 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr); 1208 } 1209 closure[0] = p; closure[1] = 0; 1210 for (t = 0; t < tmpSize; ++t, closureSize += 2) { 1211 closure[closureSize] = tmp[t]; 1212 closure[closureSize+1] = tmpO ? tmpO[t] : 0; 1213 } 1214 if (numPoints) *numPoints = closureSize/2; 1215 if (points) *points = closure; 1216 PetscFunctionReturn(0); 1217 } 1218 maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1); 1219 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr); 1220 if (*points) { 1221 closure = *points; 1222 } else { 1223 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr); 1224 } 1225 closure[0] = p; closure[1] = 0; 1226 for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) { 1227 const PetscInt cp = tmp[t]; 1228 const PetscInt co = tmpO ? tmpO[t] : 0; 1229 1230 closure[closureSize] = cp; 1231 closure[closureSize+1] = co; 1232 fifo[fifoSize] = cp; 1233 fifo[fifoSize+1] = co; 1234 } 1235 /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */ 1236 while (fifoSize - fifoStart) { 1237 const PetscInt q = fifo[fifoStart]; 1238 const PetscInt o = fifo[fifoStart+1]; 1239 const PetscInt rev = o >= 0 ? 0 : 1; 1240 const PetscInt off = rev ? -(o+1) : o; 1241 1242 if (useCone) { 1243 ierr = DMPlexGetConeSize(dm, q, &tmpSize);CHKERRQ(ierr); 1244 ierr = DMPlexGetCone(dm, q, &tmp);CHKERRQ(ierr); 1245 ierr = DMPlexGetConeOrientation(dm, q, &tmpO);CHKERRQ(ierr); 1246 } else { 1247 ierr = DMPlexGetSupportSize(dm, q, &tmpSize);CHKERRQ(ierr); 1248 ierr = DMPlexGetSupport(dm, q, &tmp);CHKERRQ(ierr); 1249 tmpO = NULL; 1250 } 1251 for (t = 0; t < tmpSize; ++t) { 1252 const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize; 1253 const PetscInt cp = tmp[i]; 1254 /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */ 1255 /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1) 1256 const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */ 1257 PetscInt co = tmpO ? tmpO[i] : 0; 1258 PetscInt c; 1259 1260 if (rev) { 1261 PetscInt childSize, coff; 1262 ierr = DMPlexGetConeSize(dm, cp, &childSize);CHKERRQ(ierr); 1263 coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i]; 1264 co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0; 1265 } 1266 /* Check for duplicate */ 1267 for (c = 0; c < closureSize; c += 2) { 1268 if (closure[c] == cp) break; 1269 } 1270 if (c == closureSize) { 1271 closure[closureSize] = cp; 1272 closure[closureSize+1] = co; 1273 fifo[fifoSize] = cp; 1274 fifo[fifoSize+1] = co; 1275 closureSize += 2; 1276 fifoSize += 2; 1277 } 1278 } 1279 fifoStart += 2; 1280 } 1281 if (numPoints) *numPoints = closureSize/2; 1282 if (points) *points = closure; 1283 ierr = DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr); 1284 PetscFunctionReturn(0); 1285 } 1286 1287 #undef __FUNCT__ 1288 #define __FUNCT__ "DMPlexGetTransitiveClosure_Internal" 1289 /*@C 1290 DMPlexGetTransitiveClosure_Internal - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG with a specified initial orientation 1291 1292 Not collective 1293 1294 Input Parameters: 1295 + mesh - The DMPlex 1296 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1297 . orientation - The orientation of the point 1298 . useCone - PETSC_TRUE for in-edges, otherwise use out-edges 1299 - points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used 1300 1301 Output Parameters: 1302 + numPoints - The number of points in the closure, so points[] is of size 2*numPoints 1303 - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...] 1304 1305 Note: 1306 If using internal storage (points is NULL on input), each call overwrites the last output. 1307 1308 Fortran Notes: 1309 Since it returns an array, this routine is only available in Fortran 90, and you must 1310 include petsc.h90 in your code. 1311 1312 The numPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1313 1314 Level: beginner 1315 1316 .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone() 1317 @*/ 1318 PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[]) 1319 { 1320 DM_Plex *mesh = (DM_Plex*) dm->data; 1321 PetscInt *closure, *fifo; 1322 const PetscInt *tmp = NULL, *tmpO = NULL; 1323 PetscInt tmpSize, t; 1324 PetscInt depth = 0, maxSize; 1325 PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0; 1326 PetscErrorCode ierr; 1327 1328 PetscFunctionBegin; 1329 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1330 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1331 /* This is only 1-level */ 1332 if (useCone) { 1333 ierr = DMPlexGetConeSize(dm, p, &tmpSize);CHKERRQ(ierr); 1334 ierr = DMPlexGetCone(dm, p, &tmp);CHKERRQ(ierr); 1335 ierr = DMPlexGetConeOrientation(dm, p, &tmpO);CHKERRQ(ierr); 1336 } else { 1337 ierr = DMPlexGetSupportSize(dm, p, &tmpSize);CHKERRQ(ierr); 1338 ierr = DMPlexGetSupport(dm, p, &tmp);CHKERRQ(ierr); 1339 } 1340 if (depth == 1) { 1341 if (*points) { 1342 closure = *points; 1343 } else { 1344 maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1); 1345 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr); 1346 } 1347 closure[0] = p; closure[1] = ornt; 1348 for (t = 0; t < tmpSize; ++t, closureSize += 2) { 1349 const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize; 1350 closure[closureSize] = tmp[i]; 1351 closure[closureSize+1] = tmpO ? tmpO[i] : 0; 1352 } 1353 if (numPoints) *numPoints = closureSize/2; 1354 if (points) *points = closure; 1355 PetscFunctionReturn(0); 1356 } 1357 maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1); 1358 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr); 1359 if (*points) { 1360 closure = *points; 1361 } else { 1362 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr); 1363 } 1364 closure[0] = p; closure[1] = ornt; 1365 for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) { 1366 const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize; 1367 const PetscInt cp = tmp[i]; 1368 PetscInt co = tmpO ? tmpO[i] : 0; 1369 1370 if (ornt < 0) { 1371 PetscInt childSize, coff; 1372 ierr = DMPlexGetConeSize(dm, cp, &childSize);CHKERRQ(ierr); 1373 coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i]; 1374 co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0; 1375 } 1376 closure[closureSize] = cp; 1377 closure[closureSize+1] = co; 1378 fifo[fifoSize] = cp; 1379 fifo[fifoSize+1] = co; 1380 } 1381 /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */ 1382 while (fifoSize - fifoStart) { 1383 const PetscInt q = fifo[fifoStart]; 1384 const PetscInt o = fifo[fifoStart+1]; 1385 const PetscInt rev = o >= 0 ? 0 : 1; 1386 const PetscInt off = rev ? -(o+1) : o; 1387 1388 if (useCone) { 1389 ierr = DMPlexGetConeSize(dm, q, &tmpSize);CHKERRQ(ierr); 1390 ierr = DMPlexGetCone(dm, q, &tmp);CHKERRQ(ierr); 1391 ierr = DMPlexGetConeOrientation(dm, q, &tmpO);CHKERRQ(ierr); 1392 } else { 1393 ierr = DMPlexGetSupportSize(dm, q, &tmpSize);CHKERRQ(ierr); 1394 ierr = DMPlexGetSupport(dm, q, &tmp);CHKERRQ(ierr); 1395 tmpO = NULL; 1396 } 1397 for (t = 0; t < tmpSize; ++t) { 1398 const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize; 1399 const PetscInt cp = tmp[i]; 1400 /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */ 1401 /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1) 1402 const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */ 1403 PetscInt co = tmpO ? tmpO[i] : 0; 1404 PetscInt c; 1405 1406 if (rev) { 1407 PetscInt childSize, coff; 1408 ierr = DMPlexGetConeSize(dm, cp, &childSize);CHKERRQ(ierr); 1409 coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i]; 1410 co = childSize ? -(((coff+childSize-1)%childSize)+1) : 0; 1411 } 1412 /* Check for duplicate */ 1413 for (c = 0; c < closureSize; c += 2) { 1414 if (closure[c] == cp) break; 1415 } 1416 if (c == closureSize) { 1417 closure[closureSize] = cp; 1418 closure[closureSize+1] = co; 1419 fifo[fifoSize] = cp; 1420 fifo[fifoSize+1] = co; 1421 closureSize += 2; 1422 fifoSize += 2; 1423 } 1424 } 1425 fifoStart += 2; 1426 } 1427 if (numPoints) *numPoints = closureSize/2; 1428 if (points) *points = closure; 1429 ierr = DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr); 1430 PetscFunctionReturn(0); 1431 } 1432 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "DMPlexRestoreTransitiveClosure" 1435 /*@C 1436 DMPlexRestoreTransitiveClosure - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG 1437 1438 Not collective 1439 1440 Input Parameters: 1441 + mesh - The DMPlex 1442 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart() 1443 . useCone - PETSC_TRUE for in-edges, otherwise use out-edges 1444 . numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit 1445 - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit 1446 1447 Note: 1448 If not using internal storage (points is not NULL on input), this call is unnecessary 1449 1450 Fortran Notes: 1451 Since it returns an array, this routine is only available in Fortran 90, and you must 1452 include petsc.h90 in your code. 1453 1454 The numPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1455 1456 Level: beginner 1457 1458 .seealso: DMPlexGetTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone() 1459 @*/ 1460 PetscErrorCode DMPlexRestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[]) 1461 { 1462 PetscErrorCode ierr; 1463 1464 PetscFunctionBegin; 1465 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1466 if (numPoints) PetscValidIntPointer(numPoints,4); 1467 if (points) PetscValidPointer(points,5); 1468 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, points);CHKERRQ(ierr); 1469 if (numPoints) *numPoints = 0; 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "DMPlexGetMaxSizes" 1475 /*@ 1476 DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG 1477 1478 Not collective 1479 1480 Input Parameter: 1481 . mesh - The DMPlex 1482 1483 Output Parameters: 1484 + maxConeSize - The maximum number of in-edges 1485 - maxSupportSize - The maximum number of out-edges 1486 1487 Level: beginner 1488 1489 .seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart() 1490 @*/ 1491 PetscErrorCode DMPlexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize) 1492 { 1493 DM_Plex *mesh = (DM_Plex*) dm->data; 1494 1495 PetscFunctionBegin; 1496 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1497 if (maxConeSize) *maxConeSize = mesh->maxConeSize; 1498 if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize; 1499 PetscFunctionReturn(0); 1500 } 1501 1502 #undef __FUNCT__ 1503 #define __FUNCT__ "DMSetUp_Plex" 1504 PetscErrorCode DMSetUp_Plex(DM dm) 1505 { 1506 DM_Plex *mesh = (DM_Plex*) dm->data; 1507 PetscInt size; 1508 PetscErrorCode ierr; 1509 1510 PetscFunctionBegin; 1511 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1512 ierr = PetscSectionSetUp(mesh->coneSection);CHKERRQ(ierr); 1513 ierr = PetscSectionGetStorageSize(mesh->coneSection, &size);CHKERRQ(ierr); 1514 ierr = PetscMalloc1(size, &mesh->cones);CHKERRQ(ierr); 1515 ierr = PetscCalloc1(size, &mesh->coneOrientations);CHKERRQ(ierr); 1516 if (mesh->maxSupportSize) { 1517 ierr = PetscSectionSetUp(mesh->supportSection);CHKERRQ(ierr); 1518 ierr = PetscSectionGetStorageSize(mesh->supportSection, &size);CHKERRQ(ierr); 1519 ierr = PetscMalloc1(size, &mesh->supports);CHKERRQ(ierr); 1520 } 1521 PetscFunctionReturn(0); 1522 } 1523 1524 #undef __FUNCT__ 1525 #define __FUNCT__ "DMCreateSubDM_Plex" 1526 PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 1527 { 1528 PetscErrorCode ierr; 1529 1530 PetscFunctionBegin; 1531 if (subdm) {ierr = DMClone(dm, subdm);CHKERRQ(ierr);} 1532 ierr = DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);CHKERRQ(ierr); 1533 PetscFunctionReturn(0); 1534 } 1535 1536 #undef __FUNCT__ 1537 #define __FUNCT__ "DMPlexSymmetrize" 1538 /*@ 1539 DMPlexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation 1540 1541 Not collective 1542 1543 Input Parameter: 1544 . mesh - The DMPlex 1545 1546 Output Parameter: 1547 1548 Note: 1549 This should be called after all calls to DMPlexSetCone() 1550 1551 Level: beginner 1552 1553 .seealso: DMPlexCreate(), DMPlexSetChart(), DMPlexSetConeSize(), DMPlexSetCone() 1554 @*/ 1555 PetscErrorCode DMPlexSymmetrize(DM dm) 1556 { 1557 DM_Plex *mesh = (DM_Plex*) dm->data; 1558 PetscInt *offsets; 1559 PetscInt supportSize; 1560 PetscInt pStart, pEnd, p; 1561 PetscErrorCode ierr; 1562 1563 PetscFunctionBegin; 1564 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1565 if (mesh->supports) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Supports were already setup in this DMPlex"); 1566 /* Calculate support sizes */ 1567 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1568 for (p = pStart; p < pEnd; ++p) { 1569 PetscInt dof, off, c; 1570 1571 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 1572 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 1573 for (c = off; c < off+dof; ++c) { 1574 ierr = PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);CHKERRQ(ierr); 1575 } 1576 } 1577 for (p = pStart; p < pEnd; ++p) { 1578 PetscInt dof; 1579 1580 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 1581 1582 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof); 1583 } 1584 ierr = PetscSectionSetUp(mesh->supportSection);CHKERRQ(ierr); 1585 /* Calculate supports */ 1586 ierr = PetscSectionGetStorageSize(mesh->supportSection, &supportSize);CHKERRQ(ierr); 1587 ierr = PetscMalloc1(supportSize, &mesh->supports);CHKERRQ(ierr); 1588 ierr = PetscCalloc1(pEnd - pStart, &offsets);CHKERRQ(ierr); 1589 for (p = pStart; p < pEnd; ++p) { 1590 PetscInt dof, off, c; 1591 1592 ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr); 1593 ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr); 1594 for (c = off; c < off+dof; ++c) { 1595 const PetscInt q = mesh->cones[c]; 1596 PetscInt offS; 1597 1598 ierr = PetscSectionGetOffset(mesh->supportSection, q, &offS);CHKERRQ(ierr); 1599 1600 mesh->supports[offS+offsets[q]] = p; 1601 ++offsets[q]; 1602 } 1603 } 1604 ierr = PetscFree(offsets);CHKERRQ(ierr); 1605 PetscFunctionReturn(0); 1606 } 1607 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "DMPlexStratify" 1610 /*@ 1611 DMPlexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and 1612 can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the 1613 same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in 1614 the DAG. 1615 1616 Not collective 1617 1618 Input Parameter: 1619 . mesh - The DMPlex 1620 1621 Output Parameter: 1622 1623 Notes: 1624 Concretely, DMPlexStratify() creates a new label named "depth" containing the dimension of each element: 0 for vertices, 1625 1 for edges, and so on. The depth label can be accessed through DMPlexGetDepthLabel() or DMPlexGetDepthStratum(), or 1626 manually via DMPlexGetLabel(). The height is defined implicitly by height = maxDimension - depth, and can be accessed 1627 via DMPlexGetHeightStratum(). For example, cells have height 0 and faces have height 1. 1628 1629 DMPlexStratify() should be called after all calls to DMPlexSymmetrize() 1630 1631 Level: beginner 1632 1633 .seealso: DMPlexCreate(), DMPlexSymmetrize() 1634 @*/ 1635 PetscErrorCode DMPlexStratify(DM dm) 1636 { 1637 DMLabel label; 1638 PetscInt pStart, pEnd, p; 1639 PetscInt numRoots = 0, numLeaves = 0; 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1644 ierr = PetscLogEventBegin(DMPLEX_Stratify,dm,0,0,0);CHKERRQ(ierr); 1645 /* Calculate depth */ 1646 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1647 ierr = DMPlexCreateLabel(dm, "depth");CHKERRQ(ierr); 1648 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 1649 /* Initialize roots and count leaves */ 1650 for (p = pStart; p < pEnd; ++p) { 1651 PetscInt coneSize, supportSize; 1652 1653 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1654 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 1655 if (!coneSize && supportSize) { 1656 ++numRoots; 1657 ierr = DMLabelSetValue(label, p, 0);CHKERRQ(ierr); 1658 } else if (!supportSize && coneSize) { 1659 ++numLeaves; 1660 } else if (!supportSize && !coneSize) { 1661 /* Isolated points */ 1662 ierr = DMLabelSetValue(label, p, 0);CHKERRQ(ierr); 1663 } 1664 } 1665 if (numRoots + numLeaves == (pEnd - pStart)) { 1666 for (p = pStart; p < pEnd; ++p) { 1667 PetscInt coneSize, supportSize; 1668 1669 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1670 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 1671 if (!supportSize && coneSize) { 1672 ierr = DMLabelSetValue(label, p, 1);CHKERRQ(ierr); 1673 } 1674 } 1675 } else { 1676 IS pointIS; 1677 PetscInt numPoints = 0, level = 0; 1678 1679 ierr = DMLabelGetStratumIS(label, level, &pointIS);CHKERRQ(ierr); 1680 if (pointIS) {ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr);} 1681 while (numPoints) { 1682 const PetscInt *points; 1683 const PetscInt newLevel = level+1; 1684 1685 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1686 for (p = 0; p < numPoints; ++p) { 1687 const PetscInt point = points[p]; 1688 const PetscInt *support; 1689 PetscInt supportSize, s; 1690 1691 ierr = DMPlexGetSupportSize(dm, point, &supportSize);CHKERRQ(ierr); 1692 ierr = DMPlexGetSupport(dm, point, &support);CHKERRQ(ierr); 1693 for (s = 0; s < supportSize; ++s) { 1694 ierr = DMLabelSetValue(label, support[s], newLevel);CHKERRQ(ierr); 1695 } 1696 } 1697 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1698 ++level; 1699 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1700 ierr = DMLabelGetStratumIS(label, level, &pointIS);CHKERRQ(ierr); 1701 if (pointIS) {ierr = ISGetLocalSize(pointIS, &numPoints);CHKERRQ(ierr);} 1702 else {numPoints = 0;} 1703 } 1704 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1705 } 1706 ierr = PetscLogEventEnd(DMPLEX_Stratify,dm,0,0,0);CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "DMPlexGetJoin" 1712 /*@C 1713 DMPlexGetJoin - Get an array for the join of the set of points 1714 1715 Not Collective 1716 1717 Input Parameters: 1718 + dm - The DMPlex object 1719 . numPoints - The number of input points for the join 1720 - points - The input points 1721 1722 Output Parameters: 1723 + numCoveredPoints - The number of points in the join 1724 - coveredPoints - The points in the join 1725 1726 Level: intermediate 1727 1728 Note: Currently, this is restricted to a single level join 1729 1730 Fortran Notes: 1731 Since it returns an array, this routine is only available in Fortran 90, and you must 1732 include petsc.h90 in your code. 1733 1734 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1735 1736 .keywords: mesh 1737 .seealso: DMPlexRestoreJoin(), DMPlexGetMeet() 1738 @*/ 1739 PetscErrorCode DMPlexGetJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints) 1740 { 1741 DM_Plex *mesh = (DM_Plex*) dm->data; 1742 PetscInt *join[2]; 1743 PetscInt joinSize, i = 0; 1744 PetscInt dof, off, p, c, m; 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1749 PetscValidPointer(points, 2); 1750 PetscValidPointer(numCoveredPoints, 3); 1751 PetscValidPointer(coveredPoints, 4); 1752 ierr = DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[0]);CHKERRQ(ierr); 1753 ierr = DMGetWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1]);CHKERRQ(ierr); 1754 /* Copy in support of first point */ 1755 ierr = PetscSectionGetDof(mesh->supportSection, points[0], &dof);CHKERRQ(ierr); 1756 ierr = PetscSectionGetOffset(mesh->supportSection, points[0], &off);CHKERRQ(ierr); 1757 for (joinSize = 0; joinSize < dof; ++joinSize) { 1758 join[i][joinSize] = mesh->supports[off+joinSize]; 1759 } 1760 /* Check each successive support */ 1761 for (p = 1; p < numPoints; ++p) { 1762 PetscInt newJoinSize = 0; 1763 1764 ierr = PetscSectionGetDof(mesh->supportSection, points[p], &dof);CHKERRQ(ierr); 1765 ierr = PetscSectionGetOffset(mesh->supportSection, points[p], &off);CHKERRQ(ierr); 1766 for (c = 0; c < dof; ++c) { 1767 const PetscInt point = mesh->supports[off+c]; 1768 1769 for (m = 0; m < joinSize; ++m) { 1770 if (point == join[i][m]) { 1771 join[1-i][newJoinSize++] = point; 1772 break; 1773 } 1774 } 1775 } 1776 joinSize = newJoinSize; 1777 i = 1-i; 1778 } 1779 *numCoveredPoints = joinSize; 1780 *coveredPoints = join[i]; 1781 ierr = DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);CHKERRQ(ierr); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 #undef __FUNCT__ 1786 #define __FUNCT__ "DMPlexRestoreJoin" 1787 /*@C 1788 DMPlexRestoreJoin - Restore an array for the join of the set of points 1789 1790 Not Collective 1791 1792 Input Parameters: 1793 + dm - The DMPlex object 1794 . numPoints - The number of input points for the join 1795 - points - The input points 1796 1797 Output Parameters: 1798 + numCoveredPoints - The number of points in the join 1799 - coveredPoints - The points in the join 1800 1801 Fortran Notes: 1802 Since it returns an array, this routine is only available in Fortran 90, and you must 1803 include petsc.h90 in your code. 1804 1805 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1806 1807 Level: intermediate 1808 1809 .keywords: mesh 1810 .seealso: DMPlexGetJoin(), DMPlexGetFullJoin(), DMPlexGetMeet() 1811 @*/ 1812 PetscErrorCode DMPlexRestoreJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints) 1813 { 1814 PetscErrorCode ierr; 1815 1816 PetscFunctionBegin; 1817 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1818 if (points) PetscValidIntPointer(points,3); 1819 if (numCoveredPoints) PetscValidIntPointer(numCoveredPoints,4); 1820 PetscValidPointer(coveredPoints, 5); 1821 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);CHKERRQ(ierr); 1822 if (numCoveredPoints) *numCoveredPoints = 0; 1823 PetscFunctionReturn(0); 1824 } 1825 1826 #undef __FUNCT__ 1827 #define __FUNCT__ "DMPlexGetFullJoin" 1828 /*@C 1829 DMPlexGetFullJoin - Get an array for the join of the set of points 1830 1831 Not Collective 1832 1833 Input Parameters: 1834 + dm - The DMPlex object 1835 . numPoints - The number of input points for the join 1836 - points - The input points 1837 1838 Output Parameters: 1839 + numCoveredPoints - The number of points in the join 1840 - coveredPoints - The points in the join 1841 1842 Fortran Notes: 1843 Since it returns an array, this routine is only available in Fortran 90, and you must 1844 include petsc.h90 in your code. 1845 1846 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1847 1848 Level: intermediate 1849 1850 .keywords: mesh 1851 .seealso: DMPlexGetJoin(), DMPlexRestoreJoin(), DMPlexGetMeet() 1852 @*/ 1853 PetscErrorCode DMPlexGetFullJoin(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints) 1854 { 1855 DM_Plex *mesh = (DM_Plex*) dm->data; 1856 PetscInt *offsets, **closures; 1857 PetscInt *join[2]; 1858 PetscInt depth = 0, maxSize, joinSize = 0, i = 0; 1859 PetscInt p, d, c, m; 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1864 PetscValidPointer(points, 2); 1865 PetscValidPointer(numCoveredPoints, 3); 1866 PetscValidPointer(coveredPoints, 4); 1867 1868 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1869 ierr = PetscCalloc1(numPoints, &closures);CHKERRQ(ierr); 1870 ierr = DMGetWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);CHKERRQ(ierr); 1871 maxSize = PetscPowInt(mesh->maxSupportSize,depth+1); 1872 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &join[0]);CHKERRQ(ierr); 1873 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &join[1]);CHKERRQ(ierr); 1874 1875 for (p = 0; p < numPoints; ++p) { 1876 PetscInt closureSize; 1877 1878 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_FALSE, &closureSize, &closures[p]);CHKERRQ(ierr); 1879 1880 offsets[p*(depth+2)+0] = 0; 1881 for (d = 0; d < depth+1; ++d) { 1882 PetscInt pStart, pEnd, i; 1883 1884 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 1885 for (i = offsets[p*(depth+2)+d]; i < closureSize; ++i) { 1886 if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) { 1887 offsets[p*(depth+2)+d+1] = i; 1888 break; 1889 } 1890 } 1891 if (i == closureSize) offsets[p*(depth+2)+d+1] = i; 1892 } 1893 if (offsets[p*(depth+2)+depth+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(depth+2)+depth+1], closureSize); 1894 } 1895 for (d = 0; d < depth+1; ++d) { 1896 PetscInt dof; 1897 1898 /* Copy in support of first point */ 1899 dof = offsets[d+1] - offsets[d]; 1900 for (joinSize = 0; joinSize < dof; ++joinSize) { 1901 join[i][joinSize] = closures[0][(offsets[d]+joinSize)*2]; 1902 } 1903 /* Check each successive cone */ 1904 for (p = 1; p < numPoints && joinSize; ++p) { 1905 PetscInt newJoinSize = 0; 1906 1907 dof = offsets[p*(depth+2)+d+1] - offsets[p*(depth+2)+d]; 1908 for (c = 0; c < dof; ++c) { 1909 const PetscInt point = closures[p][(offsets[p*(depth+2)+d]+c)*2]; 1910 1911 for (m = 0; m < joinSize; ++m) { 1912 if (point == join[i][m]) { 1913 join[1-i][newJoinSize++] = point; 1914 break; 1915 } 1916 } 1917 } 1918 joinSize = newJoinSize; 1919 i = 1-i; 1920 } 1921 if (joinSize) break; 1922 } 1923 *numCoveredPoints = joinSize; 1924 *coveredPoints = join[i]; 1925 for (p = 0; p < numPoints; ++p) { 1926 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_FALSE, NULL, &closures[p]);CHKERRQ(ierr); 1927 } 1928 ierr = PetscFree(closures);CHKERRQ(ierr); 1929 ierr = DMRestoreWorkArray(dm, numPoints*(depth+2), PETSC_INT, &offsets);CHKERRQ(ierr); 1930 ierr = DMRestoreWorkArray(dm, mesh->maxSupportSize, PETSC_INT, &join[1-i]);CHKERRQ(ierr); 1931 PetscFunctionReturn(0); 1932 } 1933 1934 #undef __FUNCT__ 1935 #define __FUNCT__ "DMPlexGetMeet" 1936 /*@C 1937 DMPlexGetMeet - Get an array for the meet of the set of points 1938 1939 Not Collective 1940 1941 Input Parameters: 1942 + dm - The DMPlex object 1943 . numPoints - The number of input points for the meet 1944 - points - The input points 1945 1946 Output Parameters: 1947 + numCoveredPoints - The number of points in the meet 1948 - coveredPoints - The points in the meet 1949 1950 Level: intermediate 1951 1952 Note: Currently, this is restricted to a single level meet 1953 1954 Fortran Notes: 1955 Since it returns an array, this routine is only available in Fortran 90, and you must 1956 include petsc.h90 in your code. 1957 1958 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 1959 1960 .keywords: mesh 1961 .seealso: DMPlexRestoreMeet(), DMPlexGetJoin() 1962 @*/ 1963 PetscErrorCode DMPlexGetMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints) 1964 { 1965 DM_Plex *mesh = (DM_Plex*) dm->data; 1966 PetscInt *meet[2]; 1967 PetscInt meetSize, i = 0; 1968 PetscInt dof, off, p, c, m; 1969 PetscErrorCode ierr; 1970 1971 PetscFunctionBegin; 1972 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1973 PetscValidPointer(points, 2); 1974 PetscValidPointer(numCoveringPoints, 3); 1975 PetscValidPointer(coveringPoints, 4); 1976 ierr = DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[0]);CHKERRQ(ierr); 1977 ierr = DMGetWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1]);CHKERRQ(ierr); 1978 /* Copy in cone of first point */ 1979 ierr = PetscSectionGetDof(mesh->coneSection, points[0], &dof);CHKERRQ(ierr); 1980 ierr = PetscSectionGetOffset(mesh->coneSection, points[0], &off);CHKERRQ(ierr); 1981 for (meetSize = 0; meetSize < dof; ++meetSize) { 1982 meet[i][meetSize] = mesh->cones[off+meetSize]; 1983 } 1984 /* Check each successive cone */ 1985 for (p = 1; p < numPoints; ++p) { 1986 PetscInt newMeetSize = 0; 1987 1988 ierr = PetscSectionGetDof(mesh->coneSection, points[p], &dof);CHKERRQ(ierr); 1989 ierr = PetscSectionGetOffset(mesh->coneSection, points[p], &off);CHKERRQ(ierr); 1990 for (c = 0; c < dof; ++c) { 1991 const PetscInt point = mesh->cones[off+c]; 1992 1993 for (m = 0; m < meetSize; ++m) { 1994 if (point == meet[i][m]) { 1995 meet[1-i][newMeetSize++] = point; 1996 break; 1997 } 1998 } 1999 } 2000 meetSize = newMeetSize; 2001 i = 1-i; 2002 } 2003 *numCoveringPoints = meetSize; 2004 *coveringPoints = meet[i]; 2005 ierr = DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);CHKERRQ(ierr); 2006 PetscFunctionReturn(0); 2007 } 2008 2009 #undef __FUNCT__ 2010 #define __FUNCT__ "DMPlexRestoreMeet" 2011 /*@C 2012 DMPlexRestoreMeet - Restore an array for the meet of the set of points 2013 2014 Not Collective 2015 2016 Input Parameters: 2017 + dm - The DMPlex object 2018 . numPoints - The number of input points for the meet 2019 - points - The input points 2020 2021 Output Parameters: 2022 + numCoveredPoints - The number of points in the meet 2023 - coveredPoints - The points in the meet 2024 2025 Level: intermediate 2026 2027 Fortran Notes: 2028 Since it returns an array, this routine is only available in Fortran 90, and you must 2029 include petsc.h90 in your code. 2030 2031 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 2032 2033 .keywords: mesh 2034 .seealso: DMPlexGetMeet(), DMPlexGetFullMeet(), DMPlexGetJoin() 2035 @*/ 2036 PetscErrorCode DMPlexRestoreMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints) 2037 { 2038 PetscErrorCode ierr; 2039 2040 PetscFunctionBegin; 2041 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2042 if (points) PetscValidIntPointer(points,3); 2043 if (numCoveredPoints) PetscValidIntPointer(numCoveredPoints,4); 2044 PetscValidPointer(coveredPoints,5); 2045 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, (void*) coveredPoints);CHKERRQ(ierr); 2046 if (numCoveredPoints) *numCoveredPoints = 0; 2047 PetscFunctionReturn(0); 2048 } 2049 2050 #undef __FUNCT__ 2051 #define __FUNCT__ "DMPlexGetFullMeet" 2052 /*@C 2053 DMPlexGetFullMeet - Get an array for the meet of the set of points 2054 2055 Not Collective 2056 2057 Input Parameters: 2058 + dm - The DMPlex object 2059 . numPoints - The number of input points for the meet 2060 - points - The input points 2061 2062 Output Parameters: 2063 + numCoveredPoints - The number of points in the meet 2064 - coveredPoints - The points in the meet 2065 2066 Level: intermediate 2067 2068 Fortran Notes: 2069 Since it returns an array, this routine is only available in Fortran 90, and you must 2070 include petsc.h90 in your code. 2071 2072 The numCoveredPoints argument is not present in the Fortran 90 binding since it is internal to the array. 2073 2074 .keywords: mesh 2075 .seealso: DMPlexGetMeet(), DMPlexRestoreMeet(), DMPlexGetJoin() 2076 @*/ 2077 PetscErrorCode DMPlexGetFullMeet(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints) 2078 { 2079 DM_Plex *mesh = (DM_Plex*) dm->data; 2080 PetscInt *offsets, **closures; 2081 PetscInt *meet[2]; 2082 PetscInt height = 0, maxSize, meetSize = 0, i = 0; 2083 PetscInt p, h, c, m; 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2088 PetscValidPointer(points, 2); 2089 PetscValidPointer(numCoveredPoints, 3); 2090 PetscValidPointer(coveredPoints, 4); 2091 2092 ierr = DMPlexGetDepth(dm, &height);CHKERRQ(ierr); 2093 ierr = PetscMalloc1(numPoints, &closures);CHKERRQ(ierr); 2094 ierr = DMGetWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);CHKERRQ(ierr); 2095 maxSize = PetscPowInt(mesh->maxConeSize,height+1); 2096 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[0]);CHKERRQ(ierr); 2097 ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &meet[1]);CHKERRQ(ierr); 2098 2099 for (p = 0; p < numPoints; ++p) { 2100 PetscInt closureSize; 2101 2102 ierr = DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closures[p]);CHKERRQ(ierr); 2103 2104 offsets[p*(height+2)+0] = 0; 2105 for (h = 0; h < height+1; ++h) { 2106 PetscInt pStart, pEnd, i; 2107 2108 ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr); 2109 for (i = offsets[p*(height+2)+h]; i < closureSize; ++i) { 2110 if ((pStart > closures[p][i*2]) || (pEnd <= closures[p][i*2])) { 2111 offsets[p*(height+2)+h+1] = i; 2112 break; 2113 } 2114 } 2115 if (i == closureSize) offsets[p*(height+2)+h+1] = i; 2116 } 2117 if (offsets[p*(height+2)+height+1] != closureSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Total size of closure %D should be %D", offsets[p*(height+2)+height+1], closureSize); 2118 } 2119 for (h = 0; h < height+1; ++h) { 2120 PetscInt dof; 2121 2122 /* Copy in cone of first point */ 2123 dof = offsets[h+1] - offsets[h]; 2124 for (meetSize = 0; meetSize < dof; ++meetSize) { 2125 meet[i][meetSize] = closures[0][(offsets[h]+meetSize)*2]; 2126 } 2127 /* Check each successive cone */ 2128 for (p = 1; p < numPoints && meetSize; ++p) { 2129 PetscInt newMeetSize = 0; 2130 2131 dof = offsets[p*(height+2)+h+1] - offsets[p*(height+2)+h]; 2132 for (c = 0; c < dof; ++c) { 2133 const PetscInt point = closures[p][(offsets[p*(height+2)+h]+c)*2]; 2134 2135 for (m = 0; m < meetSize; ++m) { 2136 if (point == meet[i][m]) { 2137 meet[1-i][newMeetSize++] = point; 2138 break; 2139 } 2140 } 2141 } 2142 meetSize = newMeetSize; 2143 i = 1-i; 2144 } 2145 if (meetSize) break; 2146 } 2147 *numCoveredPoints = meetSize; 2148 *coveredPoints = meet[i]; 2149 for (p = 0; p < numPoints; ++p) { 2150 ierr = DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, NULL, &closures[p]);CHKERRQ(ierr); 2151 } 2152 ierr = PetscFree(closures);CHKERRQ(ierr); 2153 ierr = DMRestoreWorkArray(dm, numPoints*(height+2), PETSC_INT, &offsets);CHKERRQ(ierr); 2154 ierr = DMRestoreWorkArray(dm, mesh->maxConeSize, PETSC_INT, &meet[1-i]);CHKERRQ(ierr); 2155 PetscFunctionReturn(0); 2156 } 2157 2158 #undef __FUNCT__ 2159 #define __FUNCT__ "DMPlexEqual" 2160 /*@C 2161 DMPlexEqual - Determine if two DMs have the same topology 2162 2163 Not Collective 2164 2165 Input Parameters: 2166 + dmA - A DMPlex object 2167 - dmB - A DMPlex object 2168 2169 Output Parameters: 2170 . equal - PETSC_TRUE if the topologies are identical 2171 2172 Level: intermediate 2173 2174 Notes: 2175 We are not solving graph isomorphism, so we do not permutation. 2176 2177 .keywords: mesh 2178 .seealso: DMPlexGetCone() 2179 @*/ 2180 PetscErrorCode DMPlexEqual(DM dmA, DM dmB, PetscBool *equal) 2181 { 2182 PetscInt depth, depthB, pStart, pEnd, pStartB, pEndB, p; 2183 PetscErrorCode ierr; 2184 2185 PetscFunctionBegin; 2186 PetscValidHeaderSpecific(dmA, DM_CLASSID, 1); 2187 PetscValidHeaderSpecific(dmB, DM_CLASSID, 2); 2188 PetscValidPointer(equal, 3); 2189 2190 *equal = PETSC_FALSE; 2191 ierr = DMPlexGetDepth(dmA, &depth);CHKERRQ(ierr); 2192 ierr = DMPlexGetDepth(dmB, &depthB);CHKERRQ(ierr); 2193 if (depth != depthB) PetscFunctionReturn(0); 2194 ierr = DMPlexGetChart(dmA, &pStart, &pEnd);CHKERRQ(ierr); 2195 ierr = DMPlexGetChart(dmB, &pStartB, &pEndB);CHKERRQ(ierr); 2196 if ((pStart != pStartB) || (pEnd != pEndB)) PetscFunctionReturn(0); 2197 for (p = pStart; p < pEnd; ++p) { 2198 const PetscInt *cone, *coneB, *ornt, *orntB, *support, *supportB; 2199 PetscInt coneSize, coneSizeB, c, supportSize, supportSizeB, s; 2200 2201 ierr = DMPlexGetConeSize(dmA, p, &coneSize);CHKERRQ(ierr); 2202 ierr = DMPlexGetCone(dmA, p, &cone);CHKERRQ(ierr); 2203 ierr = DMPlexGetConeOrientation(dmA, p, &ornt);CHKERRQ(ierr); 2204 ierr = DMPlexGetConeSize(dmB, p, &coneSizeB);CHKERRQ(ierr); 2205 ierr = DMPlexGetCone(dmB, p, &coneB);CHKERRQ(ierr); 2206 ierr = DMPlexGetConeOrientation(dmB, p, &orntB);CHKERRQ(ierr); 2207 if (coneSize != coneSizeB) PetscFunctionReturn(0); 2208 for (c = 0; c < coneSize; ++c) { 2209 if (cone[c] != coneB[c]) PetscFunctionReturn(0); 2210 if (ornt[c] != orntB[c]) PetscFunctionReturn(0); 2211 } 2212 ierr = DMPlexGetSupportSize(dmA, p, &supportSize);CHKERRQ(ierr); 2213 ierr = DMPlexGetSupport(dmA, p, &support);CHKERRQ(ierr); 2214 ierr = DMPlexGetSupportSize(dmB, p, &supportSizeB);CHKERRQ(ierr); 2215 ierr = DMPlexGetSupport(dmB, p, &supportB);CHKERRQ(ierr); 2216 if (supportSize != supportSizeB) PetscFunctionReturn(0); 2217 for (s = 0; s < supportSize; ++s) { 2218 if (support[s] != supportB[s]) PetscFunctionReturn(0); 2219 } 2220 } 2221 *equal = PETSC_TRUE; 2222 PetscFunctionReturn(0); 2223 } 2224 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "DMPlexGetNumFaceVertices" 2227 PetscErrorCode DMPlexGetNumFaceVertices(DM dm, PetscInt cellDim, PetscInt numCorners, PetscInt *numFaceVertices) 2228 { 2229 MPI_Comm comm; 2230 PetscErrorCode ierr; 2231 2232 PetscFunctionBegin; 2233 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 2234 PetscValidPointer(numFaceVertices,3); 2235 switch (cellDim) { 2236 case 0: 2237 *numFaceVertices = 0; 2238 break; 2239 case 1: 2240 *numFaceVertices = 1; 2241 break; 2242 case 2: 2243 switch (numCorners) { 2244 case 3: /* triangle */ 2245 *numFaceVertices = 2; /* Edge has 2 vertices */ 2246 break; 2247 case 4: /* quadrilateral */ 2248 *numFaceVertices = 2; /* Edge has 2 vertices */ 2249 break; 2250 case 6: /* quadratic triangle, tri and quad cohesive Lagrange cells */ 2251 *numFaceVertices = 3; /* Edge has 3 vertices */ 2252 break; 2253 case 9: /* quadratic quadrilateral, quadratic quad cohesive Lagrange cells */ 2254 *numFaceVertices = 3; /* Edge has 3 vertices */ 2255 break; 2256 default: 2257 SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim); 2258 } 2259 break; 2260 case 3: 2261 switch (numCorners) { 2262 case 4: /* tetradehdron */ 2263 *numFaceVertices = 3; /* Face has 3 vertices */ 2264 break; 2265 case 6: /* tet cohesive cells */ 2266 *numFaceVertices = 4; /* Face has 4 vertices */ 2267 break; 2268 case 8: /* hexahedron */ 2269 *numFaceVertices = 4; /* Face has 4 vertices */ 2270 break; 2271 case 9: /* tet cohesive Lagrange cells */ 2272 *numFaceVertices = 6; /* Face has 6 vertices */ 2273 break; 2274 case 10: /* quadratic tetrahedron */ 2275 *numFaceVertices = 6; /* Face has 6 vertices */ 2276 break; 2277 case 12: /* hex cohesive Lagrange cells */ 2278 *numFaceVertices = 6; /* Face has 6 vertices */ 2279 break; 2280 case 18: /* quadratic tet cohesive Lagrange cells */ 2281 *numFaceVertices = 6; /* Face has 6 vertices */ 2282 break; 2283 case 27: /* quadratic hexahedron, quadratic hex cohesive Lagrange cells */ 2284 *numFaceVertices = 9; /* Face has 9 vertices */ 2285 break; 2286 default: 2287 SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of face corners %d for dimension %d", numCorners, cellDim); 2288 } 2289 break; 2290 default: 2291 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid cell dimension %d", cellDim); 2292 } 2293 PetscFunctionReturn(0); 2294 } 2295 2296 #undef __FUNCT__ 2297 #define __FUNCT__ "DMPlexOrient" 2298 /* Trys to give the mesh a consistent orientation */ 2299 PetscErrorCode DMPlexOrient(DM dm) 2300 { 2301 PetscBT seenCells, flippedCells, seenFaces; 2302 PetscInt *faceFIFO, fTop, fBottom; 2303 PetscInt dim, h, cStart, cEnd, c, fStart, fEnd, face, maxConeSize, *revcone, *revconeO; 2304 PetscErrorCode ierr; 2305 2306 PetscFunctionBegin; 2307 /* Truth Table 2308 mismatch flips do action mismatch flipA ^ flipB action 2309 F 0 flips no F F F 2310 F 1 flip yes F T T 2311 F 2 flips no T F T 2312 T 0 flips yes T T F 2313 T 1 flip no 2314 T 2 flips yes 2315 */ 2316 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 2317 ierr = DMPlexGetVTKCellHeight(dm, &h);CHKERRQ(ierr); 2318 ierr = DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);CHKERRQ(ierr); 2319 ierr = DMPlexGetHeightStratum(dm, h+1, &fStart, &fEnd);CHKERRQ(ierr); 2320 ierr = PetscBTCreate(cEnd - cStart, &seenCells);CHKERRQ(ierr); 2321 ierr = PetscBTMemzero(cEnd - cStart, seenCells);CHKERRQ(ierr); 2322 ierr = PetscBTCreate(cEnd - cStart, &flippedCells);CHKERRQ(ierr); 2323 ierr = PetscBTMemzero(cEnd - cStart, flippedCells);CHKERRQ(ierr); 2324 ierr = PetscBTCreate(fEnd - fStart, &seenFaces);CHKERRQ(ierr); 2325 ierr = PetscBTMemzero(fEnd - fStart, seenFaces);CHKERRQ(ierr); 2326 ierr = PetscMalloc1((fEnd - fStart), &faceFIFO);CHKERRQ(ierr); 2327 fTop = fBottom = 0; 2328 /* Initialize FIFO with first cell */ 2329 if (cEnd > cStart) { 2330 const PetscInt *cone; 2331 PetscInt coneSize; 2332 2333 ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr); 2334 ierr = DMPlexGetCone(dm, cStart, &cone);CHKERRQ(ierr); 2335 for (c = 0; c < coneSize; ++c) { 2336 faceFIFO[fBottom++] = cone[c]; 2337 ierr = PetscBTSet(seenFaces, cone[c]-fStart);CHKERRQ(ierr); 2338 } 2339 } 2340 /* Consider each face in FIFO */ 2341 while (fTop < fBottom) { 2342 const PetscInt *support, *coneA, *coneB, *coneOA, *coneOB; 2343 PetscInt supportSize, coneSizeA, coneSizeB, posA = -1, posB = -1; 2344 PetscInt seenA, flippedA, seenB, flippedB, mismatch; 2345 2346 face = faceFIFO[fTop++]; 2347 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2348 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2349 if (supportSize < 2) continue; 2350 if (supportSize != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Faces should separate only two cells, not %d", supportSize); 2351 seenA = PetscBTLookup(seenCells, support[0]-cStart); 2352 flippedA = PetscBTLookup(flippedCells, support[0]-cStart) ? 1 : 0; 2353 seenB = PetscBTLookup(seenCells, support[1]-cStart); 2354 flippedB = PetscBTLookup(flippedCells, support[1]-cStart) ? 1 : 0; 2355 2356 ierr = DMPlexGetConeSize(dm, support[0], &coneSizeA);CHKERRQ(ierr); 2357 ierr = DMPlexGetConeSize(dm, support[1], &coneSizeB);CHKERRQ(ierr); 2358 ierr = DMPlexGetCone(dm, support[0], &coneA);CHKERRQ(ierr); 2359 ierr = DMPlexGetCone(dm, support[1], &coneB);CHKERRQ(ierr); 2360 ierr = DMPlexGetConeOrientation(dm, support[0], &coneOA);CHKERRQ(ierr); 2361 ierr = DMPlexGetConeOrientation(dm, support[1], &coneOB);CHKERRQ(ierr); 2362 for (c = 0; c < coneSizeA; ++c) { 2363 if (!PetscBTLookup(seenFaces, coneA[c]-fStart)) { 2364 faceFIFO[fBottom++] = coneA[c]; 2365 ierr = PetscBTSet(seenFaces, coneA[c]-fStart);CHKERRQ(ierr); 2366 } 2367 if (coneA[c] == face) posA = c; 2368 if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart); 2369 } 2370 if (posA < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[0]); 2371 for (c = 0; c < coneSizeB; ++c) { 2372 if (!PetscBTLookup(seenFaces, coneB[c]-fStart)) { 2373 faceFIFO[fBottom++] = coneB[c]; 2374 ierr = PetscBTSet(seenFaces, coneB[c]-fStart);CHKERRQ(ierr); 2375 } 2376 if (coneB[c] == face) posB = c; 2377 if (fBottom > fEnd-fStart) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face %d was pushed exceeding capacity %d > %d", coneA[c], fBottom, fEnd-fStart); 2378 } 2379 if (posB < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d could not be located in cell %d", face, support[1]); 2380 2381 if (dim == 1) { 2382 mismatch = posA == posB; 2383 } else { 2384 mismatch = coneOA[posA] == coneOB[posB]; 2385 } 2386 2387 if (mismatch ^ (flippedA ^ flippedB)) { 2388 if (seenA && seenB) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen cells %d and %d do not match: Fault mesh is non-orientable", support[0], support[1]); 2389 if (!seenA && !flippedA) { 2390 ierr = PetscBTSet(flippedCells, support[0]-cStart);CHKERRQ(ierr); 2391 } else if (!seenB && !flippedB) { 2392 ierr = PetscBTSet(flippedCells, support[1]-cStart);CHKERRQ(ierr); 2393 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 2394 } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 2395 ierr = PetscBTSet(seenCells, support[0]-cStart);CHKERRQ(ierr); 2396 ierr = PetscBTSet(seenCells, support[1]-cStart);CHKERRQ(ierr); 2397 } 2398 /* Now all subdomains are oriented, but we need a consistent parallel orientation */ 2399 { 2400 /* Find a representative face (edge) separating pairs of procs */ 2401 PetscSF sf; 2402 const PetscInt *lpoints; 2403 const PetscSFNode *rpoints; 2404 PetscInt *neighbors, *nranks; 2405 PetscInt numLeaves, numRoots, numNeighbors = 0, l, n; 2406 2407 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 2408 ierr = PetscSFGetGraph(sf, &numRoots, &numLeaves, &lpoints, &rpoints);CHKERRQ(ierr); 2409 if (numLeaves >= 0) { 2410 const PetscInt *cone, *ornt, *support; 2411 PetscInt coneSize, supportSize; 2412 int *rornt, *lornt; /* PetscSF cannot handle smaller than int */ 2413 PetscBool *match, flipped = PETSC_FALSE; 2414 2415 ierr = PetscMalloc1(numLeaves,&neighbors);CHKERRQ(ierr); 2416 /* I know this is p^2 time in general, but for bounded degree its alright */ 2417 for (l = 0; l < numLeaves; ++l) { 2418 const PetscInt face = lpoints[l]; 2419 if ((face >= fStart) && (face < fEnd)) { 2420 const PetscInt rank = rpoints[l].rank; 2421 for (n = 0; n < numNeighbors; ++n) if (rank == rpoints[neighbors[n]].rank) break; 2422 if (n >= numNeighbors) { 2423 PetscInt supportSize; 2424 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2425 if (supportSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Boundary faces should see one cell, not %d", supportSize); 2426 neighbors[numNeighbors++] = l; 2427 } 2428 } 2429 } 2430 ierr = PetscCalloc4(numNeighbors,&match,numNeighbors,&nranks,numRoots,&rornt,numRoots,&lornt);CHKERRQ(ierr); 2431 for (face = fStart; face < fEnd; ++face) { 2432 ierr = DMPlexGetSupportSize(dm, face, &supportSize);CHKERRQ(ierr); 2433 if (supportSize != 1) continue; 2434 ierr = DMPlexGetSupport(dm, face, &support);CHKERRQ(ierr); 2435 2436 ierr = DMPlexGetCone(dm, support[0], &cone);CHKERRQ(ierr); 2437 ierr = DMPlexGetConeSize(dm, support[0], &coneSize);CHKERRQ(ierr); 2438 ierr = DMPlexGetConeOrientation(dm, support[0], &ornt);CHKERRQ(ierr); 2439 for (c = 0; c < coneSize; ++c) if (cone[c] == face) break; 2440 if (dim == 1) { 2441 /* Use cone position instead, shifted to -1 or 1 */ 2442 rornt[face] = c*2-1; 2443 } else { 2444 if (PetscBTLookup(flippedCells, support[0]-cStart)) rornt[face] = ornt[c] < 0 ? -1 : 1; 2445 else rornt[face] = ornt[c] < 0 ? 1 : -1; 2446 } 2447 } 2448 /* Mark each edge with match or nomatch */ 2449 ierr = PetscSFBcastBegin(sf, MPI_INT, rornt, lornt);CHKERRQ(ierr); 2450 ierr = PetscSFBcastEnd(sf, MPI_INT, rornt, lornt);CHKERRQ(ierr); 2451 for (n = 0; n < numNeighbors; ++n) { 2452 const PetscInt face = lpoints[neighbors[n]]; 2453 2454 if (rornt[face]*lornt[face] < 0) match[n] = PETSC_TRUE; 2455 else match[n] = PETSC_FALSE; 2456 nranks[n] = rpoints[neighbors[n]].rank; 2457 } 2458 /* Collect the graph on 0 */ 2459 { 2460 MPI_Comm comm = PetscObjectComm((PetscObject) sf); 2461 Mat G; 2462 PetscBT seenProcs, flippedProcs; 2463 PetscInt *procFIFO, pTop, pBottom; 2464 PetscInt *adj = NULL; 2465 PetscBool *val = NULL; 2466 PetscMPIInt *recvcounts = NULL, *displs = NULL, p; 2467 PetscMPIInt N = numNeighbors, numProcs = 0, rank; 2468 PetscInt debug = 0; 2469 2470 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2471 if (!rank) {ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);} 2472 ierr = PetscCalloc2(numProcs,&recvcounts,numProcs+1,&displs);CHKERRQ(ierr); 2473 ierr = MPI_Gather(&N, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 2474 for (p = 0; p < numProcs; ++p) { 2475 displs[p+1] = displs[p] + recvcounts[p]; 2476 } 2477 if (!rank) {ierr = PetscMalloc2(displs[numProcs],&adj,displs[numProcs],&val);CHKERRQ(ierr);} 2478 ierr = MPI_Gatherv(nranks, numNeighbors, MPIU_INT, adj, recvcounts, displs, MPIU_INT, 0, comm);CHKERRQ(ierr); 2479 ierr = MPI_Gatherv(match, numNeighbors, MPIU_BOOL, val, recvcounts, displs, MPIU_BOOL, 0, comm);CHKERRQ(ierr); 2480 if (debug) { 2481 for (p = 0; p < numProcs; ++p) { 2482 ierr = PetscPrintf(comm, "Proc %d:\n", p); 2483 for (n = 0; n < recvcounts[p]; ++n) { 2484 ierr = PetscPrintf(comm, " edge %d (%d):\n", adj[displs[p]+n], val[displs[p]+n]); 2485 } 2486 } 2487 } 2488 /* Symmetrize the graph */ 2489 ierr = MatCreate(PETSC_COMM_SELF, &G);CHKERRQ(ierr); 2490 ierr = MatSetSizes(G, numProcs, numProcs, numProcs, numProcs);CHKERRQ(ierr); 2491 ierr = MatSetUp(G);CHKERRQ(ierr); 2492 for (p = 0; p < numProcs; ++p) { 2493 for (n = 0; n < recvcounts[p]; ++n) { 2494 const PetscInt r = p; 2495 const PetscInt q = adj[displs[p]+n]; 2496 const PetscScalar o = val[displs[p]+n] ? 1.0 : 0.0; 2497 2498 ierr = MatSetValues(G, 1, &r, 1, &q, &o, INSERT_VALUES);CHKERRQ(ierr); 2499 ierr = MatSetValues(G, 1, &q, 1, &r, &o, INSERT_VALUES);CHKERRQ(ierr); 2500 } 2501 } 2502 ierr = MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2503 ierr = MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2504 2505 ierr = PetscBTCreate(numProcs, &seenProcs);CHKERRQ(ierr); 2506 ierr = PetscBTMemzero(numProcs, seenProcs);CHKERRQ(ierr); 2507 ierr = PetscBTCreate(numProcs, &flippedProcs);CHKERRQ(ierr); 2508 ierr = PetscBTMemzero(numProcs, flippedProcs);CHKERRQ(ierr); 2509 ierr = PetscMalloc1(numProcs,&procFIFO);CHKERRQ(ierr); 2510 pTop = pBottom = 0; 2511 for (p = 0; p < numProcs; ++p) { 2512 if (PetscBTLookup(seenProcs, p)) continue; 2513 /* Initialize FIFO with next proc */ 2514 procFIFO[pBottom++] = p; 2515 ierr = PetscBTSet(seenProcs, p);CHKERRQ(ierr); 2516 /* Consider each proc in FIFO */ 2517 while (pTop < pBottom) { 2518 const PetscScalar *ornt; 2519 const PetscInt *neighbors; 2520 PetscInt proc, nproc, seen, flippedA, flippedB, mismatch, numNeighbors; 2521 2522 proc = procFIFO[pTop++]; 2523 flippedA = PetscBTLookup(flippedProcs, proc) ? 1 : 0; 2524 ierr = MatGetRow(G, proc, &numNeighbors, &neighbors, &ornt);CHKERRQ(ierr); 2525 /* Loop over neighboring procs */ 2526 for (n = 0; n < numNeighbors; ++n) { 2527 nproc = neighbors[n]; 2528 mismatch = PetscRealPart(ornt[n]) > 0.5 ? 0 : 1; 2529 seen = PetscBTLookup(seenProcs, nproc); 2530 flippedB = PetscBTLookup(flippedProcs, nproc) ? 1 : 0; 2531 2532 if (mismatch ^ (flippedA ^ flippedB)) { 2533 if (seen) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Previously seen procs %d and %d do not match: Fault mesh is non-orientable", proc, nproc); 2534 if (!flippedB) { 2535 ierr = PetscBTSet(flippedProcs, nproc);CHKERRQ(ierr); 2536 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Inconsistent mesh orientation: Fault mesh is non-orientable"); 2537 } else if (mismatch && flippedA && flippedB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Attempt to flip already flipped cell: Fault mesh is non-orientable"); 2538 if (!seen) { 2539 procFIFO[pBottom++] = nproc; 2540 ierr = PetscBTSet(seenProcs, nproc);CHKERRQ(ierr); 2541 } 2542 } 2543 } 2544 } 2545 ierr = PetscFree(procFIFO);CHKERRQ(ierr); 2546 ierr = MatDestroy(&G);CHKERRQ(ierr); 2547 2548 ierr = PetscFree2(recvcounts,displs);CHKERRQ(ierr); 2549 ierr = PetscFree2(adj,val);CHKERRQ(ierr); 2550 { 2551 PetscBool *flips; 2552 2553 ierr = PetscMalloc1(numProcs,&flips);CHKERRQ(ierr); 2554 for (p = 0; p < numProcs; ++p) { 2555 flips[p] = PetscBTLookup(flippedProcs, p) ? PETSC_TRUE : PETSC_FALSE; 2556 if (debug && flips[p]) {ierr = PetscPrintf(comm, "Flipping Proc %d:\n", p);} 2557 } 2558 ierr = MPI_Scatter(flips, 1, MPIU_BOOL, &flipped, 1, MPIU_BOOL, 0, comm);CHKERRQ(ierr); 2559 ierr = PetscFree(flips);CHKERRQ(ierr); 2560 } 2561 ierr = PetscBTDestroy(&seenProcs);CHKERRQ(ierr); 2562 ierr = PetscBTDestroy(&flippedProcs);CHKERRQ(ierr); 2563 } 2564 ierr = PetscFree4(match,nranks,rornt,lornt);CHKERRQ(ierr); 2565 ierr = PetscFree(neighbors);CHKERRQ(ierr); 2566 if (flipped) {for (c = cStart; c < cEnd; ++c) {ierr = PetscBTNegate(flippedCells, c-cStart);CHKERRQ(ierr);}} 2567 } 2568 } 2569 /* Reverse flipped cells in the mesh */ 2570 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRQ(ierr); 2571 ierr = DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revcone);CHKERRQ(ierr); 2572 ierr = DMGetWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);CHKERRQ(ierr); 2573 for (c = cStart; c < cEnd; ++c) { 2574 const PetscInt *cone, *coneO, *support; 2575 PetscInt coneSize, supportSize, faceSize, cp, sp; 2576 2577 if (!PetscBTLookup(flippedCells, c-cStart)) continue; 2578 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 2579 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 2580 ierr = DMPlexGetConeOrientation(dm, c, &coneO);CHKERRQ(ierr); 2581 for (cp = 0; cp < coneSize; ++cp) { 2582 const PetscInt rcp = coneSize-cp-1; 2583 2584 ierr = DMPlexGetConeSize(dm, cone[rcp], &faceSize);CHKERRQ(ierr); 2585 revcone[cp] = cone[rcp]; 2586 revconeO[cp] = coneO[rcp] >= 0 ? -(faceSize-coneO[rcp]) : faceSize+coneO[rcp]; 2587 } 2588 ierr = DMPlexSetCone(dm, c, revcone);CHKERRQ(ierr); 2589 ierr = DMPlexSetConeOrientation(dm, c, revconeO);CHKERRQ(ierr); 2590 /* Reverse orientations of support */ 2591 faceSize = coneSize; 2592 ierr = DMPlexGetSupportSize(dm, c, &supportSize);CHKERRQ(ierr); 2593 ierr = DMPlexGetSupport(dm, c, &support);CHKERRQ(ierr); 2594 for (sp = 0; sp < supportSize; ++sp) { 2595 ierr = DMPlexGetConeSize(dm, support[sp], &coneSize);CHKERRQ(ierr); 2596 ierr = DMPlexGetCone(dm, support[sp], &cone);CHKERRQ(ierr); 2597 ierr = DMPlexGetConeOrientation(dm, support[sp], &coneO);CHKERRQ(ierr); 2598 for (cp = 0; cp < coneSize; ++cp) { 2599 if (cone[cp] != c) continue; 2600 ierr = DMPlexInsertConeOrientation(dm, support[sp], cp, coneO[cp] >= 0 ? -(faceSize-coneO[cp]) : faceSize+coneO[cp]);CHKERRQ(ierr); 2601 } 2602 } 2603 } 2604 ierr = DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revcone);CHKERRQ(ierr); 2605 ierr = DMRestoreWorkArray(dm, maxConeSize, PETSC_INT, &revconeO);CHKERRQ(ierr); 2606 ierr = PetscBTDestroy(&seenCells);CHKERRQ(ierr); 2607 ierr = PetscBTDestroy(&flippedCells);CHKERRQ(ierr); 2608 ierr = PetscBTDestroy(&seenFaces);CHKERRQ(ierr); 2609 ierr = PetscFree(faceFIFO);CHKERRQ(ierr); 2610 PetscFunctionReturn(0); 2611 } 2612 2613 #undef __FUNCT__ 2614 #define __FUNCT__ "DMPlexLocalizeCoordinate_Internal" 2615 PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) 2616 { 2617 PetscInt d; 2618 2619 PetscFunctionBegin; 2620 if (!dm->maxCell) { 2621 for (d = 0; d < dim; ++d) out[d] = in[d]; 2622 } else { 2623 for (d = 0; d < dim; ++d) { 2624 if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) { 2625 out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 2626 } else { 2627 out[d] = in[d]; 2628 } 2629 } 2630 } 2631 PetscFunctionReturn(0); 2632 } 2633 2634 #undef __FUNCT__ 2635 #define __FUNCT__ "DMPlexLocalizeAddCoordinate_Internal" 2636 PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[]) 2637 { 2638 PetscInt d; 2639 2640 PetscFunctionBegin; 2641 if (!dm->maxCell) { 2642 for (d = 0; d < dim; ++d) out[d] += in[d]; 2643 } else { 2644 for (d = 0; d < dim; ++d) { 2645 if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) { 2646 out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d]; 2647 } else { 2648 out[d] += in[d]; 2649 } 2650 } 2651 } 2652 PetscFunctionReturn(0); 2653 } 2654 2655 #undef __FUNCT__ 2656 #define __FUNCT__ "DMPlexLocalizeCoordinates" 2657 PetscErrorCode DMPlexLocalizeCoordinates(DM dm) 2658 { 2659 PetscSection coordSection, cSection; 2660 Vec coordinates, cVec; 2661 PetscScalar *coords, *coords2, *anchor; 2662 PetscInt Nc, cStart, cEnd, c, vStart, vEnd, v, dof, d, off, off2, bs, coordSize; 2663 PetscErrorCode ierr; 2664 2665 PetscFunctionBegin; 2666 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2667 if (!dm->maxCell) PetscFunctionReturn(0); 2668 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 2669 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 2670 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 2671 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 2672 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);CHKERRQ(ierr); 2673 ierr = PetscSectionSetNumFields(cSection, 1);CHKERRQ(ierr); 2674 ierr = PetscSectionGetFieldComponents(coordSection, 0, &Nc);CHKERRQ(ierr); 2675 ierr = PetscSectionSetFieldComponents(cSection, 0, Nc);CHKERRQ(ierr); 2676 ierr = PetscSectionSetChart(cSection, cStart, vEnd);CHKERRQ(ierr); 2677 for (v = vStart; v < vEnd; ++v) { 2678 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 2679 ierr = PetscSectionSetDof(cSection, v, dof);CHKERRQ(ierr); 2680 ierr = PetscSectionSetFieldDof(cSection, v, 0, dof);CHKERRQ(ierr); 2681 } 2682 for (c = cStart; c < cEnd; ++c) { 2683 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, NULL);CHKERRQ(ierr); 2684 ierr = PetscSectionSetDof(cSection, c, dof);CHKERRQ(ierr); 2685 ierr = PetscSectionSetFieldDof(cSection, c, 0, dof);CHKERRQ(ierr); 2686 } 2687 ierr = PetscSectionSetUp(cSection);CHKERRQ(ierr); 2688 ierr = PetscSectionGetStorageSize(cSection, &coordSize);CHKERRQ(ierr); 2689 ierr = VecCreate(PetscObjectComm((PetscObject) dm), &cVec);CHKERRQ(ierr); 2690 ierr = VecGetBlockSize(coordinates, &bs);CHKERRQ(ierr); 2691 ierr = VecSetBlockSize(cVec, bs);CHKERRQ(ierr); 2692 ierr = VecSetSizes(cVec, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 2693 ierr = VecSetType(cVec,VECSTANDARD);CHKERRQ(ierr); 2694 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 2695 ierr = VecGetArray(cVec, &coords2);CHKERRQ(ierr); 2696 for (v = vStart; v < vEnd; ++v) { 2697 ierr = PetscSectionGetDof(coordSection, v, &dof);CHKERRQ(ierr); 2698 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 2699 ierr = PetscSectionGetOffset(cSection, v, &off2);CHKERRQ(ierr); 2700 for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d]; 2701 } 2702 ierr = DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);CHKERRQ(ierr); 2703 for (c = cStart; c < cEnd; ++c) { 2704 PetscScalar *cellCoords = NULL; 2705 PetscInt b; 2706 2707 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);CHKERRQ(ierr); 2708 ierr = PetscSectionGetOffset(cSection, c, &off2);CHKERRQ(ierr); 2709 for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b]; 2710 for (d = 0; d < dof/bs; ++d) {ierr = DMPlexLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);CHKERRQ(ierr);} 2711 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, &dof, &cellCoords);CHKERRQ(ierr); 2712 } 2713 ierr = DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);CHKERRQ(ierr); 2714 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 2715 ierr = VecRestoreArray(cVec, &coords2);CHKERRQ(ierr); 2716 ierr = DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);CHKERRQ(ierr); 2717 ierr = DMSetCoordinatesLocal(dm, cVec);CHKERRQ(ierr); 2718 ierr = VecDestroy(&cVec);CHKERRQ(ierr); 2719 ierr = PetscSectionDestroy(&cSection);CHKERRQ(ierr); 2720 PetscFunctionReturn(0); 2721 } 2722 2723 #undef __FUNCT__ 2724 #define __FUNCT__ "DMPlexGetDepthLabel" 2725 /*@ 2726 DMPlexGetDepthLabel - Get the DMLabel recording the depth of each point 2727 2728 Not Collective 2729 2730 Input Parameter: 2731 . dm - The DMPlex object 2732 2733 Output Parameter: 2734 . depthLabel - The DMLabel recording point depth 2735 2736 Level: developer 2737 2738 .keywords: mesh, points 2739 .seealso: DMPlexGetDepth(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum() 2740 @*/ 2741 PetscErrorCode DMPlexGetDepthLabel(DM dm, DMLabel *depthLabel) 2742 { 2743 DM_Plex *mesh = (DM_Plex*) dm->data; 2744 PetscErrorCode ierr; 2745 2746 PetscFunctionBegin; 2747 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2748 PetscValidPointer(depthLabel, 2); 2749 if (!mesh->depthLabel) {ierr = DMPlexGetLabel(dm, "depth", &mesh->depthLabel);CHKERRQ(ierr);} 2750 *depthLabel = mesh->depthLabel; 2751 PetscFunctionReturn(0); 2752 } 2753 2754 #undef __FUNCT__ 2755 #define __FUNCT__ "DMPlexGetDepth" 2756 /*@ 2757 DMPlexGetDepth - Get the depth of the DAG representing this mesh 2758 2759 Not Collective 2760 2761 Input Parameter: 2762 . dm - The DMPlex object 2763 2764 Output Parameter: 2765 . depth - The number of strata (breadth first levels) in the DAG 2766 2767 Level: developer 2768 2769 .keywords: mesh, points 2770 .seealso: DMPlexGetDepthLabel(), DMPlexGetHeightStratum(), DMPlexGetDepthStratum() 2771 @*/ 2772 PetscErrorCode DMPlexGetDepth(DM dm, PetscInt *depth) 2773 { 2774 DMLabel label; 2775 PetscInt d = 0; 2776 PetscErrorCode ierr; 2777 2778 PetscFunctionBegin; 2779 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2780 PetscValidPointer(depth, 2); 2781 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 2782 if (label) {ierr = DMLabelGetNumValues(label, &d);CHKERRQ(ierr);} 2783 *depth = d-1; 2784 PetscFunctionReturn(0); 2785 } 2786 2787 #undef __FUNCT__ 2788 #define __FUNCT__ "DMPlexGetDepthStratum" 2789 /*@ 2790 DMPlexGetDepthStratum - Get the bounds [start, end) for all points at a certain depth. 2791 2792 Not Collective 2793 2794 Input Parameters: 2795 + dm - The DMPlex object 2796 - stratumValue - The requested depth 2797 2798 Output Parameters: 2799 + start - The first point at this depth 2800 - end - One beyond the last point at this depth 2801 2802 Level: developer 2803 2804 .keywords: mesh, points 2805 .seealso: DMPlexGetHeightStratum(), DMPlexGetDepth() 2806 @*/ 2807 PetscErrorCode DMPlexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) 2808 { 2809 DMLabel label; 2810 PetscInt pStart, pEnd; 2811 PetscErrorCode ierr; 2812 2813 PetscFunctionBegin; 2814 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2815 if (start) {PetscValidPointer(start, 3); *start = 0;} 2816 if (end) {PetscValidPointer(end, 4); *end = 0;} 2817 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2818 if (pStart == pEnd) PetscFunctionReturn(0); 2819 if (stratumValue < 0) { 2820 if (start) *start = pStart; 2821 if (end) *end = pEnd; 2822 PetscFunctionReturn(0); 2823 } 2824 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 2825 if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found"); 2826 ierr = DMLabelGetStratumBounds(label, stratumValue, start, end);CHKERRQ(ierr); 2827 PetscFunctionReturn(0); 2828 } 2829 2830 #undef __FUNCT__ 2831 #define __FUNCT__ "DMPlexGetHeightStratum" 2832 /*@ 2833 DMPlexGetHeightStratum - Get the bounds [start, end) for all points at a certain height. 2834 2835 Not Collective 2836 2837 Input Parameters: 2838 + dm - The DMPlex object 2839 - stratumValue - The requested height 2840 2841 Output Parameters: 2842 + start - The first point at this height 2843 - end - One beyond the last point at this height 2844 2845 Level: developer 2846 2847 .keywords: mesh, points 2848 .seealso: DMPlexGetDepthStratum(), DMPlexGetDepth() 2849 @*/ 2850 PetscErrorCode DMPlexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) 2851 { 2852 DMLabel label; 2853 PetscInt depth, pStart, pEnd; 2854 PetscErrorCode ierr; 2855 2856 PetscFunctionBegin; 2857 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2858 if (start) {PetscValidPointer(start, 3); *start = 0;} 2859 if (end) {PetscValidPointer(end, 4); *end = 0;} 2860 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2861 if (pStart == pEnd) PetscFunctionReturn(0); 2862 if (stratumValue < 0) { 2863 if (start) *start = pStart; 2864 if (end) *end = pEnd; 2865 PetscFunctionReturn(0); 2866 } 2867 ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr); 2868 if (!label) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "No label named depth was found");CHKERRQ(ierr); 2869 ierr = DMLabelGetNumValues(label, &depth);CHKERRQ(ierr); 2870 ierr = DMLabelGetStratumBounds(label, depth-1-stratumValue, start, end);CHKERRQ(ierr); 2871 PetscFunctionReturn(0); 2872 } 2873 2874 #undef __FUNCT__ 2875 #define __FUNCT__ "DMPlexCreateSectionInitial" 2876 /* Set the number of dof on each point and separate by fields */ 2877 PetscErrorCode DMPlexCreateSectionInitial(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscSection *section) 2878 { 2879 PetscInt *numDofTot; 2880 PetscInt depth, pStart = 0, pEnd = 0; 2881 PetscInt p, d, dep, f; 2882 PetscErrorCode ierr; 2883 2884 PetscFunctionBegin; 2885 ierr = PetscMalloc1((dim+1), &numDofTot);CHKERRQ(ierr); 2886 for (d = 0; d <= dim; ++d) { 2887 numDofTot[d] = 0; 2888 for (f = 0; f < numFields; ++f) numDofTot[d] += numDof[f*(dim+1)+d]; 2889 } 2890 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);CHKERRQ(ierr); 2891 if (numFields > 0) { 2892 ierr = PetscSectionSetNumFields(*section, numFields);CHKERRQ(ierr); 2893 if (numComp) { 2894 for (f = 0; f < numFields; ++f) { 2895 ierr = PetscSectionSetFieldComponents(*section, f, numComp[f]);CHKERRQ(ierr); 2896 } 2897 } 2898 } 2899 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 2900 ierr = PetscSectionSetChart(*section, pStart, pEnd);CHKERRQ(ierr); 2901 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 2902 for (dep = 0; dep <= depth; ++dep) { 2903 d = dim == depth ? dep : (!dep ? 0 : dim); 2904 ierr = DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);CHKERRQ(ierr); 2905 for (p = pStart; p < pEnd; ++p) { 2906 for (f = 0; f < numFields; ++f) { 2907 ierr = PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);CHKERRQ(ierr); 2908 } 2909 ierr = PetscSectionSetDof(*section, p, numDofTot[d]);CHKERRQ(ierr); 2910 } 2911 } 2912 ierr = PetscFree(numDofTot);CHKERRQ(ierr); 2913 PetscFunctionReturn(0); 2914 } 2915 2916 #undef __FUNCT__ 2917 #define __FUNCT__ "DMPlexCreateSectionBCDof" 2918 /* Set the number of dof on each point and separate by fields 2919 If constDof is PETSC_DETERMINE, constrain every dof on the point 2920 */ 2921 PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], PetscInt constDof, PetscSection section) 2922 { 2923 PetscInt numFields; 2924 PetscInt bc; 2925 PetscErrorCode ierr; 2926 2927 PetscFunctionBegin; 2928 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 2929 for (bc = 0; bc < numBC; ++bc) { 2930 PetscInt field = 0; 2931 const PetscInt *idx; 2932 PetscInt n, i; 2933 2934 if (numFields) field = bcField[bc]; 2935 ierr = ISGetLocalSize(bcPoints[bc], &n);CHKERRQ(ierr); 2936 ierr = ISGetIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 2937 for (i = 0; i < n; ++i) { 2938 const PetscInt p = idx[i]; 2939 PetscInt numConst = constDof; 2940 2941 /* Constrain every dof on the point */ 2942 if (numConst < 0) { 2943 if (numFields) { 2944 ierr = PetscSectionGetFieldDof(section, p, field, &numConst);CHKERRQ(ierr); 2945 } else { 2946 ierr = PetscSectionGetDof(section, p, &numConst);CHKERRQ(ierr); 2947 } 2948 } 2949 if (numFields) { 2950 ierr = PetscSectionAddFieldConstraintDof(section, p, field, numConst);CHKERRQ(ierr); 2951 } 2952 ierr = PetscSectionAddConstraintDof(section, p, numConst);CHKERRQ(ierr); 2953 } 2954 ierr = ISRestoreIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 2955 } 2956 PetscFunctionReturn(0); 2957 } 2958 2959 #undef __FUNCT__ 2960 #define __FUNCT__ "DMPlexCreateSectionBCIndicesAll" 2961 /* Set the constrained indices on each point and separate by fields */ 2962 PetscErrorCode DMPlexCreateSectionBCIndicesAll(DM dm, PetscSection section) 2963 { 2964 PetscInt *maxConstraints; 2965 PetscInt numFields, f, pStart = 0, pEnd = 0, p; 2966 PetscErrorCode ierr; 2967 2968 PetscFunctionBegin; 2969 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 2970 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 2971 ierr = PetscMalloc1((numFields+1), &maxConstraints);CHKERRQ(ierr); 2972 for (f = 0; f <= numFields; ++f) maxConstraints[f] = 0; 2973 for (p = pStart; p < pEnd; ++p) { 2974 PetscInt cdof; 2975 2976 if (numFields) { 2977 for (f = 0; f < numFields; ++f) { 2978 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &cdof);CHKERRQ(ierr); 2979 maxConstraints[f] = PetscMax(maxConstraints[f], cdof); 2980 } 2981 } else { 2982 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 2983 maxConstraints[0] = PetscMax(maxConstraints[0], cdof); 2984 } 2985 } 2986 for (f = 0; f < numFields; ++f) { 2987 maxConstraints[numFields] += maxConstraints[f]; 2988 } 2989 if (maxConstraints[numFields]) { 2990 PetscInt *indices; 2991 2992 ierr = PetscMalloc1(maxConstraints[numFields], &indices);CHKERRQ(ierr); 2993 for (p = pStart; p < pEnd; ++p) { 2994 PetscInt cdof, d; 2995 2996 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 2997 if (cdof) { 2998 if (cdof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %D cDof %D > maxConstraints %D", p, cdof, maxConstraints[numFields]); 2999 if (numFields) { 3000 PetscInt numConst = 0, foff = 0; 3001 3002 for (f = 0; f < numFields; ++f) { 3003 PetscInt cfdof, fdof; 3004 3005 ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr); 3006 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);CHKERRQ(ierr); 3007 /* Change constraint numbering from absolute local dof number to field relative local dof number */ 3008 for (d = 0; d < cfdof; ++d) indices[numConst+d] = d; 3009 ierr = PetscSectionSetFieldConstraintIndices(section, p, f, &indices[numConst]);CHKERRQ(ierr); 3010 for (d = 0; d < cfdof; ++d) indices[numConst+d] += foff; 3011 numConst += cfdof; 3012 foff += fdof; 3013 } 3014 if (cdof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof); 3015 } else { 3016 for (d = 0; d < cdof; ++d) indices[d] = d; 3017 } 3018 ierr = PetscSectionSetConstraintIndices(section, p, indices);CHKERRQ(ierr); 3019 } 3020 } 3021 ierr = PetscFree(indices);CHKERRQ(ierr); 3022 } 3023 ierr = PetscFree(maxConstraints);CHKERRQ(ierr); 3024 PetscFunctionReturn(0); 3025 } 3026 3027 #undef __FUNCT__ 3028 #define __FUNCT__ "DMPlexCreateSectionBCIndicesField" 3029 /* Set the constrained field indices on each point */ 3030 PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt field, IS bcPoints, IS constraintIndices, PetscSection section) 3031 { 3032 const PetscInt *points, *indices; 3033 PetscInt numFields, maxDof, numPoints, p, numConstraints; 3034 PetscErrorCode ierr; 3035 3036 PetscFunctionBegin; 3037 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 3038 if ((field < 0) || (field >= numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, numFields); 3039 3040 ierr = ISGetLocalSize(bcPoints, &numPoints);CHKERRQ(ierr); 3041 ierr = ISGetIndices(bcPoints, &points);CHKERRQ(ierr); 3042 if (!constraintIndices) { 3043 PetscInt *idx, i; 3044 3045 ierr = PetscSectionGetMaxDof(section, &maxDof);CHKERRQ(ierr); 3046 ierr = PetscMalloc1(maxDof, &idx);CHKERRQ(ierr); 3047 for (i = 0; i < maxDof; ++i) idx[i] = i; 3048 for (p = 0; p < numPoints; ++p) { 3049 ierr = PetscSectionSetFieldConstraintIndices(section, points[p], field, idx);CHKERRQ(ierr); 3050 } 3051 ierr = PetscFree(idx);CHKERRQ(ierr); 3052 } else { 3053 ierr = ISGetLocalSize(constraintIndices, &numConstraints);CHKERRQ(ierr); 3054 ierr = ISGetIndices(constraintIndices, &indices);CHKERRQ(ierr); 3055 for (p = 0; p < numPoints; ++p) { 3056 PetscInt fcdof; 3057 3058 ierr = PetscSectionGetFieldConstraintDof(section, points[p], field, &fcdof);CHKERRQ(ierr); 3059 if (fcdof != numConstraints) SETERRQ4(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Section point %d field %d has %d constraints, but yo ugave %d indices", p, field, fcdof, numConstraints); 3060 ierr = PetscSectionSetFieldConstraintIndices(section, points[p], field, indices);CHKERRQ(ierr); 3061 } 3062 ierr = ISRestoreIndices(constraintIndices, &indices);CHKERRQ(ierr); 3063 } 3064 ierr = ISRestoreIndices(bcPoints, &points);CHKERRQ(ierr); 3065 PetscFunctionReturn(0); 3066 } 3067 3068 #undef __FUNCT__ 3069 #define __FUNCT__ "DMPlexCreateSectionBCIndices" 3070 /* Set the constrained indices on each point and separate by fields */ 3071 PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section) 3072 { 3073 PetscInt *indices; 3074 PetscInt numFields, maxDof, f, pStart = 0, pEnd = 0, p; 3075 PetscErrorCode ierr; 3076 3077 PetscFunctionBegin; 3078 ierr = PetscSectionGetMaxDof(section, &maxDof);CHKERRQ(ierr); 3079 ierr = PetscMalloc1(maxDof, &indices);CHKERRQ(ierr); 3080 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 3081 if (!numFields) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This function only works after users have set field constraint indices."); 3082 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3083 for (p = pStart; p < pEnd; ++p) { 3084 PetscInt cdof, d; 3085 3086 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 3087 if (cdof) { 3088 PetscInt numConst = 0, foff = 0; 3089 3090 for (f = 0; f < numFields; ++f) { 3091 const PetscInt *fcind; 3092 PetscInt fdof, fcdof; 3093 3094 ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr); 3095 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 3096 if (fcdof) {ierr = PetscSectionGetFieldConstraintIndices(section, p, f, &fcind);CHKERRQ(ierr);} 3097 /* Change constraint numbering from field relative local dof number to absolute local dof number */ 3098 for (d = 0; d < fcdof; ++d) indices[numConst+d] = fcind[d]+foff; 3099 foff += fdof; 3100 numConst += fcdof; 3101 } 3102 if (cdof != numConst) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cdof); 3103 ierr = PetscSectionSetConstraintIndices(section, p, indices);CHKERRQ(ierr); 3104 } 3105 } 3106 ierr = PetscFree(indices);CHKERRQ(ierr); 3107 PetscFunctionReturn(0); 3108 } 3109 3110 #undef __FUNCT__ 3111 #define __FUNCT__ "DMPlexCreateSection" 3112 /*@C 3113 DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided. 3114 3115 Not Collective 3116 3117 Input Parameters: 3118 + dm - The DMPlex object 3119 . dim - The spatial dimension of the problem 3120 . numFields - The number of fields in the problem 3121 . numComp - An array of size numFields that holds the number of components for each field 3122 . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d 3123 . numBC - The number of boundary conditions 3124 . bcField - An array of size numBC giving the field number for each boundry condition 3125 . bcPoints - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies 3126 - perm - Optional permutation of the chart, or NULL 3127 3128 Output Parameter: 3129 . section - The PetscSection object 3130 3131 Notes: numDof[f*(dim+1)+d] gives the number of dof for field f on sieve points of dimension d. For instance, numDof[1] is the 3132 number of dof for field 0 on each edge. 3133 3134 The chart permutation is the same one set using PetscSectionSetPermutation() 3135 3136 Level: developer 3137 3138 Fortran Notes: 3139 A Fortran 90 version is available as DMPlexCreateSectionF90() 3140 3141 .keywords: mesh, elements 3142 .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation() 3143 @*/ 3144 PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], IS perm, PetscSection *section) 3145 { 3146 PetscErrorCode ierr; 3147 3148 PetscFunctionBegin; 3149 ierr = DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);CHKERRQ(ierr); 3150 ierr = DMPlexCreateSectionBCDof(dm, numBC, bcField, bcPoints, PETSC_DETERMINE, *section);CHKERRQ(ierr); 3151 if (perm) {ierr = PetscSectionSetPermutation(*section, perm);CHKERRQ(ierr);} 3152 ierr = PetscSectionSetUp(*section);CHKERRQ(ierr); 3153 if (numBC) {ierr = DMPlexCreateSectionBCIndicesAll(dm, *section);CHKERRQ(ierr);} 3154 ierr = PetscSectionViewFromOptions(*section,NULL,"-section_view");CHKERRQ(ierr); 3155 PetscFunctionReturn(0); 3156 } 3157 3158 #undef __FUNCT__ 3159 #define __FUNCT__ "DMCreateCoordinateDM_Plex" 3160 PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm) 3161 { 3162 PetscSection section; 3163 PetscErrorCode ierr; 3164 3165 PetscFunctionBegin; 3166 ierr = DMClone(dm, cdm);CHKERRQ(ierr); 3167 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 3168 ierr = DMSetDefaultSection(*cdm, section);CHKERRQ(ierr); 3169 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 3170 PetscFunctionReturn(0); 3171 } 3172 3173 #undef __FUNCT__ 3174 #define __FUNCT__ "DMPlexGetConeSection" 3175 PetscErrorCode DMPlexGetConeSection(DM dm, PetscSection *section) 3176 { 3177 DM_Plex *mesh = (DM_Plex*) dm->data; 3178 3179 PetscFunctionBegin; 3180 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3181 if (section) *section = mesh->coneSection; 3182 PetscFunctionReturn(0); 3183 } 3184 3185 #undef __FUNCT__ 3186 #define __FUNCT__ "DMPlexGetSupportSection" 3187 PetscErrorCode DMPlexGetSupportSection(DM dm, PetscSection *section) 3188 { 3189 DM_Plex *mesh = (DM_Plex*) dm->data; 3190 3191 PetscFunctionBegin; 3192 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3193 if (section) *section = mesh->supportSection; 3194 PetscFunctionReturn(0); 3195 } 3196 3197 #undef __FUNCT__ 3198 #define __FUNCT__ "DMPlexGetCones" 3199 PetscErrorCode DMPlexGetCones(DM dm, PetscInt *cones[]) 3200 { 3201 DM_Plex *mesh = (DM_Plex*) dm->data; 3202 3203 PetscFunctionBegin; 3204 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3205 if (cones) *cones = mesh->cones; 3206 PetscFunctionReturn(0); 3207 } 3208 3209 #undef __FUNCT__ 3210 #define __FUNCT__ "DMPlexGetConeOrientations" 3211 PetscErrorCode DMPlexGetConeOrientations(DM dm, PetscInt *coneOrientations[]) 3212 { 3213 DM_Plex *mesh = (DM_Plex*) dm->data; 3214 3215 PetscFunctionBegin; 3216 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3217 if (coneOrientations) *coneOrientations = mesh->coneOrientations; 3218 PetscFunctionReturn(0); 3219 } 3220 3221 /******************************** FEM Support **********************************/ 3222 3223 #undef __FUNCT__ 3224 #define __FUNCT__ "DMPlexVecGetClosure_Depth1_Static" 3225 PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Depth1_Static(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[]) 3226 { 3227 PetscScalar *array, *vArray; 3228 const PetscInt *cone, *coneO; 3229 PetscInt pStart, pEnd, p, numPoints, size = 0, offset = 0; 3230 PetscErrorCode ierr; 3231 3232 PetscFunctionBeginHot; 3233 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3234 ierr = DMPlexGetConeSize(dm, point, &numPoints);CHKERRQ(ierr); 3235 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3236 ierr = DMPlexGetConeOrientation(dm, point, &coneO);CHKERRQ(ierr); 3237 if (!values || !*values) { 3238 if ((point >= pStart) && (point < pEnd)) { 3239 PetscInt dof; 3240 3241 ierr = PetscSectionGetDof(section, point, &dof);CHKERRQ(ierr); 3242 size += dof; 3243 } 3244 for (p = 0; p < numPoints; ++p) { 3245 const PetscInt cp = cone[p]; 3246 PetscInt dof; 3247 3248 if ((cp < pStart) || (cp >= pEnd)) continue; 3249 ierr = PetscSectionGetDof(section, cp, &dof);CHKERRQ(ierr); 3250 size += dof; 3251 } 3252 if (!values) { 3253 if (csize) *csize = size; 3254 PetscFunctionReturn(0); 3255 } 3256 ierr = DMGetWorkArray(dm, size, PETSC_SCALAR, &array);CHKERRQ(ierr); 3257 } else { 3258 array = *values; 3259 } 3260 size = 0; 3261 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 3262 if ((point >= pStart) && (point < pEnd)) { 3263 PetscInt dof, off, d; 3264 PetscScalar *varr; 3265 3266 ierr = PetscSectionGetDof(section, point, &dof);CHKERRQ(ierr); 3267 ierr = PetscSectionGetOffset(section, point, &off);CHKERRQ(ierr); 3268 varr = &vArray[off]; 3269 for (d = 0; d < dof; ++d, ++offset) { 3270 array[offset] = varr[d]; 3271 } 3272 size += dof; 3273 } 3274 for (p = 0; p < numPoints; ++p) { 3275 const PetscInt cp = cone[p]; 3276 PetscInt o = coneO[p]; 3277 PetscInt dof, off, d; 3278 PetscScalar *varr; 3279 3280 if ((cp < pStart) || (cp >= pEnd)) continue; 3281 ierr = PetscSectionGetDof(section, cp, &dof);CHKERRQ(ierr); 3282 ierr = PetscSectionGetOffset(section, cp, &off);CHKERRQ(ierr); 3283 varr = &vArray[off]; 3284 if (o >= 0) { 3285 for (d = 0; d < dof; ++d, ++offset) { 3286 array[offset] = varr[d]; 3287 } 3288 } else { 3289 for (d = dof-1; d >= 0; --d, ++offset) { 3290 array[offset] = varr[d]; 3291 } 3292 } 3293 size += dof; 3294 } 3295 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 3296 if (!*values) { 3297 if (csize) *csize = size; 3298 *values = array; 3299 } else { 3300 if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size); 3301 *csize = size; 3302 } 3303 PetscFunctionReturn(0); 3304 } 3305 3306 #undef __FUNCT__ 3307 #define __FUNCT__ "DMPlexVecGetClosure_Static" 3308 PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], const PetscScalar vArray[], PetscInt *size, PetscScalar array[]) 3309 { 3310 PetscInt offset = 0, p; 3311 PetscErrorCode ierr; 3312 3313 PetscFunctionBeginHot; 3314 *size = 0; 3315 for (p = 0; p < numPoints*2; p += 2) { 3316 const PetscInt point = points[p]; 3317 const PetscInt o = points[p+1]; 3318 PetscInt dof, off, d; 3319 const PetscScalar *varr; 3320 3321 ierr = PetscSectionGetDof(section, point, &dof);CHKERRQ(ierr); 3322 ierr = PetscSectionGetOffset(section, point, &off);CHKERRQ(ierr); 3323 varr = &vArray[off]; 3324 if (o >= 0) { 3325 for (d = 0; d < dof; ++d, ++offset) array[offset] = varr[d]; 3326 } else { 3327 for (d = dof-1; d >= 0; --d, ++offset) array[offset] = varr[d]; 3328 } 3329 } 3330 *size = offset; 3331 PetscFunctionReturn(0); 3332 } 3333 3334 #undef __FUNCT__ 3335 #define __FUNCT__ "DMPlexVecGetClosure_Fields_Static" 3336 PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[]) 3337 { 3338 PetscInt offset = 0, f; 3339 PetscErrorCode ierr; 3340 3341 PetscFunctionBeginHot; 3342 *size = 0; 3343 for (f = 0; f < numFields; ++f) { 3344 PetscInt fcomp, p; 3345 3346 ierr = PetscSectionGetFieldComponents(section, f, &fcomp);CHKERRQ(ierr); 3347 for (p = 0; p < numPoints*2; p += 2) { 3348 const PetscInt point = points[p]; 3349 const PetscInt o = points[p+1]; 3350 PetscInt fdof, foff, d, c; 3351 const PetscScalar *varr; 3352 3353 ierr = PetscSectionGetFieldDof(section, point, f, &fdof);CHKERRQ(ierr); 3354 ierr = PetscSectionGetFieldOffset(section, point, f, &foff);CHKERRQ(ierr); 3355 varr = &vArray[foff]; 3356 if (o >= 0) { 3357 for (d = 0; d < fdof; ++d, ++offset) array[offset] = varr[d]; 3358 } else { 3359 for (d = fdof/fcomp-1; d >= 0; --d) { 3360 for (c = 0; c < fcomp; ++c, ++offset) { 3361 array[offset] = varr[d*fcomp+c]; 3362 } 3363 } 3364 } 3365 } 3366 } 3367 *size = offset; 3368 PetscFunctionReturn(0); 3369 } 3370 3371 #undef __FUNCT__ 3372 #define __FUNCT__ "DMPlexVecGetClosure" 3373 /*@C 3374 DMPlexVecGetClosure - Get an array of the values on the closure of 'point' 3375 3376 Not collective 3377 3378 Input Parameters: 3379 + dm - The DM 3380 . section - The section describing the layout in v, or NULL to use the default section 3381 . v - The local vector 3382 - point - The sieve point in the DM 3383 3384 Output Parameters: 3385 + csize - The number of values in the closure, or NULL 3386 - values - The array of values, which is a borrowed array and should not be freed 3387 3388 Fortran Notes: 3389 Since it returns an array, this routine is only available in Fortran 90, and you must 3390 include petsc.h90 in your code. 3391 3392 The csize argument is not present in the Fortran 90 binding since it is internal to the array. 3393 3394 Level: intermediate 3395 3396 .seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure() 3397 @*/ 3398 PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[]) 3399 { 3400 PetscSection clSection; 3401 IS clPoints; 3402 PetscScalar *array, *vArray; 3403 PetscInt *points = NULL; 3404 const PetscInt *clp; 3405 PetscInt depth, numFields, numPoints, size; 3406 PetscErrorCode ierr; 3407 3408 PetscFunctionBeginHot; 3409 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3410 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 3411 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 3412 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 3413 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3414 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 3415 if (depth == 1 && numFields < 2) { 3416 ierr = DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);CHKERRQ(ierr); 3417 PetscFunctionReturn(0); 3418 } 3419 /* Get points */ 3420 ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr); 3421 if (!clPoints) { 3422 PetscInt pStart, pEnd, p, q; 3423 3424 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3425 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 3426 /* Compress out points not in the section */ 3427 for (p = 0, q = 0; p < numPoints*2; p += 2) { 3428 if ((points[p] >= pStart) && (points[p] < pEnd)) { 3429 points[q*2] = points[p]; 3430 points[q*2+1] = points[p+1]; 3431 ++q; 3432 } 3433 } 3434 numPoints = q; 3435 } else { 3436 PetscInt dof, off; 3437 3438 ierr = PetscSectionGetDof(clSection, point, &dof);CHKERRQ(ierr); 3439 ierr = PetscSectionGetOffset(clSection, point, &off);CHKERRQ(ierr); 3440 ierr = ISGetIndices(clPoints, &clp);CHKERRQ(ierr); 3441 numPoints = dof/2; 3442 points = (PetscInt *) &clp[off]; 3443 } 3444 /* Get array */ 3445 if (!values || !*values) { 3446 PetscInt asize = 0, dof, p; 3447 3448 for (p = 0; p < numPoints*2; p += 2) { 3449 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3450 asize += dof; 3451 } 3452 if (!values) { 3453 if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);} 3454 else {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);} 3455 if (csize) *csize = asize; 3456 PetscFunctionReturn(0); 3457 } 3458 ierr = DMGetWorkArray(dm, asize, PETSC_SCALAR, &array);CHKERRQ(ierr); 3459 } else { 3460 array = *values; 3461 } 3462 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 3463 /* Get values */ 3464 if (numFields > 0) {ierr = DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);CHKERRQ(ierr);} 3465 else {ierr = DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);CHKERRQ(ierr);} 3466 /* Cleanup points */ 3467 if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);} 3468 else {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);} 3469 /* Cleanup array */ 3470 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 3471 if (!*values) { 3472 if (csize) *csize = size; 3473 *values = array; 3474 } else { 3475 if (size > *csize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Size of input array %d < actual size %d", *csize, size); 3476 *csize = size; 3477 } 3478 PetscFunctionReturn(0); 3479 } 3480 3481 #undef __FUNCT__ 3482 #define __FUNCT__ "DMPlexVecRestoreClosure" 3483 /*@C 3484 DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point' 3485 3486 Not collective 3487 3488 Input Parameters: 3489 + dm - The DM 3490 . section - The section describing the layout in v, or NULL to use the default section 3491 . v - The local vector 3492 . point - The sieve point in the DM 3493 . csize - The number of values in the closure, or NULL 3494 - values - The array of values, which is a borrowed array and should not be freed 3495 3496 Fortran Notes: 3497 Since it returns an array, this routine is only available in Fortran 90, and you must 3498 include petsc.h90 in your code. 3499 3500 The csize argument is not present in the Fortran 90 binding since it is internal to the array. 3501 3502 Level: intermediate 3503 3504 .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure() 3505 @*/ 3506 PetscErrorCode DMPlexVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[]) 3507 { 3508 PetscInt size = 0; 3509 PetscErrorCode ierr; 3510 3511 PetscFunctionBegin; 3512 /* Should work without recalculating size */ 3513 ierr = DMRestoreWorkArray(dm, size, PETSC_SCALAR, (void*) values);CHKERRQ(ierr); 3514 PetscFunctionReturn(0); 3515 } 3516 3517 PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;} 3518 PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;} 3519 3520 #undef __FUNCT__ 3521 #define __FUNCT__ "updatePoint_private" 3522 PETSC_STATIC_INLINE PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[]) 3523 { 3524 PetscInt cdof; /* The number of constraints on this point */ 3525 const PetscInt *cdofs; /* The indices of the constrained dofs on this point */ 3526 PetscScalar *a; 3527 PetscInt off, cind = 0, k; 3528 PetscErrorCode ierr; 3529 3530 PetscFunctionBegin; 3531 ierr = PetscSectionGetConstraintDof(section, point, &cdof);CHKERRQ(ierr); 3532 ierr = PetscSectionGetOffset(section, point, &off);CHKERRQ(ierr); 3533 a = &array[off]; 3534 if (!cdof || setBC) { 3535 if (orientation >= 0) { 3536 for (k = 0; k < dof; ++k) { 3537 fuse(&a[k], values[k]); 3538 } 3539 } else { 3540 for (k = 0; k < dof; ++k) { 3541 fuse(&a[k], values[dof-k-1]); 3542 } 3543 } 3544 } else { 3545 ierr = PetscSectionGetConstraintIndices(section, point, &cdofs);CHKERRQ(ierr); 3546 if (orientation >= 0) { 3547 for (k = 0; k < dof; ++k) { 3548 if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;} 3549 fuse(&a[k], values[k]); 3550 } 3551 } else { 3552 for (k = 0; k < dof; ++k) { 3553 if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;} 3554 fuse(&a[k], values[dof-k-1]); 3555 } 3556 } 3557 } 3558 PetscFunctionReturn(0); 3559 } 3560 3561 #undef __FUNCT__ 3562 #define __FUNCT__ "updatePointBC_private" 3563 PETSC_STATIC_INLINE PetscErrorCode updatePointBC_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar*, PetscScalar), PetscInt orientation, const PetscScalar values[], PetscScalar array[]) 3564 { 3565 PetscInt cdof; /* The number of constraints on this point */ 3566 const PetscInt *cdofs; /* The indices of the constrained dofs on this point */ 3567 PetscScalar *a; 3568 PetscInt off, cind = 0, k; 3569 PetscErrorCode ierr; 3570 3571 PetscFunctionBegin; 3572 ierr = PetscSectionGetConstraintDof(section, point, &cdof);CHKERRQ(ierr); 3573 ierr = PetscSectionGetOffset(section, point, &off);CHKERRQ(ierr); 3574 a = &array[off]; 3575 if (cdof) { 3576 ierr = PetscSectionGetConstraintIndices(section, point, &cdofs);CHKERRQ(ierr); 3577 if (orientation >= 0) { 3578 for (k = 0; k < dof; ++k) { 3579 if ((cind < cdof) && (k == cdofs[cind])) { 3580 fuse(&a[k], values[k]); 3581 ++cind; 3582 } 3583 } 3584 } else { 3585 for (k = 0; k < dof; ++k) { 3586 if ((cind < cdof) && (k == cdofs[cind])) { 3587 fuse(&a[k], values[dof-k-1]); 3588 ++cind; 3589 } 3590 } 3591 } 3592 } 3593 PetscFunctionReturn(0); 3594 } 3595 3596 #undef __FUNCT__ 3597 #define __FUNCT__ "updatePointFields_private" 3598 PETSC_STATIC_INLINE PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), PetscBool setBC, const PetscScalar values[], PetscInt *offset, PetscScalar array[]) 3599 { 3600 PetscScalar *a; 3601 PetscInt fdof, foff, fcdof, foffset = *offset; 3602 const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */ 3603 PetscInt cind = 0, k, c; 3604 PetscErrorCode ierr; 3605 3606 PetscFunctionBegin; 3607 ierr = PetscSectionGetFieldDof(section, point, f, &fdof);CHKERRQ(ierr); 3608 ierr = PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);CHKERRQ(ierr); 3609 ierr = PetscSectionGetFieldOffset(section, point, f, &foff);CHKERRQ(ierr); 3610 a = &array[foff]; 3611 if (!fcdof || setBC) { 3612 if (o >= 0) { 3613 for (k = 0; k < fdof; ++k) fuse(&a[k], values[foffset+k]); 3614 } else { 3615 for (k = fdof/fcomp-1; k >= 0; --k) { 3616 for (c = 0; c < fcomp; ++c) { 3617 fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]); 3618 } 3619 } 3620 } 3621 } else { 3622 ierr = PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);CHKERRQ(ierr); 3623 if (o >= 0) { 3624 for (k = 0; k < fdof; ++k) { 3625 if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;} 3626 fuse(&a[k], values[foffset+k]); 3627 } 3628 } else { 3629 for (k = fdof/fcomp-1; k >= 0; --k) { 3630 for (c = 0; c < fcomp; ++c) { 3631 if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;} 3632 fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]); 3633 } 3634 } 3635 } 3636 } 3637 *offset += fdof; 3638 PetscFunctionReturn(0); 3639 } 3640 3641 #undef __FUNCT__ 3642 #define __FUNCT__ "updatePointFieldsBC_private" 3643 PETSC_STATIC_INLINE PetscErrorCode updatePointFieldsBC_private(PetscSection section, PetscInt point, PetscInt o, PetscInt f, PetscInt fcomp, void (*fuse)(PetscScalar*, PetscScalar), const PetscScalar values[], PetscInt *offset, PetscScalar array[]) 3644 { 3645 PetscScalar *a; 3646 PetscInt fdof, foff, fcdof, foffset = *offset; 3647 const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */ 3648 PetscInt cind = 0, k, c; 3649 PetscErrorCode ierr; 3650 3651 PetscFunctionBegin; 3652 ierr = PetscSectionGetFieldDof(section, point, f, &fdof);CHKERRQ(ierr); 3653 ierr = PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);CHKERRQ(ierr); 3654 ierr = PetscSectionGetFieldOffset(section, point, f, &foff);CHKERRQ(ierr); 3655 a = &array[foff]; 3656 if (fcdof) { 3657 ierr = PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);CHKERRQ(ierr); 3658 if (o >= 0) { 3659 for (k = 0; k < fdof; ++k) { 3660 if ((cind < fcdof) && (k == fcdofs[cind])) { 3661 fuse(&a[k], values[foffset+k]); 3662 ++cind; 3663 } 3664 } 3665 } else { 3666 for (k = fdof/fcomp-1; k >= 0; --k) { 3667 for (c = 0; c < fcomp; ++c) { 3668 if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) { 3669 fuse(&a[(fdof/fcomp-1-k)*fcomp+c], values[foffset+k*fcomp+c]); 3670 ++cind; 3671 } 3672 } 3673 } 3674 } 3675 } 3676 *offset += fdof; 3677 PetscFunctionReturn(0); 3678 } 3679 3680 #undef __FUNCT__ 3681 #define __FUNCT__ "DMPlexVecSetClosure_Static" 3682 PETSC_STATIC_INLINE PetscErrorCode DMPlexVecSetClosure_Static(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode) 3683 { 3684 PetscScalar *array; 3685 const PetscInt *cone, *coneO; 3686 PetscInt pStart, pEnd, p, numPoints, off, dof; 3687 PetscErrorCode ierr; 3688 3689 PetscFunctionBeginHot; 3690 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3691 ierr = DMPlexGetConeSize(dm, point, &numPoints);CHKERRQ(ierr); 3692 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 3693 ierr = DMPlexGetConeOrientation(dm, point, &coneO);CHKERRQ(ierr); 3694 ierr = VecGetArray(v, &array);CHKERRQ(ierr); 3695 for (p = 0, off = 0; p <= numPoints; ++p, off += dof) { 3696 const PetscInt cp = !p ? point : cone[p-1]; 3697 const PetscInt o = !p ? 0 : coneO[p-1]; 3698 3699 if ((cp < pStart) || (cp >= pEnd)) {dof = 0; continue;} 3700 ierr = PetscSectionGetDof(section, cp, &dof);CHKERRQ(ierr); 3701 /* ADD_VALUES */ 3702 { 3703 const PetscInt *cdofs; /* The indices of the constrained dofs on this point */ 3704 PetscScalar *a; 3705 PetscInt cdof, coff, cind = 0, k; 3706 3707 ierr = PetscSectionGetConstraintDof(section, cp, &cdof);CHKERRQ(ierr); 3708 ierr = PetscSectionGetOffset(section, cp, &coff);CHKERRQ(ierr); 3709 a = &array[coff]; 3710 if (!cdof) { 3711 if (o >= 0) { 3712 for (k = 0; k < dof; ++k) { 3713 a[k] += values[off+k]; 3714 } 3715 } else { 3716 for (k = 0; k < dof; ++k) { 3717 a[k] += values[off+dof-k-1]; 3718 } 3719 } 3720 } else { 3721 ierr = PetscSectionGetConstraintIndices(section, cp, &cdofs);CHKERRQ(ierr); 3722 if (o >= 0) { 3723 for (k = 0; k < dof; ++k) { 3724 if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;} 3725 a[k] += values[off+k]; 3726 } 3727 } else { 3728 for (k = 0; k < dof; ++k) { 3729 if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;} 3730 a[k] += values[off+dof-k-1]; 3731 } 3732 } 3733 } 3734 } 3735 } 3736 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr); 3737 PetscFunctionReturn(0); 3738 } 3739 3740 #undef __FUNCT__ 3741 #define __FUNCT__ "DMPlexVecSetClosure" 3742 /*@C 3743 DMPlexVecSetClosure - Set an array of the values on the closure of 'point' 3744 3745 Not collective 3746 3747 Input Parameters: 3748 + dm - The DM 3749 . section - The section describing the layout in v, or NULL to use the default section 3750 . v - The local vector 3751 . point - The sieve point in the DM 3752 . values - The array of values 3753 - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions 3754 3755 Fortran Notes: 3756 This routine is only available in Fortran 90, and you must include petsc.h90 in your code. 3757 3758 Level: intermediate 3759 3760 .seealso DMPlexVecGetClosure(), DMPlexMatSetClosure() 3761 @*/ 3762 PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode) 3763 { 3764 PetscSection clSection; 3765 IS clPoints; 3766 PetscScalar *array; 3767 PetscInt *points = NULL; 3768 const PetscInt *clp; 3769 PetscInt depth, numFields, numPoints, p; 3770 PetscErrorCode ierr; 3771 3772 PetscFunctionBeginHot; 3773 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3774 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 3775 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 3776 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 3777 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 3778 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 3779 if (depth == 1 && numFields < 2 && mode == ADD_VALUES) { 3780 ierr = DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);CHKERRQ(ierr); 3781 PetscFunctionReturn(0); 3782 } 3783 /* Get points */ 3784 ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr); 3785 if (!clPoints) { 3786 PetscInt pStart, pEnd, q; 3787 3788 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3789 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 3790 /* Compress out points not in the section */ 3791 for (p = 0, q = 0; p < numPoints*2; p += 2) { 3792 if ((points[p] >= pStart) && (points[p] < pEnd)) { 3793 points[q*2] = points[p]; 3794 points[q*2+1] = points[p+1]; 3795 ++q; 3796 } 3797 } 3798 numPoints = q; 3799 } else { 3800 PetscInt dof, off; 3801 3802 ierr = PetscSectionGetDof(clSection, point, &dof);CHKERRQ(ierr); 3803 ierr = PetscSectionGetOffset(clSection, point, &off);CHKERRQ(ierr); 3804 ierr = ISGetIndices(clPoints, &clp);CHKERRQ(ierr); 3805 numPoints = dof/2; 3806 points = (PetscInt *) &clp[off]; 3807 } 3808 /* Get array */ 3809 ierr = VecGetArray(v, &array);CHKERRQ(ierr); 3810 /* Get values */ 3811 if (numFields > 0) { 3812 PetscInt offset = 0, fcomp, f; 3813 for (f = 0; f < numFields; ++f) { 3814 ierr = PetscSectionGetFieldComponents(section, f, &fcomp);CHKERRQ(ierr); 3815 switch (mode) { 3816 case INSERT_VALUES: 3817 for (p = 0; p < numPoints*2; p += 2) { 3818 const PetscInt point = points[p]; 3819 const PetscInt o = points[p+1]; 3820 updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array); 3821 } break; 3822 case INSERT_ALL_VALUES: 3823 for (p = 0; p < numPoints*2; p += 2) { 3824 const PetscInt point = points[p]; 3825 const PetscInt o = points[p+1]; 3826 updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array); 3827 } break; 3828 case INSERT_BC_VALUES: 3829 for (p = 0; p < numPoints*2; p += 2) { 3830 const PetscInt point = points[p]; 3831 const PetscInt o = points[p+1]; 3832 updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array); 3833 } break; 3834 case ADD_VALUES: 3835 for (p = 0; p < numPoints*2; p += 2) { 3836 const PetscInt point = points[p]; 3837 const PetscInt o = points[p+1]; 3838 updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array); 3839 } break; 3840 case ADD_ALL_VALUES: 3841 for (p = 0; p < numPoints*2; p += 2) { 3842 const PetscInt point = points[p]; 3843 const PetscInt o = points[p+1]; 3844 updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array); 3845 } break; 3846 default: 3847 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode); 3848 } 3849 } 3850 } else { 3851 PetscInt dof, off; 3852 3853 switch (mode) { 3854 case INSERT_VALUES: 3855 for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) { 3856 PetscInt o = points[p+1]; 3857 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3858 updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array); 3859 } break; 3860 case INSERT_ALL_VALUES: 3861 for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) { 3862 PetscInt o = points[p+1]; 3863 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3864 updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array); 3865 } break; 3866 case INSERT_BC_VALUES: 3867 for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) { 3868 PetscInt o = points[p+1]; 3869 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3870 updatePointBC_private(section, points[p], dof, insert, o, &values[off], array); 3871 } break; 3872 case ADD_VALUES: 3873 for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) { 3874 PetscInt o = points[p+1]; 3875 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3876 updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array); 3877 } break; 3878 case ADD_ALL_VALUES: 3879 for (p = 0, off = 0; p < numPoints*2; p += 2, off += dof) { 3880 PetscInt o = points[p+1]; 3881 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 3882 updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array); 3883 } break; 3884 default: 3885 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode); 3886 } 3887 } 3888 /* Cleanup points */ 3889 if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);} 3890 else {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);} 3891 /* Cleanup array */ 3892 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr); 3893 PetscFunctionReturn(0); 3894 } 3895 3896 #undef __FUNCT__ 3897 #define __FUNCT__ "DMPlexVecSetFieldClosure_Internal" 3898 PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM dm, PetscSection section, Vec v, PetscBool fieldActive[], PetscInt point, const PetscScalar values[], InsertMode mode) 3899 { 3900 PetscSection clSection; 3901 IS clPoints; 3902 PetscScalar *array; 3903 PetscInt *points = NULL; 3904 const PetscInt *clp; 3905 PetscInt numFields, numPoints, p; 3906 PetscInt offset = 0, fcomp, f; 3907 PetscErrorCode ierr; 3908 3909 PetscFunctionBeginHot; 3910 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3911 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 3912 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 3913 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 3914 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 3915 /* Get points */ 3916 ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr); 3917 if (!clPoints) { 3918 PetscInt pStart, pEnd, q; 3919 3920 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 3921 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 3922 /* Compress out points not in the section */ 3923 for (p = 0, q = 0; p < numPoints*2; p += 2) { 3924 if ((points[p] >= pStart) && (points[p] < pEnd)) { 3925 points[q*2] = points[p]; 3926 points[q*2+1] = points[p+1]; 3927 ++q; 3928 } 3929 } 3930 numPoints = q; 3931 } else { 3932 PetscInt dof, off; 3933 3934 ierr = PetscSectionGetDof(clSection, point, &dof);CHKERRQ(ierr); 3935 ierr = PetscSectionGetOffset(clSection, point, &off);CHKERRQ(ierr); 3936 ierr = ISGetIndices(clPoints, &clp);CHKERRQ(ierr); 3937 numPoints = dof/2; 3938 points = (PetscInt *) &clp[off]; 3939 } 3940 /* Get array */ 3941 ierr = VecGetArray(v, &array);CHKERRQ(ierr); 3942 /* Get values */ 3943 for (f = 0; f < numFields; ++f) { 3944 ierr = PetscSectionGetFieldComponents(section, f, &fcomp);CHKERRQ(ierr); 3945 if (!fieldActive[f]) { 3946 for (p = 0; p < numPoints*2; p += 2) { 3947 PetscInt fdof; 3948 ierr = PetscSectionGetFieldDof(section, points[p], f, &fdof);CHKERRQ(ierr); 3949 offset += fdof; 3950 } 3951 continue; 3952 } 3953 switch (mode) { 3954 case INSERT_VALUES: 3955 for (p = 0; p < numPoints*2; p += 2) { 3956 const PetscInt point = points[p]; 3957 const PetscInt o = points[p+1]; 3958 updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_FALSE, values, &offset, array); 3959 } break; 3960 case INSERT_ALL_VALUES: 3961 for (p = 0; p < numPoints*2; p += 2) { 3962 const PetscInt point = points[p]; 3963 const PetscInt o = points[p+1]; 3964 updatePointFields_private(section, point, o, f, fcomp, insert, PETSC_TRUE, values, &offset, array); 3965 } break; 3966 case INSERT_BC_VALUES: 3967 for (p = 0; p < numPoints*2; p += 2) { 3968 const PetscInt point = points[p]; 3969 const PetscInt o = points[p+1]; 3970 updatePointFieldsBC_private(section, point, o, f, fcomp, insert, values, &offset, array); 3971 } break; 3972 case ADD_VALUES: 3973 for (p = 0; p < numPoints*2; p += 2) { 3974 const PetscInt point = points[p]; 3975 const PetscInt o = points[p+1]; 3976 updatePointFields_private(section, point, o, f, fcomp, add, PETSC_FALSE, values, &offset, array); 3977 } break; 3978 case ADD_ALL_VALUES: 3979 for (p = 0; p < numPoints*2; p += 2) { 3980 const PetscInt point = points[p]; 3981 const PetscInt o = points[p+1]; 3982 updatePointFields_private(section, point, o, f, fcomp, add, PETSC_TRUE, values, &offset, array); 3983 } break; 3984 default: 3985 SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode); 3986 } 3987 } 3988 /* Cleanup points */ 3989 if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);} 3990 else {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);} 3991 /* Cleanup array */ 3992 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr); 3993 PetscFunctionReturn(0); 3994 } 3995 3996 #undef __FUNCT__ 3997 #define __FUNCT__ "DMPlexPrintMatSetValues" 3998 PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numRIndices, const PetscInt rindices[], PetscInt numCIndices, const PetscInt cindices[], const PetscScalar values[]) 3999 { 4000 PetscMPIInt rank; 4001 PetscInt i, j; 4002 PetscErrorCode ierr; 4003 4004 PetscFunctionBegin; 4005 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr); 4006 ierr = PetscViewerASCIIPrintf(viewer, "[%D]mat for sieve point %D\n", rank, point);CHKERRQ(ierr); 4007 for (i = 0; i < numRIndices; i++) {ierr = PetscViewerASCIIPrintf(viewer, "[%D]mat row indices[%D] = %D\n", rank, i, rindices[i]);CHKERRQ(ierr);} 4008 for (i = 0; i < numCIndices; i++) {ierr = PetscViewerASCIIPrintf(viewer, "[%D]mat col indices[%D] = %D\n", rank, i, cindices[i]);CHKERRQ(ierr);} 4009 numCIndices = numCIndices ? numCIndices : numRIndices; 4010 for (i = 0; i < numRIndices; i++) { 4011 ierr = PetscViewerASCIIPrintf(viewer, "[%D]", rank);CHKERRQ(ierr); 4012 for (j = 0; j < numCIndices; j++) { 4013 #if defined(PETSC_USE_COMPLEX) 4014 ierr = PetscViewerASCIIPrintf(viewer, " (%g,%g)", (double)PetscRealPart(values[i*numCIndices+j]), (double)PetscImaginaryPart(values[i*numCIndices+j]));CHKERRQ(ierr); 4015 #else 4016 ierr = PetscViewerASCIIPrintf(viewer, " %g", (double)values[i*numCIndices+j]);CHKERRQ(ierr); 4017 #endif 4018 } 4019 ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); 4020 } 4021 PetscFunctionReturn(0); 4022 } 4023 4024 #undef __FUNCT__ 4025 #define __FUNCT__ "indicesPoint_private" 4026 /* . off - The global offset of this point */ 4027 PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt off, PetscInt *loff, PetscBool setBC, PetscInt orientation, PetscInt indices[]) 4028 { 4029 PetscInt dof; /* The number of unknowns on this point */ 4030 PetscInt cdof; /* The number of constraints on this point */ 4031 const PetscInt *cdofs; /* The indices of the constrained dofs on this point */ 4032 PetscInt cind = 0, k; 4033 PetscErrorCode ierr; 4034 4035 PetscFunctionBegin; 4036 ierr = PetscSectionGetDof(section, point, &dof);CHKERRQ(ierr); 4037 ierr = PetscSectionGetConstraintDof(section, point, &cdof);CHKERRQ(ierr); 4038 if (!cdof || setBC) { 4039 if (orientation >= 0) { 4040 for (k = 0; k < dof; ++k) indices[*loff+k] = off+k; 4041 } else { 4042 for (k = 0; k < dof; ++k) indices[*loff+dof-k-1] = off+k; 4043 } 4044 } else { 4045 ierr = PetscSectionGetConstraintIndices(section, point, &cdofs);CHKERRQ(ierr); 4046 if (orientation >= 0) { 4047 for (k = 0; k < dof; ++k) { 4048 if ((cind < cdof) && (k == cdofs[cind])) { 4049 /* Insert check for returning constrained indices */ 4050 indices[*loff+k] = -(off+k+1); 4051 ++cind; 4052 } else { 4053 indices[*loff+k] = off+k-cind; 4054 } 4055 } 4056 } else { 4057 for (k = 0; k < dof; ++k) { 4058 if ((cind < cdof) && (k == cdofs[cind])) { 4059 /* Insert check for returning constrained indices */ 4060 indices[*loff+dof-k-1] = -(off+k+1); 4061 ++cind; 4062 } else { 4063 indices[*loff+dof-k-1] = off+k-cind; 4064 } 4065 } 4066 } 4067 } 4068 *loff += dof; 4069 PetscFunctionReturn(0); 4070 } 4071 4072 #undef __FUNCT__ 4073 #define __FUNCT__ "indicesPointFields_private" 4074 /* . off - The global offset of this point */ 4075 PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[]) 4076 { 4077 PetscInt numFields, foff, f; 4078 PetscErrorCode ierr; 4079 4080 PetscFunctionBegin; 4081 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 4082 for (f = 0, foff = 0; f < numFields; ++f) { 4083 PetscInt fdof, fcomp, cfdof; 4084 const PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */ 4085 PetscInt cind = 0, k, c; 4086 4087 ierr = PetscSectionGetFieldComponents(section, f, &fcomp);CHKERRQ(ierr); 4088 ierr = PetscSectionGetFieldDof(section, point, f, &fdof);CHKERRQ(ierr); 4089 ierr = PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);CHKERRQ(ierr); 4090 if (!cfdof || setBC) { 4091 if (orientation >= 0) { 4092 for (k = 0; k < fdof; ++k) indices[foffs[f]+k] = off+foff+k; 4093 } else { 4094 for (k = fdof/fcomp-1; k >= 0; --k) { 4095 for (c = 0; c < fcomp; ++c) { 4096 indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c; 4097 } 4098 } 4099 } 4100 } else { 4101 ierr = PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);CHKERRQ(ierr); 4102 if (orientation >= 0) { 4103 for (k = 0; k < fdof; ++k) { 4104 if ((cind < cfdof) && (k == fcdofs[cind])) { 4105 indices[foffs[f]+k] = -(off+foff+k+1); 4106 ++cind; 4107 } else { 4108 indices[foffs[f]+k] = off+foff+k-cind; 4109 } 4110 } 4111 } else { 4112 for (k = fdof/fcomp-1; k >= 0; --k) { 4113 for (c = 0; c < fcomp; ++c) { 4114 if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) { 4115 indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1); 4116 ++cind; 4117 } else { 4118 indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c-cind; 4119 } 4120 } 4121 } 4122 } 4123 } 4124 foff += fdof - cfdof; 4125 foffs[f] += fdof; 4126 } 4127 PetscFunctionReturn(0); 4128 } 4129 4130 #undef __FUNCT__ 4131 #define __FUNCT__ "DMPlexMatSetClosure" 4132 /*@C 4133 DMPlexMatSetClosure - Set an array of the values on the closure of 'point' 4134 4135 Not collective 4136 4137 Input Parameters: 4138 + dm - The DM 4139 . section - The section describing the layout in v, or NULL to use the default section 4140 . globalSection - The section describing the layout in v, or NULL to use the default global section 4141 . A - The matrix 4142 . point - The sieve point in the DM 4143 . values - The array of values 4144 - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions 4145 4146 Fortran Notes: 4147 This routine is only available in Fortran 90, and you must include petsc.h90 in your code. 4148 4149 Level: intermediate 4150 4151 .seealso DMPlexVecGetClosure(), DMPlexVecSetClosure() 4152 @*/ 4153 PetscErrorCode DMPlexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode) 4154 { 4155 DM_Plex *mesh = (DM_Plex*) dm->data; 4156 PetscSection clSection; 4157 IS clPoints; 4158 PetscInt *points = NULL; 4159 const PetscInt *clp; 4160 PetscInt *indices; 4161 PetscInt offsets[32]; 4162 PetscInt numFields, numPoints, numIndices, dof, off, globalOff, pStart, pEnd, p, q, f; 4163 PetscErrorCode ierr; 4164 4165 PetscFunctionBegin; 4166 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4167 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 4168 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 4169 if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);} 4170 PetscValidHeaderSpecific(globalSection, PETSC_SECTION_CLASSID, 3); 4171 PetscValidHeaderSpecific(A, MAT_CLASSID, 4); 4172 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 4173 if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields); 4174 ierr = PetscMemzero(offsets, 32 * sizeof(PetscInt));CHKERRQ(ierr); 4175 ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr); 4176 if (!clPoints) { 4177 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 4178 /* Compress out points not in the section */ 4179 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 4180 for (p = 0, q = 0; p < numPoints*2; p += 2) { 4181 if ((points[p] >= pStart) && (points[p] < pEnd)) { 4182 points[q*2] = points[p]; 4183 points[q*2+1] = points[p+1]; 4184 ++q; 4185 } 4186 } 4187 numPoints = q; 4188 } else { 4189 PetscInt dof, off; 4190 4191 ierr = PetscSectionGetDof(clSection, point, &dof);CHKERRQ(ierr); 4192 numPoints = dof/2; 4193 ierr = PetscSectionGetOffset(clSection, point, &off);CHKERRQ(ierr); 4194 ierr = ISGetIndices(clPoints, &clp);CHKERRQ(ierr); 4195 points = (PetscInt *) &clp[off]; 4196 } 4197 for (p = 0, numIndices = 0; p < numPoints*2; p += 2) { 4198 PetscInt fdof; 4199 4200 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 4201 for (f = 0; f < numFields; ++f) { 4202 ierr = PetscSectionGetFieldDof(section, points[p], f, &fdof);CHKERRQ(ierr); 4203 offsets[f+1] += fdof; 4204 } 4205 numIndices += dof; 4206 } 4207 for (f = 1; f < numFields; ++f) offsets[f+1] += offsets[f]; 4208 4209 if (numFields && offsets[numFields] != numIndices) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", offsets[numFields], numIndices); 4210 ierr = DMGetWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr); 4211 if (numFields) { 4212 for (p = 0; p < numPoints*2; p += 2) { 4213 PetscInt o = points[p+1]; 4214 ierr = PetscSectionGetOffset(globalSection, points[p], &globalOff);CHKERRQ(ierr); 4215 indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices); 4216 } 4217 } else { 4218 for (p = 0, off = 0; p < numPoints*2; p += 2) { 4219 PetscInt o = points[p+1]; 4220 ierr = PetscSectionGetOffset(globalSection, points[p], &globalOff);CHKERRQ(ierr); 4221 indicesPoint_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, indices); 4222 } 4223 } 4224 if (mesh->printSetValues) {ierr = DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr);} 4225 ierr = MatSetValues(A, numIndices, indices, numIndices, indices, values, mode); 4226 if (ierr) { 4227 PetscMPIInt rank; 4228 PetscErrorCode ierr2; 4229 4230 ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2); 4231 ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2); 4232 ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numIndices, indices, 0, NULL, values);CHKERRQ(ierr2); 4233 ierr2 = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr2); 4234 CHKERRQ(ierr); 4235 } 4236 if (!clPoints) { 4237 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 4238 } else { 4239 ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr); 4240 } 4241 ierr = DMRestoreWorkArray(dm, numIndices, PETSC_INT, &indices);CHKERRQ(ierr); 4242 PetscFunctionReturn(0); 4243 } 4244 4245 #undef __FUNCT__ 4246 #define __FUNCT__ "DMPlexMatSetClosureRefined" 4247 PetscErrorCode DMPlexMatSetClosureRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, Mat A, PetscInt point, const PetscScalar values[], InsertMode mode) 4248 { 4249 DM_Plex *mesh = (DM_Plex*) dmf->data; 4250 PetscInt *fpoints = NULL, *ftotpoints = NULL; 4251 PetscInt *cpoints = NULL; 4252 PetscInt *findices, *cindices; 4253 PetscInt foffsets[32], coffsets[32]; 4254 CellRefiner cellRefiner; 4255 PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f; 4256 PetscErrorCode ierr; 4257 4258 PetscFunctionBegin; 4259 PetscValidHeaderSpecific(dmf, DM_CLASSID, 1); 4260 PetscValidHeaderSpecific(dmc, DM_CLASSID, 4); 4261 if (!fsection) {ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);} 4262 PetscValidHeaderSpecific(fsection, PETSC_SECTION_CLASSID, 2); 4263 if (!csection) {ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);} 4264 PetscValidHeaderSpecific(csection, PETSC_SECTION_CLASSID, 5); 4265 if (!globalFSection) {ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr);} 4266 PetscValidHeaderSpecific(globalFSection, PETSC_SECTION_CLASSID, 3); 4267 if (!globalCSection) {ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr);} 4268 PetscValidHeaderSpecific(globalCSection, PETSC_SECTION_CLASSID, 6); 4269 PetscValidHeaderSpecific(A, MAT_CLASSID, 7); 4270 ierr = PetscSectionGetNumFields(fsection, &numFields);CHKERRQ(ierr); 4271 if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields); 4272 ierr = PetscMemzero(foffsets, 32 * sizeof(PetscInt));CHKERRQ(ierr); 4273 ierr = PetscMemzero(coffsets, 32 * sizeof(PetscInt));CHKERRQ(ierr); 4274 /* Column indices */ 4275 ierr = DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);CHKERRQ(ierr); 4276 maxFPoints = numCPoints; 4277 /* Compress out points not in the section */ 4278 /* TODO: Squeeze out points with 0 dof as well */ 4279 ierr = PetscSectionGetChart(csection, &pStart, &pEnd);CHKERRQ(ierr); 4280 for (p = 0, q = 0; p < numCPoints*2; p += 2) { 4281 if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) { 4282 cpoints[q*2] = cpoints[p]; 4283 cpoints[q*2+1] = cpoints[p+1]; 4284 ++q; 4285 } 4286 } 4287 numCPoints = q; 4288 for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) { 4289 PetscInt fdof; 4290 4291 ierr = PetscSectionGetDof(csection, cpoints[p], &dof);CHKERRQ(ierr); 4292 if (!dof) continue; 4293 for (f = 0; f < numFields; ++f) { 4294 ierr = PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);CHKERRQ(ierr); 4295 coffsets[f+1] += fdof; 4296 } 4297 numCIndices += dof; 4298 } 4299 for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f]; 4300 /* Row indices */ 4301 ierr = DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);CHKERRQ(ierr); 4302 ierr = CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);CHKERRQ(ierr); 4303 ierr = DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);CHKERRQ(ierr); 4304 for (r = 0, q = 0; r < numSubcells; ++r) { 4305 /* TODO Map from coarse to fine cells */ 4306 ierr = DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);CHKERRQ(ierr); 4307 /* Compress out points not in the section */ 4308 ierr = PetscSectionGetChart(fsection, &pStart, &pEnd);CHKERRQ(ierr); 4309 for (p = 0; p < numFPoints*2; p += 2) { 4310 if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) { 4311 ierr = PetscSectionGetDof(fsection, fpoints[p], &dof);CHKERRQ(ierr); 4312 if (!dof) continue; 4313 for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break; 4314 if (s < q) continue; 4315 ftotpoints[q*2] = fpoints[p]; 4316 ftotpoints[q*2+1] = fpoints[p+1]; 4317 ++q; 4318 } 4319 } 4320 ierr = DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);CHKERRQ(ierr); 4321 } 4322 numFPoints = q; 4323 for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) { 4324 PetscInt fdof; 4325 4326 ierr = PetscSectionGetDof(fsection, ftotpoints[p], &dof);CHKERRQ(ierr); 4327 if (!dof) continue; 4328 for (f = 0; f < numFields; ++f) { 4329 ierr = PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);CHKERRQ(ierr); 4330 foffsets[f+1] += fdof; 4331 } 4332 numFIndices += dof; 4333 } 4334 for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f]; 4335 4336 if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices); 4337 if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices); 4338 ierr = DMGetWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr); 4339 ierr = DMGetWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr); 4340 if (numFields) { 4341 for (p = 0; p < numFPoints*2; p += 2) { 4342 PetscInt o = ftotpoints[p+1]; 4343 ierr = PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);CHKERRQ(ierr); 4344 indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices); 4345 } 4346 for (p = 0; p < numCPoints*2; p += 2) { 4347 PetscInt o = cpoints[p+1]; 4348 ierr = PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);CHKERRQ(ierr); 4349 indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices); 4350 } 4351 } else { 4352 for (p = 0, off = 0; p < numFPoints*2; p += 2) { 4353 PetscInt o = ftotpoints[p+1]; 4354 ierr = PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);CHKERRQ(ierr); 4355 indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices); 4356 } 4357 for (p = 0, off = 0; p < numCPoints*2; p += 2) { 4358 PetscInt o = cpoints[p+1]; 4359 ierr = PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);CHKERRQ(ierr); 4360 indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices); 4361 } 4362 } 4363 if (mesh->printSetValues) {ierr = DMPlexPrintMatSetValues(PETSC_VIEWER_STDOUT_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr);} 4364 ierr = MatSetValues(A, numFIndices, findices, numCIndices, cindices, values, mode); 4365 if (ierr) { 4366 PetscMPIInt rank; 4367 PetscErrorCode ierr2; 4368 4369 ierr2 = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank);CHKERRQ(ierr2); 4370 ierr2 = (*PetscErrorPrintf)("[%D]ERROR in DMPlexMatSetClosure\n", rank);CHKERRQ(ierr2); 4371 ierr2 = DMPlexPrintMatSetValues(PETSC_VIEWER_STDERR_SELF, A, point, numFIndices, findices, numCIndices, cindices, values);CHKERRQ(ierr2); 4372 ierr2 = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr2); 4373 ierr2 = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr2); 4374 CHKERRQ(ierr); 4375 } 4376 ierr = DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);CHKERRQ(ierr); 4377 ierr = DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);CHKERRQ(ierr); 4378 ierr = DMRestoreWorkArray(dmf, numFIndices, PETSC_INT, &findices);CHKERRQ(ierr); 4379 ierr = DMRestoreWorkArray(dmc, numCIndices, PETSC_INT, &cindices);CHKERRQ(ierr); 4380 PetscFunctionReturn(0); 4381 } 4382 4383 #undef __FUNCT__ 4384 #define __FUNCT__ "DMPlexMatGetClosureIndicesRefined" 4385 PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[]) 4386 { 4387 PetscInt *fpoints = NULL, *ftotpoints = NULL; 4388 PetscInt *cpoints = NULL; 4389 PetscInt foffsets[32], coffsets[32]; 4390 CellRefiner cellRefiner; 4391 PetscInt numFields, numSubcells, maxFPoints, numFPoints, numCPoints, numFIndices, numCIndices, dof, off, globalOff, pStart, pEnd, p, q, r, s, f; 4392 PetscErrorCode ierr; 4393 4394 PetscFunctionBegin; 4395 PetscValidHeaderSpecific(dmf, DM_CLASSID, 1); 4396 PetscValidHeaderSpecific(dmc, DM_CLASSID, 4); 4397 if (!fsection) {ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);} 4398 PetscValidHeaderSpecific(fsection, PETSC_SECTION_CLASSID, 2); 4399 if (!csection) {ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);} 4400 PetscValidHeaderSpecific(csection, PETSC_SECTION_CLASSID, 5); 4401 if (!globalFSection) {ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr);} 4402 PetscValidHeaderSpecific(globalFSection, PETSC_SECTION_CLASSID, 3); 4403 if (!globalCSection) {ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr);} 4404 PetscValidHeaderSpecific(globalCSection, PETSC_SECTION_CLASSID, 6); 4405 ierr = PetscSectionGetNumFields(fsection, &numFields);CHKERRQ(ierr); 4406 if (numFields > 31) SETERRQ1(PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 31", numFields); 4407 ierr = PetscMemzero(foffsets, 32 * sizeof(PetscInt));CHKERRQ(ierr); 4408 ierr = PetscMemzero(coffsets, 32 * sizeof(PetscInt));CHKERRQ(ierr); 4409 /* Column indices */ 4410 ierr = DMPlexGetTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);CHKERRQ(ierr); 4411 maxFPoints = numCPoints; 4412 /* Compress out points not in the section */ 4413 /* TODO: Squeeze out points with 0 dof as well */ 4414 ierr = PetscSectionGetChart(csection, &pStart, &pEnd);CHKERRQ(ierr); 4415 for (p = 0, q = 0; p < numCPoints*2; p += 2) { 4416 if ((cpoints[p] >= pStart) && (cpoints[p] < pEnd)) { 4417 cpoints[q*2] = cpoints[p]; 4418 cpoints[q*2+1] = cpoints[p+1]; 4419 ++q; 4420 } 4421 } 4422 numCPoints = q; 4423 for (p = 0, numCIndices = 0; p < numCPoints*2; p += 2) { 4424 PetscInt fdof; 4425 4426 ierr = PetscSectionGetDof(csection, cpoints[p], &dof);CHKERRQ(ierr); 4427 if (!dof) continue; 4428 for (f = 0; f < numFields; ++f) { 4429 ierr = PetscSectionGetFieldDof(csection, cpoints[p], f, &fdof);CHKERRQ(ierr); 4430 coffsets[f+1] += fdof; 4431 } 4432 numCIndices += dof; 4433 } 4434 for (f = 1; f < numFields; ++f) coffsets[f+1] += coffsets[f]; 4435 /* Row indices */ 4436 ierr = DMPlexGetCellRefiner_Internal(dmc, &cellRefiner);CHKERRQ(ierr); 4437 ierr = CellRefinerGetAffineTransforms_Internal(cellRefiner, &numSubcells, NULL, NULL, NULL);CHKERRQ(ierr); 4438 ierr = DMGetWorkArray(dmf, maxFPoints*2*numSubcells, PETSC_INT, &ftotpoints);CHKERRQ(ierr); 4439 for (r = 0, q = 0; r < numSubcells; ++r) { 4440 /* TODO Map from coarse to fine cells */ 4441 ierr = DMPlexGetTransitiveClosure(dmf, point*numSubcells + r, PETSC_TRUE, &numFPoints, &fpoints);CHKERRQ(ierr); 4442 /* Compress out points not in the section */ 4443 ierr = PetscSectionGetChart(fsection, &pStart, &pEnd);CHKERRQ(ierr); 4444 for (p = 0; p < numFPoints*2; p += 2) { 4445 if ((fpoints[p] >= pStart) && (fpoints[p] < pEnd)) { 4446 ierr = PetscSectionGetDof(fsection, fpoints[p], &dof);CHKERRQ(ierr); 4447 if (!dof) continue; 4448 for (s = 0; s < q; ++s) if (fpoints[p] == ftotpoints[s*2]) break; 4449 if (s < q) continue; 4450 ftotpoints[q*2] = fpoints[p]; 4451 ftotpoints[q*2+1] = fpoints[p+1]; 4452 ++q; 4453 } 4454 } 4455 ierr = DMPlexRestoreTransitiveClosure(dmf, point, PETSC_TRUE, &numFPoints, &fpoints);CHKERRQ(ierr); 4456 } 4457 numFPoints = q; 4458 for (p = 0, numFIndices = 0; p < numFPoints*2; p += 2) { 4459 PetscInt fdof; 4460 4461 ierr = PetscSectionGetDof(fsection, ftotpoints[p], &dof);CHKERRQ(ierr); 4462 if (!dof) continue; 4463 for (f = 0; f < numFields; ++f) { 4464 ierr = PetscSectionGetFieldDof(fsection, ftotpoints[p], f, &fdof);CHKERRQ(ierr); 4465 foffsets[f+1] += fdof; 4466 } 4467 numFIndices += dof; 4468 } 4469 for (f = 1; f < numFields; ++f) foffsets[f+1] += foffsets[f]; 4470 4471 if (numFields && foffsets[numFields] != numFIndices) SETERRQ2(PetscObjectComm((PetscObject)dmf), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", foffsets[numFields], numFIndices); 4472 if (numFields && coffsets[numFields] != numCIndices) SETERRQ2(PetscObjectComm((PetscObject)dmc), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", coffsets[numFields], numCIndices); 4473 if (numFields) { 4474 for (p = 0; p < numFPoints*2; p += 2) { 4475 PetscInt o = ftotpoints[p+1]; 4476 ierr = PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);CHKERRQ(ierr); 4477 indicesPointFields_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, foffsets, PETSC_FALSE, o, findices); 4478 } 4479 for (p = 0; p < numCPoints*2; p += 2) { 4480 PetscInt o = cpoints[p+1]; 4481 ierr = PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);CHKERRQ(ierr); 4482 indicesPointFields_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, coffsets, PETSC_FALSE, o, cindices); 4483 } 4484 } else { 4485 for (p = 0, off = 0; p < numFPoints*2; p += 2) { 4486 PetscInt o = ftotpoints[p+1]; 4487 ierr = PetscSectionGetOffset(globalFSection, ftotpoints[p], &globalOff);CHKERRQ(ierr); 4488 indicesPoint_private(fsection, ftotpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, findices); 4489 } 4490 for (p = 0, off = 0; p < numCPoints*2; p += 2) { 4491 PetscInt o = cpoints[p+1]; 4492 ierr = PetscSectionGetOffset(globalCSection, cpoints[p], &globalOff);CHKERRQ(ierr); 4493 indicesPoint_private(csection, cpoints[p], globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE, o, cindices); 4494 } 4495 } 4496 ierr = DMRestoreWorkArray(dmf, numCPoints*2*4, PETSC_INT, &ftotpoints);CHKERRQ(ierr); 4497 ierr = DMPlexRestoreTransitiveClosure(dmc, point, PETSC_TRUE, &numCPoints, &cpoints);CHKERRQ(ierr); 4498 PetscFunctionReturn(0); 4499 } 4500 4501 #undef __FUNCT__ 4502 #define __FUNCT__ "DMPlexGetHybridBounds" 4503 /*@ 4504 DMPlexGetHybridBounds - Get the first mesh point of each dimension which is a hybrid 4505 4506 Input Parameter: 4507 . dm - The DMPlex object 4508 4509 Output Parameters: 4510 + cMax - The first hybrid cell 4511 . fMax - The first hybrid face 4512 . eMax - The first hybrid edge 4513 - vMax - The first hybrid vertex 4514 4515 Level: developer 4516 4517 .seealso DMPlexCreateHybridMesh(), DMPlexSetHybridBounds() 4518 @*/ 4519 PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax) 4520 { 4521 DM_Plex *mesh = (DM_Plex*) dm->data; 4522 PetscInt dim; 4523 PetscErrorCode ierr; 4524 4525 PetscFunctionBegin; 4526 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4527 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 4528 if (cMax) *cMax = mesh->hybridPointMax[dim]; 4529 if (fMax) *fMax = mesh->hybridPointMax[dim-1]; 4530 if (eMax) *eMax = mesh->hybridPointMax[1]; 4531 if (vMax) *vMax = mesh->hybridPointMax[0]; 4532 PetscFunctionReturn(0); 4533 } 4534 4535 #undef __FUNCT__ 4536 #define __FUNCT__ "DMPlexSetHybridBounds" 4537 /*@ 4538 DMPlexSetHybridBounds - Set the first mesh point of each dimension which is a hybrid 4539 4540 Input Parameters: 4541 . dm - The DMPlex object 4542 . cMax - The first hybrid cell 4543 . fMax - The first hybrid face 4544 . eMax - The first hybrid edge 4545 - vMax - The first hybrid vertex 4546 4547 Level: developer 4548 4549 .seealso DMPlexCreateHybridMesh(), DMPlexGetHybridBounds() 4550 @*/ 4551 PetscErrorCode DMPlexSetHybridBounds(DM dm, PetscInt cMax, PetscInt fMax, PetscInt eMax, PetscInt vMax) 4552 { 4553 DM_Plex *mesh = (DM_Plex*) dm->data; 4554 PetscInt dim; 4555 PetscErrorCode ierr; 4556 4557 PetscFunctionBegin; 4558 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4559 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 4560 if (cMax >= 0) mesh->hybridPointMax[dim] = cMax; 4561 if (fMax >= 0) mesh->hybridPointMax[dim-1] = fMax; 4562 if (eMax >= 0) mesh->hybridPointMax[1] = eMax; 4563 if (vMax >= 0) mesh->hybridPointMax[0] = vMax; 4564 PetscFunctionReturn(0); 4565 } 4566 4567 #undef __FUNCT__ 4568 #define __FUNCT__ "DMPlexGetVTKCellHeight" 4569 PetscErrorCode DMPlexGetVTKCellHeight(DM dm, PetscInt *cellHeight) 4570 { 4571 DM_Plex *mesh = (DM_Plex*) dm->data; 4572 4573 PetscFunctionBegin; 4574 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4575 PetscValidPointer(cellHeight, 2); 4576 *cellHeight = mesh->vtkCellHeight; 4577 PetscFunctionReturn(0); 4578 } 4579 4580 #undef __FUNCT__ 4581 #define __FUNCT__ "DMPlexSetVTKCellHeight" 4582 PetscErrorCode DMPlexSetVTKCellHeight(DM dm, PetscInt cellHeight) 4583 { 4584 DM_Plex *mesh = (DM_Plex*) dm->data; 4585 4586 PetscFunctionBegin; 4587 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4588 mesh->vtkCellHeight = cellHeight; 4589 PetscFunctionReturn(0); 4590 } 4591 4592 #undef __FUNCT__ 4593 #define __FUNCT__ "DMPlexCreateNumbering_Private" 4594 /* We can easily have a form that takes an IS instead */ 4595 static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering) 4596 { 4597 PetscSection section, globalSection; 4598 PetscInt *numbers, p; 4599 PetscErrorCode ierr; 4600 4601 PetscFunctionBegin; 4602 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 4603 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 4604 for (p = pStart; p < pEnd; ++p) { 4605 ierr = PetscSectionSetDof(section, p, 1);CHKERRQ(ierr); 4606 } 4607 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 4608 ierr = PetscSectionCreateGlobalSection(section, sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr); 4609 ierr = PetscMalloc1((pEnd - pStart), &numbers);CHKERRQ(ierr); 4610 for (p = pStart; p < pEnd; ++p) { 4611 ierr = PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);CHKERRQ(ierr); 4612 if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift; 4613 else numbers[p-pStart] += shift; 4614 } 4615 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);CHKERRQ(ierr); 4616 if (globalSize) { 4617 PetscLayout layout; 4618 ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);CHKERRQ(ierr); 4619 ierr = PetscLayoutGetSize(layout, globalSize);CHKERRQ(ierr); 4620 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 4621 } 4622 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 4623 ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr); 4624 PetscFunctionReturn(0); 4625 } 4626 4627 #undef __FUNCT__ 4628 #define __FUNCT__ "DMPlexGetCellNumbering" 4629 PetscErrorCode DMPlexGetCellNumbering(DM dm, IS *globalCellNumbers) 4630 { 4631 DM_Plex *mesh = (DM_Plex*) dm->data; 4632 PetscInt cellHeight, cStart, cEnd, cMax; 4633 PetscErrorCode ierr; 4634 4635 PetscFunctionBegin; 4636 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4637 if (!mesh->globalCellNumbers) { 4638 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 4639 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 4640 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 4641 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 4642 ierr = DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);CHKERRQ(ierr); 4643 } 4644 *globalCellNumbers = mesh->globalCellNumbers; 4645 PetscFunctionReturn(0); 4646 } 4647 4648 #undef __FUNCT__ 4649 #define __FUNCT__ "DMPlexGetVertexNumbering" 4650 PetscErrorCode DMPlexGetVertexNumbering(DM dm, IS *globalVertexNumbers) 4651 { 4652 DM_Plex *mesh = (DM_Plex*) dm->data; 4653 PetscInt vStart, vEnd, vMax; 4654 PetscErrorCode ierr; 4655 4656 PetscFunctionBegin; 4657 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4658 if (!mesh->globalVertexNumbers) { 4659 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 4660 ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr); 4661 if (vMax >= 0) vEnd = PetscMin(vEnd, vMax); 4662 ierr = DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);CHKERRQ(ierr); 4663 } 4664 *globalVertexNumbers = mesh->globalVertexNumbers; 4665 PetscFunctionReturn(0); 4666 } 4667 4668 #undef __FUNCT__ 4669 #define __FUNCT__ "DMPlexCreatePointNumbering" 4670 PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers) 4671 { 4672 IS nums[4]; 4673 PetscInt depths[4]; 4674 PetscInt depth, d, shift = 0; 4675 PetscErrorCode ierr; 4676 4677 PetscFunctionBegin; 4678 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4679 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 4680 depths[0] = depth; depths[1] = 0; 4681 for (d = 2; d <= depth; ++d) depths[d] = depth-d+1; 4682 for (d = 0; d <= depth; ++d) { 4683 PetscInt pStart, pEnd, gsize; 4684 4685 ierr = DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);CHKERRQ(ierr); 4686 ierr = DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);CHKERRQ(ierr); 4687 shift += gsize; 4688 } 4689 ierr = ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers); 4690 for (d = 0; d <= depth; ++d) {ierr = ISDestroy(&nums[d]);CHKERRQ(ierr);} 4691 PetscFunctionReturn(0); 4692 } 4693 4694 4695 #undef __FUNCT__ 4696 #define __FUNCT__ "PetscSectionCreateGlobalSectionLabel" 4697 /*@C 4698 PetscSectionCreateGlobalSectionLabel - Create a section describing the global field layout using 4699 the local section and an SF describing the section point overlap. 4700 4701 Input Parameters: 4702 + s - The PetscSection for the local field layout 4703 . sf - The SF describing parallel layout of the section points 4704 . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs 4705 . label - The label specifying the points 4706 - labelValue - The label stratum specifying the points 4707 4708 Output Parameter: 4709 . gsection - The PetscSection for the global field layout 4710 4711 Note: This gives negative sizes and offsets to points not owned by this process 4712 4713 Level: developer 4714 4715 .seealso: PetscSectionCreate() 4716 @*/ 4717 PetscErrorCode PetscSectionCreateGlobalSectionLabel(PetscSection s, PetscSF sf, PetscBool includeConstraints, DMLabel label, PetscInt labelValue, PetscSection *gsection) 4718 { 4719 PetscInt *neg = NULL, *tmpOff = NULL; 4720 PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots; 4721 PetscErrorCode ierr; 4722 4723 PetscFunctionBegin; 4724 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);CHKERRQ(ierr); 4725 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 4726 ierr = PetscSectionSetChart(*gsection, pStart, pEnd);CHKERRQ(ierr); 4727 ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 4728 if (nroots >= 0) { 4729 if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart); 4730 ierr = PetscCalloc1(nroots, &neg);CHKERRQ(ierr); 4731 if (nroots > pEnd-pStart) { 4732 ierr = PetscCalloc1(nroots, &tmpOff);CHKERRQ(ierr); 4733 } else { 4734 tmpOff = &(*gsection)->atlasDof[-pStart]; 4735 } 4736 } 4737 /* Mark ghost points with negative dof */ 4738 for (p = pStart; p < pEnd; ++p) { 4739 PetscInt value; 4740 4741 ierr = DMLabelGetValue(label, p, &value);CHKERRQ(ierr); 4742 if (value != labelValue) continue; 4743 ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr); 4744 ierr = PetscSectionSetDof(*gsection, p, dof);CHKERRQ(ierr); 4745 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 4746 if (!includeConstraints && cdof > 0) {ierr = PetscSectionSetConstraintDof(*gsection, p, cdof);CHKERRQ(ierr);} 4747 if (neg) neg[p] = -(dof+1); 4748 } 4749 ierr = PetscSectionSetUpBC(*gsection);CHKERRQ(ierr); 4750 if (nroots >= 0) { 4751 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr); 4752 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr); 4753 if (nroots > pEnd-pStart) { 4754 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];} 4755 } 4756 } 4757 /* Calculate new sizes, get proccess offset, and calculate point offsets */ 4758 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 4759 cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0; 4760 (*gsection)->atlasOff[p] = off; 4761 off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0; 4762 } 4763 ierr = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));CHKERRQ(ierr); 4764 globalOff -= off; 4765 for (p = 0, off = 0; p < pEnd-pStart; ++p) { 4766 (*gsection)->atlasOff[p] += globalOff; 4767 if (neg) neg[p] = -((*gsection)->atlasOff[p]+1); 4768 } 4769 /* Put in negative offsets for ghost points */ 4770 if (nroots >= 0) { 4771 ierr = PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr); 4772 ierr = PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);CHKERRQ(ierr); 4773 if (nroots > pEnd-pStart) { 4774 for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];} 4775 } 4776 } 4777 if (nroots >= 0 && nroots > pEnd-pStart) {ierr = PetscFree(tmpOff);CHKERRQ(ierr);} 4778 ierr = PetscFree(neg);CHKERRQ(ierr); 4779 PetscFunctionReturn(0); 4780 } 4781 4782 #undef __FUNCT__ 4783 #define __FUNCT__ "DMPlexCheckSymmetry" 4784 /*@ 4785 DMPlexCheckSymmetry - Check that the adjacency information in the mesh is symmetric. 4786 4787 Input Parameters: 4788 + dm - The DMPlex object 4789 4790 Note: This is a useful diagnostic when creating meshes programmatically. 4791 4792 Level: developer 4793 4794 .seealso: DMCreate(), DMCheckSkeleton(), DMCheckFaces() 4795 @*/ 4796 PetscErrorCode DMPlexCheckSymmetry(DM dm) 4797 { 4798 PetscSection coneSection, supportSection; 4799 const PetscInt *cone, *support; 4800 PetscInt coneSize, c, supportSize, s; 4801 PetscInt pStart, pEnd, p, csize, ssize; 4802 PetscErrorCode ierr; 4803 4804 PetscFunctionBegin; 4805 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4806 ierr = DMPlexGetConeSection(dm, &coneSection);CHKERRQ(ierr); 4807 ierr = DMPlexGetSupportSection(dm, &supportSection);CHKERRQ(ierr); 4808 /* Check that point p is found in the support of its cone points, and vice versa */ 4809 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 4810 for (p = pStart; p < pEnd; ++p) { 4811 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 4812 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 4813 for (c = 0; c < coneSize; ++c) { 4814 PetscBool dup = PETSC_FALSE; 4815 PetscInt d; 4816 for (d = c-1; d >= 0; --d) { 4817 if (cone[c] == cone[d]) {dup = PETSC_TRUE; break;} 4818 } 4819 ierr = DMPlexGetSupportSize(dm, cone[c], &supportSize);CHKERRQ(ierr); 4820 ierr = DMPlexGetSupport(dm, cone[c], &support);CHKERRQ(ierr); 4821 for (s = 0; s < supportSize; ++s) { 4822 if (support[s] == p) break; 4823 } 4824 if ((s >= supportSize) || (dup && (support[s+1] != p))) { 4825 ierr = PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", p); 4826 for (s = 0; s < coneSize; ++s) { 4827 ierr = PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[s]); 4828 } 4829 ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); 4830 ierr = PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", cone[c]); 4831 for (s = 0; s < supportSize; ++s) { 4832 ierr = PetscPrintf(PETSC_COMM_SELF, "%d, ", support[s]); 4833 } 4834 ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); 4835 if (dup) { 4836 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not repeatedly found in support of repeated cone point %d", p, cone[c]); 4837 } else { 4838 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in support of cone point %d", p, cone[c]); 4839 } 4840 } 4841 } 4842 ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr); 4843 ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr); 4844 for (s = 0; s < supportSize; ++s) { 4845 ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr); 4846 ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr); 4847 for (c = 0; c < coneSize; ++c) { 4848 if (cone[c] == p) break; 4849 } 4850 if (c >= coneSize) { 4851 ierr = PetscPrintf(PETSC_COMM_SELF, "p: %d support: ", p); 4852 for (c = 0; c < supportSize; ++c) { 4853 ierr = PetscPrintf(PETSC_COMM_SELF, "%d, ", support[c]); 4854 } 4855 ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); 4856 ierr = PetscPrintf(PETSC_COMM_SELF, "p: %d cone: ", support[s]); 4857 for (c = 0; c < coneSize; ++c) { 4858 ierr = PetscPrintf(PETSC_COMM_SELF, "%d, ", cone[c]); 4859 } 4860 ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); 4861 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point %d not found in cone of support point %d", p, support[s]); 4862 } 4863 } 4864 } 4865 ierr = PetscSectionGetStorageSize(coneSection, &csize);CHKERRQ(ierr); 4866 ierr = PetscSectionGetStorageSize(supportSection, &ssize);CHKERRQ(ierr); 4867 if (csize != ssize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total cone size %d != Total support size %d", csize, ssize); 4868 PetscFunctionReturn(0); 4869 } 4870 4871 #undef __FUNCT__ 4872 #define __FUNCT__ "DMPlexCheckSkeleton" 4873 /*@ 4874 DMPlexCheckSkeleton - Check that each cell has the correct number of vertices 4875 4876 Input Parameters: 4877 + dm - The DMPlex object 4878 . isSimplex - Are the cells simplices or tensor products 4879 - cellHeight - Normally 0 4880 4881 Note: This is a useful diagnostic when creating meshes programmatically. 4882 4883 Level: developer 4884 4885 .seealso: DMCreate(), DMCheckSymmetry(), DMCheckFaces() 4886 @*/ 4887 PetscErrorCode DMPlexCheckSkeleton(DM dm, PetscBool isSimplex, PetscInt cellHeight) 4888 { 4889 PetscInt dim, numCorners, numHybridCorners, vStart, vEnd, cStart, cEnd, cMax, c; 4890 PetscErrorCode ierr; 4891 4892 PetscFunctionBegin; 4893 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4894 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 4895 switch (dim) { 4896 case 1: numCorners = isSimplex ? 2 : 2; numHybridCorners = isSimplex ? 2 : 2; break; 4897 case 2: numCorners = isSimplex ? 3 : 4; numHybridCorners = isSimplex ? 4 : 4; break; 4898 case 3: numCorners = isSimplex ? 4 : 8; numHybridCorners = isSimplex ? 6 : 8; break; 4899 default: 4900 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle meshes of dimension %d", dim); 4901 } 4902 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 4903 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 4904 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 4905 cMax = cMax >= 0 ? cMax : cEnd; 4906 for (c = cStart; c < cMax; ++c) { 4907 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 4908 4909 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4910 for (cl = 0; cl < closureSize*2; cl += 2) { 4911 const PetscInt p = closure[cl]; 4912 if ((p >= vStart) && (p < vEnd)) ++coneSize; 4913 } 4914 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4915 if (coneSize != numCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d vertices != %d", c, coneSize, numCorners); 4916 } 4917 for (c = cMax; c < cEnd; ++c) { 4918 PetscInt *closure = NULL, closureSize, cl, coneSize = 0; 4919 4920 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4921 for (cl = 0; cl < closureSize*2; cl += 2) { 4922 const PetscInt p = closure[cl]; 4923 if ((p >= vStart) && (p < vEnd)) ++coneSize; 4924 } 4925 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4926 if (coneSize > numHybridCorners) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hybrid cell %d has %d vertices > %d", c, coneSize, numHybridCorners); 4927 } 4928 PetscFunctionReturn(0); 4929 } 4930 4931 #undef __FUNCT__ 4932 #define __FUNCT__ "DMPlexCheckFaces" 4933 /*@ 4934 DMPlexCheckFaces - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type 4935 4936 Input Parameters: 4937 + dm - The DMPlex object 4938 . isSimplex - Are the cells simplices or tensor products 4939 - cellHeight - Normally 0 4940 4941 Note: This is a useful diagnostic when creating meshes programmatically. 4942 4943 Level: developer 4944 4945 .seealso: DMCreate(), DMCheckSymmetry(), DMCheckSkeleton() 4946 @*/ 4947 PetscErrorCode DMPlexCheckFaces(DM dm, PetscBool isSimplex, PetscInt cellHeight) 4948 { 4949 PetscInt pMax[4]; 4950 PetscInt dim, vStart, vEnd, cStart, cEnd, c, h; 4951 PetscErrorCode ierr; 4952 4953 PetscFunctionBegin; 4954 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4955 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 4956 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 4957 ierr = DMPlexGetHybridBounds(dm, &pMax[dim], &pMax[dim-1], &pMax[1], &pMax[0]);CHKERRQ(ierr); 4958 for (h = cellHeight; h < dim; ++h) { 4959 ierr = DMPlexGetHeightStratum(dm, h, &cStart, &cEnd);CHKERRQ(ierr); 4960 for (c = cStart; c < cEnd; ++c) { 4961 const PetscInt *cone, *ornt, *faces; 4962 PetscInt numFaces, faceSize, coneSize,f; 4963 PetscInt *closure = NULL, closureSize, cl, numCorners = 0; 4964 4965 if (pMax[dim-h] >= 0 && c >= pMax[dim-h]) continue; 4966 ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr); 4967 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 4968 ierr = DMPlexGetConeOrientation(dm, c, &ornt);CHKERRQ(ierr); 4969 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4970 for (cl = 0; cl < closureSize*2; cl += 2) { 4971 const PetscInt p = closure[cl]; 4972 if ((p >= vStart) && (p < vEnd)) closure[numCorners++] = p; 4973 } 4974 ierr = DMPlexGetRawFaces_Internal(dm, dim-h, numCorners, closure, &numFaces, &faceSize, &faces);CHKERRQ(ierr); 4975 if (coneSize != numFaces) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %d has %d faces but should have %d", c, coneSize, numFaces); 4976 for (f = 0; f < numFaces; ++f) { 4977 PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v; 4978 4979 ierr = DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);CHKERRQ(ierr); 4980 for (cl = 0; cl < fclosureSize*2; cl += 2) { 4981 const PetscInt p = fclosure[cl]; 4982 if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p; 4983 } 4984 if (fnumCorners != faceSize) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d has %d vertices but should have %d", cone[f], f, c, fnumCorners, faceSize); 4985 for (v = 0; v < fnumCorners; ++v) { 4986 if (fclosure[v] != faces[f*faceSize+v]) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d (%d) of cell %d vertex %d, %d != %d", cone[f], f, c, v, fclosure[v], faces[f*faceSize+v]); 4987 } 4988 ierr = DMPlexRestoreTransitiveClosure(dm, cone[f], PETSC_TRUE, &fclosureSize, &fclosure);CHKERRQ(ierr); 4989 } 4990 ierr = DMPlexRestoreFaces_Internal(dm, dim, c, &numFaces, &faceSize, &faces);CHKERRQ(ierr); 4991 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 4992 } 4993 } 4994 PetscFunctionReturn(0); 4995 } 4996 4997 #undef __FUNCT__ 4998 #define __FUNCT__ "DMCreateInterpolation_Plex" 4999 /* Pointwise interpolation 5000 Just code FEM for now 5001 u^f = I u^c 5002 sum_k u^f_k phi^f_k = I sum_j u^c_j phi^c_j 5003 u^f_i = sum_j psi^f_i I phi^c_j u^c_j 5004 I_{ij} = psi^f_i phi^c_j 5005 */ 5006 PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling) 5007 { 5008 PetscSection gsc, gsf; 5009 PetscInt m, n; 5010 void *ctx; 5011 PetscErrorCode ierr; 5012 5013 PetscFunctionBegin; 5014 /* 5015 Loop over coarse cells 5016 Loop over coarse basis functions 5017 Loop over fine cells in coarse cell 5018 Loop over fine dual basis functions 5019 Evaluate coarse basis on fine dual basis quad points 5020 Sum 5021 Update local element matrix 5022 Accumulate to interpolation matrix 5023 5024 Can extend PetscFEIntegrateJacobian_Basic() to do a specialized cell loop 5025 */ 5026 ierr = DMGetDefaultGlobalSection(dmFine, &gsf);CHKERRQ(ierr); 5027 ierr = PetscSectionGetConstrainedStorageSize(gsf, &m);CHKERRQ(ierr); 5028 ierr = DMGetDefaultGlobalSection(dmCoarse, &gsc);CHKERRQ(ierr); 5029 ierr = PetscSectionGetConstrainedStorageSize(gsc, &n);CHKERRQ(ierr); 5030 /* We need to preallocate properly */ 5031 ierr = MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);CHKERRQ(ierr); 5032 ierr = MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 5033 ierr = MatSetType(*interpolation, dmCoarse->mattype);CHKERRQ(ierr); 5034 ierr = DMGetApplicationContext(dmFine, &ctx);CHKERRQ(ierr); 5035 ierr = DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);CHKERRQ(ierr); 5036 /* Use naive scaling */ 5037 ierr = DMCreateInterpolationScale(dmCoarse, dmFine, *interpolation, scaling);CHKERRQ(ierr); 5038 PetscFunctionReturn(0); 5039 } 5040 5041 #undef __FUNCT__ 5042 #define __FUNCT__ "DMCreateInjection_Plex" 5043 PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, VecScatter *ctx) 5044 { 5045 PetscErrorCode ierr; 5046 5047 PetscFunctionBegin; 5048 ierr = DMPlexComputeInjectorFEM(dmCoarse, dmFine, ctx, NULL);CHKERRQ(ierr); 5049 PetscFunctionReturn(0); 5050 } 5051 5052 #undef __FUNCT__ 5053 #define __FUNCT__ "DMCreateDefaultSection_Plex" 5054 /* Pointwise interpolation 5055 Just code FEM for now 5056 u^f = I u^c 5057 sum_k u^f_k phi^f_k = I sum_l u^c_l phi^c_l 5058 u^f_i = sum_l int psi^f_i I phi^c_l u^c_l 5059 I_{ij} = int psi^f_i phi^c_j 5060 */ 5061 PetscErrorCode DMCreateDefaultSection_Plex(DM dm) 5062 { 5063 PetscSection section; 5064 IS *bcPoints; 5065 PetscInt *bcFields, *numComp, *numDof; 5066 PetscInt depth, dim, numBd, numBC = 0, numFields, bd, bc, f; 5067 PetscErrorCode ierr; 5068 5069 PetscFunctionBegin; 5070 /* Handle boundary conditions */ 5071 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 5072 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 5073 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 5074 for (bd = 0; bd < numBd; ++bd) { 5075 PetscBool isEssential; 5076 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 5077 if (isEssential) ++numBC; 5078 } 5079 ierr = PetscMalloc2(numBC,&bcFields,numBC,&bcPoints);CHKERRQ(ierr); 5080 for (bd = 0, bc = 0; bd < numBd; ++bd) { 5081 const char *bdLabel; 5082 DMLabel label; 5083 const PetscInt *values; 5084 PetscInt bd2, field, numValues; 5085 PetscBool isEssential, duplicate = PETSC_FALSE; 5086 5087 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 5088 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 5089 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 5090 /* Only want to do this for FEM, and only once */ 5091 for (bd2 = 0; bd2 < bd; ++bd2) { 5092 const char *bdname; 5093 ierr = DMPlexGetBoundary(dm, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 5094 ierr = PetscStrcmp(bdname, bdLabel, &duplicate);CHKERRQ(ierr); 5095 if (duplicate) break; 5096 } 5097 if (!duplicate) { 5098 ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr); 5099 ierr = DMPlexLabelAddCells(dm, label);CHKERRQ(ierr); 5100 } 5101 /* Filter out cells, if you actually want to constraint cells you need to do things by hand right now */ 5102 if (isEssential) { 5103 IS tmp; 5104 PetscInt *newidx; 5105 const PetscInt *idx; 5106 PetscInt cStart, cEnd, n, p, newn = 0; 5107 5108 bcFields[bc] = field; 5109 ierr = DMPlexGetStratumIS(dm, bdLabel, values[0], &tmp);CHKERRQ(ierr); 5110 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 5111 ierr = ISGetLocalSize(tmp, &n);CHKERRQ(ierr); 5112 ierr = ISGetIndices(tmp, &idx);CHKERRQ(ierr); 5113 for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn; 5114 ierr = PetscMalloc1(newn,&newidx);CHKERRQ(ierr); 5115 newn = 0; 5116 for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p]; 5117 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);CHKERRQ(ierr); 5118 ierr = ISRestoreIndices(tmp, &idx);CHKERRQ(ierr); 5119 ierr = ISDestroy(&tmp);CHKERRQ(ierr); 5120 } 5121 } 5122 /* Handle discretization */ 5123 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 5124 ierr = PetscMalloc2(numFields,&numComp,numFields*(dim+1),&numDof);CHKERRQ(ierr); 5125 for (f = 0; f < numFields; ++f) { 5126 PetscFE fe; 5127 const PetscInt *numFieldDof; 5128 PetscInt d; 5129 5130 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 5131 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 5132 ierr = PetscFEGetNumDof(fe, &numFieldDof);CHKERRQ(ierr); 5133 for (d = 0; d < dim+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d]; 5134 } 5135 for (f = 0; f < numFields; ++f) { 5136 PetscInt d; 5137 for (d = 1; d < dim; ++d) { 5138 if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces."); 5139 } 5140 } 5141 ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, §ion);CHKERRQ(ierr); 5142 for (f = 0; f < numFields; ++f) { 5143 PetscFE fe; 5144 const char *name; 5145 5146 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 5147 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 5148 ierr = PetscSectionSetFieldName(section, f, name);CHKERRQ(ierr); 5149 } 5150 ierr = DMSetDefaultSection(dm, section);CHKERRQ(ierr); 5151 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 5152 for (bc = 0; bc < numBC; ++bc) {ierr = ISDestroy(&bcPoints[bc]);CHKERRQ(ierr);} 5153 ierr = PetscFree2(bcFields,bcPoints);CHKERRQ(ierr); 5154 ierr = PetscFree2(numComp,numDof);CHKERRQ(ierr); 5155 PetscFunctionReturn(0); 5156 } 5157 5158 #undef __FUNCT__ 5159 #define __FUNCT__ "DMPlexGetCoarseDM" 5160 /*@ 5161 DMPlexGetCoarseDM - Get the coarse mesh from which this was obtained by refinement 5162 5163 Input Parameter: 5164 . dm - The DMPlex object 5165 5166 Output Parameter: 5167 . cdm - The coarse DM 5168 5169 Level: intermediate 5170 5171 .seealso: DMPlexSetCoarseDM() 5172 @*/ 5173 PetscErrorCode DMPlexGetCoarseDM(DM dm, DM *cdm) 5174 { 5175 PetscFunctionBegin; 5176 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 5177 PetscValidPointer(cdm, 2); 5178 *cdm = ((DM_Plex *) dm->data)->coarseMesh; 5179 PetscFunctionReturn(0); 5180 } 5181 5182 #undef __FUNCT__ 5183 #define __FUNCT__ "DMPlexSetCoarseDM" 5184 /*@ 5185 DMPlexSetCoarseDM - Set the coarse mesh from which this was obtained by refinement 5186 5187 Input Parameters: 5188 + dm - The DMPlex object 5189 - cdm - The coarse DM 5190 5191 Level: intermediate 5192 5193 .seealso: DMPlexGetCoarseDM() 5194 @*/ 5195 PetscErrorCode DMPlexSetCoarseDM(DM dm, DM cdm) 5196 { 5197 DM_Plex *mesh; 5198 PetscErrorCode ierr; 5199 5200 PetscFunctionBegin; 5201 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 5202 if (cdm) PetscValidHeaderSpecific(cdm, DM_CLASSID, 2); 5203 mesh = (DM_Plex *) dm->data; 5204 ierr = DMDestroy(&mesh->coarseMesh);CHKERRQ(ierr); 5205 mesh->coarseMesh = cdm; 5206 ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr); 5207 PetscFunctionReturn(0); 5208 } 5209