1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petsc-private/petscfeimpl.h> 4 #include <petscfv.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMPlexGetScale" 8 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 9 { 10 DM_Plex *mesh = (DM_Plex*) dm->data; 11 12 PetscFunctionBegin; 13 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 14 PetscValidPointer(scale, 3); 15 *scale = mesh->scale[unit]; 16 PetscFunctionReturn(0); 17 } 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "DMPlexSetScale" 21 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 22 { 23 DM_Plex *mesh = (DM_Plex*) dm->data; 24 25 PetscFunctionBegin; 26 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 27 mesh->scale[unit] = scale; 28 PetscFunctionReturn(0); 29 } 30 31 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 32 { 33 switch (i) { 34 case 0: 35 switch (j) { 36 case 0: return 0; 37 case 1: 38 switch (k) { 39 case 0: return 0; 40 case 1: return 0; 41 case 2: return 1; 42 } 43 case 2: 44 switch (k) { 45 case 0: return 0; 46 case 1: return -1; 47 case 2: return 0; 48 } 49 } 50 case 1: 51 switch (j) { 52 case 0: 53 switch (k) { 54 case 0: return 0; 55 case 1: return 0; 56 case 2: return -1; 57 } 58 case 1: return 0; 59 case 2: 60 switch (k) { 61 case 0: return 1; 62 case 1: return 0; 63 case 2: return 0; 64 } 65 } 66 case 2: 67 switch (j) { 68 case 0: 69 switch (k) { 70 case 0: return 0; 71 case 1: return 1; 72 case 2: return 0; 73 } 74 case 1: 75 switch (k) { 76 case 0: return -1; 77 case 1: return 0; 78 case 2: return 0; 79 } 80 case 2: return 0; 81 } 82 } 83 return 0; 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "DMPlexCreateRigidBody" 88 /*@C 89 DMPlexCreateRigidBody - create rigid body modes from coordinates 90 91 Collective on DM 92 93 Input Arguments: 94 + dm - the DM 95 . section - the local section associated with the rigid field, or NULL for the default section 96 - globalSection - the global section associated with the rigid field, or NULL for the default section 97 98 Output Argument: 99 . sp - the null space 100 101 Note: This is necessary to take account of Dirichlet conditions on the displacements 102 103 Level: advanced 104 105 .seealso: MatNullSpaceCreate() 106 @*/ 107 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp) 108 { 109 MPI_Comm comm; 110 Vec coordinates, localMode, mode[6]; 111 PetscSection coordSection; 112 PetscScalar *coords; 113 PetscInt dim, vStart, vEnd, v, n, m, d, i, j; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 118 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 119 if (dim == 1) { 120 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 124 if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);} 125 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 126 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 127 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 128 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 129 m = (dim*(dim+1))/2; 130 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 131 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 132 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 133 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 134 /* Assume P1 */ 135 ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr); 136 for (d = 0; d < dim; ++d) { 137 PetscScalar values[3] = {0.0, 0.0, 0.0}; 138 139 values[d] = 1.0; 140 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 141 for (v = vStart; v < vEnd; ++v) { 142 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 143 } 144 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 145 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 146 } 147 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 148 for (d = dim; d < dim*(dim+1)/2; ++d) { 149 PetscInt i, j, k = dim > 2 ? d - dim : d; 150 151 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 152 for (v = vStart; v < vEnd; ++v) { 153 PetscScalar values[3] = {0.0, 0.0, 0.0}; 154 PetscInt off; 155 156 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 157 for (i = 0; i < dim; ++i) { 158 for (j = 0; j < dim; ++j) { 159 values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]); 160 } 161 } 162 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 163 } 164 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 165 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 166 } 167 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 168 ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr); 169 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 170 /* Orthonormalize system */ 171 for (i = dim; i < m; ++i) { 172 PetscScalar dots[6]; 173 174 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 175 for (j = 0; j < i; ++j) dots[j] *= -1.0; 176 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 177 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 178 } 179 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 180 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "DMPlexProjectFunctionLabelLocal" 186 PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 187 { 188 PetscDualSpace *sp; 189 PetscSection section; 190 PetscScalar *values; 191 PetscReal *v0, *J, detJ; 192 PetscBool *fieldActive; 193 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, d, v, i, comp; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 198 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 199 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 200 ierr = PetscMalloc3(numFields,&sp,dim,&v0,dim*dim,&J);CHKERRQ(ierr); 201 for (f = 0; f < numFields; ++f) { 202 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 203 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 204 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 205 totDim += spDim*numComp; 206 } 207 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 208 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 209 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 210 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 211 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 212 for (f = 0; f < numFields; ++f) fieldActive[f] = funcs[f] ? PETSC_TRUE : PETSC_FALSE; 213 for (i = 0; i < numIds; ++i) { 214 IS pointIS; 215 const PetscInt *points; 216 PetscInt n, p; 217 218 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 219 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 220 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 221 for (p = 0; p < n; ++p) { 222 const PetscInt point = points[p]; 223 PetscCellGeometry geom; 224 225 if ((point < cStart) || (point >= cEnd)) continue; 226 ierr = DMPlexComputeCellGeometry(dm, point, v0, J, NULL, &detJ);CHKERRQ(ierr); 227 geom.v0 = v0; 228 geom.J = J; 229 geom.detJ = &detJ; 230 for (f = 0, v = 0; f < numFields; ++f) { 231 void * const ctx = ctxs ? ctxs[f] : NULL; 232 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 233 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 234 for (d = 0; d < spDim; ++d) { 235 if (funcs[f]) { 236 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 237 } else { 238 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 239 } 240 v += numComp; 241 } 242 } 243 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 244 } 245 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 246 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 247 } 248 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 249 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 250 ierr = PetscFree3(sp,v0,J);CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNCT__ 255 #define __FUNCT__ "DMPlexProjectFunctionLocal" 256 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 257 { 258 PetscDualSpace *sp; 259 PetscSection section; 260 PetscScalar *values; 261 PetscReal *v0, *J, detJ; 262 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 267 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 268 ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr); 269 for (f = 0; f < numFields; ++f) { 270 PetscFE fe; 271 272 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 273 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 274 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 275 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 276 totDim += spDim*numComp; 277 } 278 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 279 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 280 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 281 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 282 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 283 ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr); 284 for (c = cStart; c < cEnd; ++c) { 285 PetscCellGeometry geom; 286 287 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 288 geom.v0 = v0; 289 geom.J = J; 290 geom.detJ = &detJ; 291 for (f = 0, v = 0; f < numFields; ++f) { 292 PetscFE fe; 293 void * const ctx = ctxs ? ctxs[f] : NULL; 294 295 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 296 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 297 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 298 for (d = 0; d < spDim; ++d) { 299 if (funcs[f]) { 300 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 301 } else { 302 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 303 } 304 v += numComp; 305 } 306 } 307 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 308 } 309 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 310 ierr = PetscFree2(v0,J);CHKERRQ(ierr); 311 ierr = PetscFree(sp);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "DMPlexProjectFunction" 317 /*@C 318 DMPlexProjectFunction - This projects the given function into the function space provided. 319 320 Input Parameters: 321 + dm - The DM 322 . funcs - The coordinate functions to evaluate, one per field 323 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 324 - mode - The insertion mode for values 325 326 Output Parameter: 327 . X - vector 328 329 Level: developer 330 331 .seealso: DMPlexComputeL2Diff() 332 @*/ 333 PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 334 { 335 Vec localX; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 340 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 341 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);CHKERRQ(ierr); 342 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 343 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 344 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "DMPlexProjectFieldLocal" 350 PetscErrorCode DMPlexProjectFieldLocal(DM dm, Vec localU, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec localX) 351 { 352 DM dmAux; 353 PetscProblem prob, probAux; 354 Vec A; 355 PetscSection section, sectionAux; 356 PetscScalar *values, *u, *gradU, *a, *gradA, *coefficients, *coefficientsAux; 357 PetscReal *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux; 358 PetscInt Nf, NfAux, cOffset, cOffsetAux, dim, spDim, totDim, numValues, cStart, cEnd, c, f, d, v, comp; 359 PetscErrorCode ierr; 360 361 PetscFunctionBegin; 362 ierr = DMGetProblem(dm, &prob);CHKERRQ(ierr); 363 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 364 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 365 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 366 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 367 ierr = PetscProblemGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 368 ierr = PetscProblemGetTabulation(prob, &basisField, &basisFieldDer);CHKERRQ(ierr); 369 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 370 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 371 if (dmAux) { 372 ierr = DMGetProblem(dmAux, &probAux);CHKERRQ(ierr); 373 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 374 ierr = PetscProblemGetTabulation(prob, &basisFieldAux, &basisFieldDerAux);CHKERRQ(ierr); 375 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 376 } 377 ierr = DMPlexInsertBoundaryValuesFEM(dm, localU);CHKERRQ(ierr); 378 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 379 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 380 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 381 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 382 for (c = cStart; c < cEnd; ++c) { 383 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 384 for (f = 0, v = 0; f < Nf; ++f) { 385 PetscFE fe; 386 PetscDualSpace sp; 387 PetscInt Ncf; 388 389 ierr = PetscProblemGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 390 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 391 ierr = PetscFEGetNumComponents(fe, &Ncf);CHKERRQ(ierr); 392 ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr); 393 for (d = 0; d < spDim; ++d) { 394 PetscQuadrature quad; 395 const PetscReal *points, *weights; 396 PetscInt numPoints, q; 397 398 if (funcs[f]) { 399 ierr = PetscDualSpaceGetFunctional(sp, d, &quad);CHKERRQ(ierr); 400 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 401 ierr = PetscFEGetTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr); 402 for (q = 0; q < numPoints; ++q) { 403 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 404 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], NULL, u, gradU, NULL);CHKERRQ(ierr); 405 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, gradA, NULL);CHKERRQ(ierr); 406 (*funcs[f])(u, gradU, a, gradA, x, &values[v]); 407 } 408 ierr = PetscFERestoreTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr); 409 } else { 410 for (comp = 0; comp < Ncf; ++comp) values[v+comp] = 0.0; 411 } 412 v += Ncf; 413 } 414 } 415 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 416 } 417 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 418 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 419 PetscFunctionReturn(0); 420 } 421 422 #undef __FUNCT__ 423 #define __FUNCT__ "DMPlexProjectField" 424 /*@C 425 DMPlexProjectField - This projects the given function of the fields into the function space provided. 426 427 Input Parameters: 428 + dm - The DM 429 . fe - The PetscFE associated with the field 430 . U - The input field vector 431 . funcs - The functions to evaluate, one per field 432 - mode - The insertion mode for values 433 434 Output Parameter: 435 . X - The output vector 436 437 Level: developer 438 439 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 440 @*/ 441 PetscErrorCode DMPlexProjectField(DM dm, PetscFE fe[], PetscFE feAux[], Vec U, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec X) 442 { 443 Vec localX, localU; 444 PetscErrorCode ierr; 445 446 PetscFunctionBegin; 447 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 448 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 449 ierr = DMGetLocalVector(dm, &localU);CHKERRQ(ierr); 450 ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr); 451 ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr); 452 ierr = DMPlexProjectFieldLocal(dm, localU, funcs, mode, localX);CHKERRQ(ierr); 453 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 454 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 455 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 456 ierr = DMRestoreLocalVector(dm, &localU);CHKERRQ(ierr); 457 PetscFunctionReturn(0); 458 } 459 460 #undef __FUNCT__ 461 #define __FUNCT__ "DMPlexInsertBoundaryValuesFEM" 462 PetscErrorCode DMPlexInsertBoundaryValuesFEM(DM dm, Vec localX) 463 { 464 void (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx); 465 void **ctxs; 466 PetscFE *fe; 467 PetscInt numFields, f, numBd, b; 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 472 PetscValidHeaderSpecific(localX, VEC_CLASSID, 2); 473 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 474 ierr = PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 475 for (f = 0; f < numFields; ++f) {ierr = DMGetField(dm, f, (PetscObject *) &fe[f]);CHKERRQ(ierr);} 476 /* OPT: Could attempt to do multiple BCs at once */ 477 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 478 for (b = 0; b < numBd; ++b) { 479 DMLabel label; 480 const PetscInt *ids; 481 const char *labelname; 482 PetscInt numids, field; 483 PetscBool isEssential; 484 void (*func)(); 485 void *ctx; 486 487 ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 488 ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr); 489 for (f = 0; f < numFields; ++f) { 490 funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL; 491 ctxs[f] = field == f ? ctx : NULL; 492 } 493 ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 494 } 495 ierr = PetscFree3(fe,funcs,ctxs);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 #undef __FUNCT__ 500 #define __FUNCT__ "DMPlexComputeL2Diff" 501 /*@C 502 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 503 504 Input Parameters: 505 + dm - The DM 506 . funcs - The functions to evaluate for each field component 507 . ctxs - Optional array of contexts to pass to each function, or NULL. 508 - X - The coefficient vector u_h 509 510 Output Parameter: 511 . diff - The diff ||u - u_h||_2 512 513 Level: developer 514 515 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff() 516 @*/ 517 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 518 { 519 const PetscInt debug = 0; 520 PetscSection section; 521 PetscQuadrature quad; 522 Vec localX; 523 PetscScalar *funcVal; 524 PetscReal *coords, *v0, *J, *invJ, detJ; 525 PetscReal localDiff = 0.0; 526 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 527 PetscErrorCode ierr; 528 529 PetscFunctionBegin; 530 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 531 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 532 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 533 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 534 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 535 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 536 for (field = 0; field < numFields; ++field) { 537 PetscFE fe; 538 PetscInt Nc; 539 540 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 541 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 542 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 543 numComponents += Nc; 544 } 545 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 546 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 547 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 548 for (c = cStart; c < cEnd; ++c) { 549 PetscScalar *x = NULL; 550 PetscReal elemDiff = 0.0; 551 552 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 553 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 554 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 555 556 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 557 PetscFE fe; 558 void * const ctx = ctxs ? ctxs[field] : NULL; 559 const PetscReal *quadPoints, *quadWeights; 560 PetscReal *basis; 561 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 562 563 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 564 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 565 ierr = PetscFEGetDimension(fe, &numBasisFuncs);CHKERRQ(ierr); 566 ierr = PetscFEGetNumComponents(fe, &numBasisComps);CHKERRQ(ierr); 567 ierr = PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);CHKERRQ(ierr); 568 if (debug) { 569 char title[1024]; 570 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 571 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 572 } 573 for (q = 0; q < numQuadPoints; ++q) { 574 for (d = 0; d < dim; d++) { 575 coords[d] = v0[d]; 576 for (e = 0; e < dim; e++) { 577 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 578 } 579 } 580 (*funcs[field])(coords, funcVal, ctx); 581 for (fc = 0; fc < numBasisComps; ++fc) { 582 PetscScalar interpolant = 0.0; 583 584 for (f = 0; f < numBasisFuncs; ++f) { 585 const PetscInt fidx = f*numBasisComps+fc; 586 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 587 } 588 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 589 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 590 } 591 } 592 comp += numBasisComps; 593 fieldOffset += numBasisFuncs*numBasisComps; 594 } 595 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 596 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 597 localDiff += elemDiff; 598 } 599 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 600 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 601 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 602 *diff = PetscSqrtReal(*diff); 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 608 /*@C 609 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 610 611 Input Parameters: 612 + dm - The DM 613 . funcs - The gradient functions to evaluate for each field component 614 . ctxs - Optional array of contexts to pass to each function, or NULL. 615 . X - The coefficient vector u_h 616 - n - The vector to project along 617 618 Output Parameter: 619 . diff - The diff ||(grad u - grad u_h) . n||_2 620 621 Level: developer 622 623 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 624 @*/ 625 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 626 { 627 const PetscInt debug = 0; 628 PetscSection section; 629 PetscQuadrature quad; 630 Vec localX; 631 PetscScalar *funcVal, *interpolantVec; 632 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 633 PetscReal localDiff = 0.0; 634 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 635 PetscErrorCode ierr; 636 637 PetscFunctionBegin; 638 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 639 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 640 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 641 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 642 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 643 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 644 for (field = 0; field < numFields; ++field) { 645 PetscFE fe; 646 PetscInt Nc; 647 648 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 649 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 650 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 651 numComponents += Nc; 652 } 653 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 654 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 655 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 656 for (c = cStart; c < cEnd; ++c) { 657 PetscScalar *x = NULL; 658 PetscReal elemDiff = 0.0; 659 660 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 661 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 662 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 663 664 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 665 PetscFE fe; 666 void * const ctx = ctxs ? ctxs[field] : NULL; 667 const PetscReal *quadPoints, *quadWeights; 668 PetscReal *basisDer; 669 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 670 671 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 672 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 673 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 674 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 675 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 676 if (debug) { 677 char title[1024]; 678 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 679 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 680 } 681 for (q = 0; q < numQuadPoints; ++q) { 682 for (d = 0; d < dim; d++) { 683 coords[d] = v0[d]; 684 for (e = 0; e < dim; e++) { 685 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 686 } 687 } 688 (*funcs[field])(coords, n, funcVal, ctx); 689 for (fc = 0; fc < Ncomp; ++fc) { 690 PetscScalar interpolant = 0.0; 691 692 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 693 for (f = 0; f < Nb; ++f) { 694 const PetscInt fidx = f*Ncomp+fc; 695 696 for (d = 0; d < dim; ++d) { 697 realSpaceDer[d] = 0.0; 698 for (g = 0; g < dim; ++g) { 699 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 700 } 701 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 702 } 703 } 704 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 705 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 706 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 707 } 708 } 709 comp += Ncomp; 710 fieldOffset += Nb*Ncomp; 711 } 712 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 713 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 714 localDiff += elemDiff; 715 } 716 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 717 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 718 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 719 *diff = PetscSqrtReal(*diff); 720 PetscFunctionReturn(0); 721 } 722 723 #undef __FUNCT__ 724 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 725 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 726 { 727 const PetscInt debug = 0; 728 PetscSection section; 729 PetscQuadrature quad; 730 Vec localX; 731 PetscScalar *funcVal; 732 PetscReal *coords, *v0, *J, *invJ, detJ; 733 PetscReal *localDiff; 734 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 735 PetscErrorCode ierr; 736 737 PetscFunctionBegin; 738 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 739 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 740 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 741 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 742 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 743 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 744 for (field = 0; field < numFields; ++field) { 745 PetscFE fe; 746 PetscInt Nc; 747 748 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 749 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 750 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 751 numComponents += Nc; 752 } 753 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 754 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 755 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 756 for (c = cStart; c < cEnd; ++c) { 757 PetscScalar *x = NULL; 758 759 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 760 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 761 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 762 763 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 764 PetscFE fe; 765 void * const ctx = ctxs ? ctxs[field] : NULL; 766 const PetscReal *quadPoints, *quadWeights; 767 PetscReal *basis, elemDiff = 0.0; 768 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 769 770 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 771 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 772 ierr = PetscFEGetDimension(fe, &numBasisFuncs);CHKERRQ(ierr); 773 ierr = PetscFEGetNumComponents(fe, &numBasisComps);CHKERRQ(ierr); 774 ierr = PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);CHKERRQ(ierr); 775 if (debug) { 776 char title[1024]; 777 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 778 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 779 } 780 for (q = 0; q < numQuadPoints; ++q) { 781 for (d = 0; d < dim; d++) { 782 coords[d] = v0[d]; 783 for (e = 0; e < dim; e++) { 784 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 785 } 786 } 787 (*funcs[field])(coords, funcVal, ctx); 788 for (fc = 0; fc < numBasisComps; ++fc) { 789 PetscScalar interpolant = 0.0; 790 791 for (f = 0; f < numBasisFuncs; ++f) { 792 const PetscInt fidx = f*numBasisComps+fc; 793 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 794 } 795 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 796 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 797 } 798 } 799 comp += numBasisComps; 800 fieldOffset += numBasisFuncs*numBasisComps; 801 localDiff[field] += elemDiff; 802 } 803 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 804 } 805 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 806 ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 807 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 808 ierr = PetscFree6(localDiff,funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 809 PetscFunctionReturn(0); 810 } 811 812 #undef __FUNCT__ 813 #define __FUNCT__ "DMPlexComputeIntegralFEM" 814 /*@ 815 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 816 817 Input Parameters: 818 + dm - The mesh 819 . X - Local input vector 820 - user - The user context 821 822 Output Parameter: 823 . integral - Local integral for each field 824 825 Level: developer 826 827 .seealso: DMPlexComputeResidualFEM() 828 @*/ 829 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 830 { 831 DM_Plex *mesh = (DM_Plex *) dm->data; 832 DM dmAux; 833 Vec localX, A; 834 PetscProblem prob, probAux = NULL; 835 PetscQuadrature q; 836 PetscCellGeometry geom; 837 PetscSection section, sectionAux; 838 PetscReal *v0, *J, *invJ, *detJ; 839 PetscScalar *u, *a = NULL; 840 PetscInt dim, Nf, f, numCells, cStart, cEnd, c; 841 PetscInt totDim, totDimAux; 842 PetscErrorCode ierr; 843 844 PetscFunctionBegin; 845 /*ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 846 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 847 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 848 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 849 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 850 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 851 ierr = DMGetProblem(dm, &prob);CHKERRQ(ierr); 852 ierr = PetscProblemGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 853 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 854 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 855 numCells = cEnd - cStart; 856 for (f = 0; f < Nf; ++f) {integral[f] = 0.0;} 857 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 858 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 859 if (dmAux) { 860 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 861 ierr = DMGetProblem(dmAux, &probAux);CHKERRQ(ierr); 862 ierr = PetscProblemGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 863 } 864 ierr = DMPlexInsertBoundaryValuesFEM(dm, localX);CHKERRQ(ierr); 865 ierr = PetscMalloc5(numCells*totDim,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ);CHKERRQ(ierr); 866 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 867 for (c = cStart; c < cEnd; ++c) { 868 PetscScalar *x = NULL; 869 PetscInt i; 870 871 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 872 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 873 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 874 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 875 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 876 if (dmAux) { 877 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 878 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 879 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 880 } 881 } 882 for (f = 0; f < Nf; ++f) { 883 PetscFE fe; 884 PetscInt numQuadPoints, Nb; 885 /* Conforming batches */ 886 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 887 /* Remainder */ 888 PetscInt Nr, offset; 889 890 ierr = PetscProblemGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 891 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 892 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 893 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 894 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 895 blockSize = Nb*numQuadPoints; 896 batchSize = numBlocks * blockSize; 897 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 898 numChunks = numCells / (numBatches*batchSize); 899 Ne = numChunks*numBatches*batchSize; 900 Nr = numCells % (numBatches*batchSize); 901 offset = numCells - Nr; 902 geom.v0 = v0; 903 geom.J = J; 904 geom.invJ = invJ; 905 geom.detJ = detJ; 906 ierr = PetscFEIntegrate(fe, prob, f, Ne, geom, u, probAux, a, integral);CHKERRQ(ierr); 907 geom.v0 = &v0[offset*dim]; 908 geom.J = &J[offset*dim*dim]; 909 geom.invJ = &invJ[offset*dim*dim]; 910 geom.detJ = &detJ[offset]; 911 ierr = PetscFEIntegrate(fe, prob, f, Nr, geom, &u[offset*totDim], probAux, &a[offset*totDimAux], integral);CHKERRQ(ierr); 912 } 913 ierr = PetscFree5(u,v0,J,invJ,detJ);CHKERRQ(ierr); 914 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 915 if (mesh->printFEM) { 916 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 917 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 918 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 919 } 920 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 921 /* TODO: Allreduce for integral */ 922 /*ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 923 PetscFunctionReturn(0); 924 } 925 926 #undef __FUNCT__ 927 #define __FUNCT__ "DMPlexComputeResidualFEM_Internal" 928 PetscErrorCode DMPlexComputeResidualFEM_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user) 929 { 930 DM_Plex *mesh = (DM_Plex *) dm->data; 931 const char *name = "Residual"; 932 DM dmAux; 933 DMLabel depth; 934 Vec A; 935 PetscProblem prob, probAux = NULL; 936 PetscQuadrature q; 937 PetscCellGeometry geom; 938 PetscSection section, sectionAux; 939 PetscReal *v0, *J, *invJ, *detJ; 940 PetscScalar *elemVec, *u, *u_t, *a = NULL; 941 PetscInt dim, Nf, f, numCells, cStart, cEnd, c, numBd, bd; 942 PetscInt totDim, totDimBd, totDimAux; 943 PetscErrorCode ierr; 944 945 PetscFunctionBegin; 946 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 947 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 948 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 949 ierr = DMGetProblem(dm, &prob);CHKERRQ(ierr); 950 ierr = PetscProblemGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 951 ierr = PetscProblemGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr); 952 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 953 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 954 numCells = cEnd - cStart; 955 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 956 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 957 if (dmAux) { 958 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 959 ierr = DMGetProblem(dmAux, &probAux);CHKERRQ(ierr); 960 ierr = PetscProblemGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 961 } 962 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 963 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 964 ierr = PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim,&elemVec);CHKERRQ(ierr); 965 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 966 for (c = cStart; c < cEnd; ++c) { 967 PetscScalar *x = NULL, *x_t = NULL; 968 PetscInt i; 969 970 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 971 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 972 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 973 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 974 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 975 if (X_t) { 976 ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 977 for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i]; 978 ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 979 } 980 if (dmAux) { 981 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 982 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 983 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 984 } 985 } 986 for (f = 0; f < Nf; ++f) { 987 PetscFE fe; 988 PetscInt numQuadPoints, Nb; 989 /* Conforming batches */ 990 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 991 /* Remainder */ 992 PetscInt Nr, offset; 993 994 ierr = PetscProblemGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 995 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 996 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 997 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 998 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 999 blockSize = Nb*numQuadPoints; 1000 batchSize = numBlocks * blockSize; 1001 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1002 numChunks = numCells / (numBatches*batchSize); 1003 Ne = numChunks*numBatches*batchSize; 1004 Nr = numCells % (numBatches*batchSize); 1005 offset = numCells - Nr; 1006 geom.v0 = v0; 1007 geom.J = J; 1008 geom.invJ = invJ; 1009 geom.detJ = detJ; 1010 ierr = PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);CHKERRQ(ierr); 1011 geom.v0 = &v0[offset*dim]; 1012 geom.J = &J[offset*dim*dim]; 1013 geom.invJ = &invJ[offset*dim*dim]; 1014 geom.detJ = &detJ[offset]; 1015 ierr = PetscFEIntegrateResidual(fe, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);CHKERRQ(ierr); 1016 } 1017 for (c = cStart; c < cEnd; ++c) { 1018 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, totDim, &elemVec[c*totDim]);CHKERRQ(ierr);} 1019 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);CHKERRQ(ierr); 1020 } 1021 ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1022 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1023 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1024 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1025 for (bd = 0; bd < numBd; ++bd) { 1026 const char *bdLabel; 1027 DMLabel label; 1028 IS pointIS; 1029 const PetscInt *points; 1030 const PetscInt *values; 1031 PetscReal *n; 1032 PetscInt field, numValues, numPoints, p, dep, numFaces; 1033 PetscBool isEssential; 1034 1035 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 1036 if (isEssential) continue; 1037 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 1038 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 1039 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1040 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1041 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1042 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1043 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1044 if (dep == dim-1) ++numFaces; 1045 } 1046 ierr = PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*totDimBd,&elemVec);CHKERRQ(ierr); 1047 if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);} 1048 for (p = 0, f = 0; p < numPoints; ++p) { 1049 const PetscInt point = points[p]; 1050 PetscScalar *x = NULL; 1051 PetscInt i; 1052 1053 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1054 if (dep != dim-1) continue; 1055 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1056 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1057 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1058 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1059 for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i]; 1060 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1061 if (X_t) { 1062 ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1063 for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i]; 1064 ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1065 } 1066 ++f; 1067 } 1068 for (f = 0; f < Nf; ++f) { 1069 PetscFE fe; 1070 PetscInt numQuadPoints, Nb; 1071 /* Conforming batches */ 1072 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1073 /* Remainder */ 1074 PetscInt Nr, offset; 1075 1076 ierr = PetscProblemGetBdDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1077 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1078 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1079 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1080 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1081 blockSize = Nb*numQuadPoints; 1082 batchSize = numBlocks * blockSize; 1083 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1084 numChunks = numFaces / (numBatches*batchSize); 1085 Ne = numChunks*numBatches*batchSize; 1086 Nr = numFaces % (numBatches*batchSize); 1087 offset = numFaces - Nr; 1088 geom.v0 = v0; 1089 geom.n = n; 1090 geom.J = J; 1091 geom.invJ = invJ; 1092 geom.detJ = detJ; 1093 ierr = PetscFEIntegrateBdResidual(fe, prob, f, Ne, geom, u, u_t, NULL, NULL, elemVec);CHKERRQ(ierr); 1094 geom.v0 = &v0[offset*dim]; 1095 geom.n = &n[offset*dim]; 1096 geom.J = &J[offset*dim*dim]; 1097 geom.invJ = &invJ[offset*dim*dim]; 1098 geom.detJ = &detJ[offset]; 1099 ierr = PetscFEIntegrateBdResidual(fe, prob, f, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);CHKERRQ(ierr); 1100 } 1101 for (p = 0, f = 0; p < numPoints; ++p) { 1102 const PetscInt point = points[p]; 1103 1104 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1105 if (dep != dim-1) continue; 1106 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);CHKERRQ(ierr);} 1107 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*totDimBd], ADD_VALUES);CHKERRQ(ierr); 1108 ++f; 1109 } 1110 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1111 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1112 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1113 if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);} 1114 } 1115 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 1116 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1117 PetscFunctionReturn(0); 1118 } 1119 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "DMPlexSNESComputeResidualFEM" 1122 /*@ 1123 DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 1124 1125 Input Parameters: 1126 + dm - The mesh 1127 . X - Local solution 1128 - user - The user context 1129 1130 Output Parameter: 1131 . F - Local output vector 1132 1133 Level: developer 1134 1135 .seealso: DMPlexComputeJacobianActionFEM() 1136 @*/ 1137 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 1138 { 1139 PetscErrorCode ierr; 1140 1141 PetscFunctionBegin; 1142 ierr = DMPlexComputeResidualFEM_Internal(dm, X, NULL, F, user);CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 1148 /*@ 1149 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 1150 1151 Input Parameters: 1152 + dm - The mesh 1153 . t - The time 1154 . X - Local solution 1155 . X_t - Local solution time derivative, or NULL 1156 - user - The user context 1157 1158 Output Parameter: 1159 . F - Local output vector 1160 1161 Level: developer 1162 1163 .seealso: DMPlexComputeJacobianActionFEM() 1164 @*/ 1165 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user) 1166 { 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = DMPlexComputeResidualFEM_Internal(dm, X, X_t, F, user);CHKERRQ(ierr); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "DMPlexComputeJacobianFEM_Internal" 1176 PetscErrorCode DMPlexComputeJacobianFEM_Internal(DM dm, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user) 1177 { 1178 DM_Plex *mesh = (DM_Plex *) dm->data; 1179 const char *name = "Jacobian"; 1180 DM dmAux; 1181 DMLabel depth; 1182 Vec A; 1183 PetscProblem prob, probAux = NULL; 1184 PetscQuadrature quad; 1185 PetscCellGeometry geom; 1186 PetscSection section, globalSection, sectionAux; 1187 PetscReal *v0, *J, *invJ, *detJ; 1188 PetscScalar *elemMat, *u, *u_t, *a = NULL; 1189 PetscInt dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c; 1190 PetscInt totDim, totDimBd, totDimAux, numBd, bd; 1191 PetscBool isShell; 1192 PetscErrorCode ierr; 1193 1194 PetscFunctionBegin; 1195 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1196 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1197 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1198 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 1199 ierr = DMGetProblem(dm, &prob);CHKERRQ(ierr); 1200 ierr = PetscProblemGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1201 ierr = PetscProblemGetTotalBdDimension(prob, &totDimBd);CHKERRQ(ierr); 1202 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1203 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1204 numCells = cEnd - cStart; 1205 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1206 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1207 if (dmAux) { 1208 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1209 ierr = DMGetProblem(dmAux, &probAux);CHKERRQ(ierr); 1210 ierr = PetscProblemGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1211 } 1212 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 1213 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 1214 ierr = PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim*totDim,&elemMat);CHKERRQ(ierr); 1215 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1216 for (c = cStart; c < cEnd; ++c) { 1217 PetscScalar *x = NULL, *x_t = NULL; 1218 PetscInt i; 1219 1220 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1221 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1222 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1223 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1224 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1225 if (X_t) { 1226 ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 1227 for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i]; 1228 ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 1229 } 1230 if (dmAux) { 1231 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1232 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1233 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1234 } 1235 } 1236 ierr = PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));CHKERRQ(ierr); 1237 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1238 PetscFE fe; 1239 PetscInt numQuadPoints, Nb; 1240 /* Conforming batches */ 1241 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1242 /* Remainder */ 1243 PetscInt Nr, offset; 1244 1245 ierr = PetscProblemGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1246 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1247 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1248 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1249 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1250 blockSize = Nb*numQuadPoints; 1251 batchSize = numBlocks * blockSize; 1252 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1253 numChunks = numCells / (numBatches*batchSize); 1254 Ne = numChunks*numBatches*batchSize; 1255 Nr = numCells % (numBatches*batchSize); 1256 offset = numCells - Nr; 1257 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1258 geom.v0 = v0; 1259 geom.J = J; 1260 geom.invJ = invJ; 1261 geom.detJ = detJ; 1262 ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, probAux, a, elemMat);CHKERRQ(ierr); 1263 geom.v0 = &v0[offset*dim]; 1264 geom.J = &J[offset*dim*dim]; 1265 geom.invJ = &invJ[offset*dim*dim]; 1266 geom.detJ = &detJ[offset]; 1267 ierr = PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);CHKERRQ(ierr); 1268 } 1269 } 1270 for (c = cStart; c < cEnd; ++c) { 1271 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[c*totDim*totDim]);CHKERRQ(ierr);} 1272 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*totDim*totDim], ADD_VALUES);CHKERRQ(ierr); 1273 } 1274 ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1275 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1276 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1277 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1278 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1279 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1280 for (bd = 0; bd < numBd; ++bd) { 1281 const char *bdLabel; 1282 DMLabel label; 1283 IS pointIS; 1284 const PetscInt *points; 1285 const PetscInt *values; 1286 PetscReal *n; 1287 PetscInt field, numValues, numPoints, p, dep, numFaces; 1288 PetscBool isEssential; 1289 1290 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 1291 if (isEssential) continue; 1292 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 1293 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 1294 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1295 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1296 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1297 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1298 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1299 if (dep == dim-1) ++numFaces; 1300 } 1301 ierr = PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*totDimBd*totDimBd,&elemMat);CHKERRQ(ierr); 1302 if (X_t) {ierr = PetscMalloc1(numFaces*totDimBd,&u_t);CHKERRQ(ierr);} 1303 for (p = 0, f = 0; p < numPoints; ++p) { 1304 const PetscInt point = points[p]; 1305 PetscScalar *x = NULL; 1306 PetscInt i; 1307 1308 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1309 if (dep != dim-1) continue; 1310 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1311 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1312 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1313 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1314 for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i]; 1315 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1316 if (X_t) { 1317 ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1318 for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i]; 1319 ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1320 } 1321 ++f; 1322 } 1323 ierr = PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));CHKERRQ(ierr); 1324 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1325 PetscFE fe; 1326 PetscInt numQuadPoints, Nb; 1327 /* Conforming batches */ 1328 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1329 /* Remainder */ 1330 PetscInt Nr, offset; 1331 1332 ierr = PetscProblemGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr); 1333 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 1334 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1335 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1336 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1337 blockSize = Nb*numQuadPoints; 1338 batchSize = numBlocks * blockSize; 1339 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1340 numChunks = numFaces / (numBatches*batchSize); 1341 Ne = numChunks*numBatches*batchSize; 1342 Nr = numFaces % (numBatches*batchSize); 1343 offset = numFaces - Nr; 1344 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1345 geom.v0 = v0; 1346 geom.n = n; 1347 geom.J = J; 1348 geom.invJ = invJ; 1349 geom.detJ = detJ; 1350 ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, NULL, NULL, elemMat);CHKERRQ(ierr); 1351 geom.v0 = &v0[offset*dim]; 1352 geom.n = &n[offset*dim]; 1353 geom.J = &J[offset*dim*dim]; 1354 geom.invJ = &invJ[offset*dim*dim]; 1355 geom.detJ = &detJ[offset]; 1356 ierr = PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);CHKERRQ(ierr); 1357 } 1358 } 1359 for (p = 0, f = 0; p < numPoints; ++p) { 1360 const PetscInt point = points[p]; 1361 1362 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1363 if (dep != dim-1) continue; 1364 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);CHKERRQ(ierr);} 1365 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);CHKERRQ(ierr); 1366 ++f; 1367 } 1368 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1369 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1370 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1371 if (X_t) {ierr = PetscFree(u_t);CHKERRQ(ierr);} 1372 } 1373 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1374 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1375 if (mesh->printFEM) { 1376 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1377 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 1378 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1379 } 1380 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1381 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr); 1382 if (isShell) { 1383 JacActionCtx *jctx; 1384 1385 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1386 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 1387 } 1388 PetscFunctionReturn(0); 1389 } 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "DMPlexSNESComputeJacobianFEM" 1393 /*@ 1394 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 1395 1396 Input Parameters: 1397 + dm - The mesh 1398 . X - Local input vector 1399 - user - The user context 1400 1401 Output Parameter: 1402 . Jac - Jacobian matrix 1403 1404 Note: 1405 The first member of the user context must be an FEMContext. 1406 1407 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1408 like a GPU, or vectorize on a multicore machine. 1409 1410 Level: developer 1411 1412 .seealso: FormFunctionLocal() 1413 @*/ 1414 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user) 1415 { 1416 PetscErrorCode ierr; 1417 1418 PetscFunctionBegin; 1419 ierr = DMPlexComputeJacobianFEM_Internal(dm, X, NULL, Jac, JacP, user);CHKERRQ(ierr); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 #undef __FUNCT__ 1424 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1425 /*@ 1426 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1427 1428 Input Parameters: 1429 + dmf - The fine mesh 1430 . dmc - The coarse mesh 1431 - user - The user context 1432 1433 Output Parameter: 1434 . In - The interpolation matrix 1435 1436 Note: 1437 The first member of the user context must be an FEMContext. 1438 1439 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1440 like a GPU, or vectorize on a multicore machine. 1441 1442 Level: developer 1443 1444 .seealso: DMPlexComputeJacobianFEM() 1445 @*/ 1446 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1447 { 1448 DM_Plex *mesh = (DM_Plex *) dmc->data; 1449 const char *name = "Interpolator"; 1450 PetscProblem prob; 1451 PetscFE *feRef; 1452 PetscSection fsection, fglobalSection; 1453 PetscSection csection, cglobalSection; 1454 PetscScalar *elemMat; 1455 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 1456 PetscInt cTotDim, rTotDim = 0; 1457 PetscErrorCode ierr; 1458 1459 PetscFunctionBegin; 1460 #if 0 1461 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1462 #endif 1463 ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr); 1464 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1465 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1466 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1467 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1468 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1469 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1470 ierr = DMGetProblem(dmf, &prob);CHKERRQ(ierr); 1471 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1472 for (f = 0; f < Nf; ++f) { 1473 PetscFE fe; 1474 PetscInt rNb, Nc; 1475 1476 ierr = PetscProblemGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1477 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1478 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1479 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1480 rTotDim += rNb*Nc; 1481 } 1482 ierr = PetscProblemGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1483 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1484 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1485 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1486 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1487 PetscDualSpace Qref; 1488 PetscQuadrature f; 1489 const PetscReal *qpoints, *qweights; 1490 PetscReal *points; 1491 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1492 1493 /* Compose points from all dual basis functionals */ 1494 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1495 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1496 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1497 for (i = 0; i < fpdim; ++i) { 1498 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1499 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1500 npoints += Np; 1501 } 1502 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1503 for (i = 0, k = 0; i < fpdim; ++i) { 1504 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1505 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1506 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1507 } 1508 1509 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1510 PetscFE fe; 1511 PetscReal *B; 1512 PetscInt NcJ, cpdim, j; 1513 1514 /* Evaluate basis at points */ 1515 ierr = PetscProblemGetDiscretization(prob, fieldJ, (PetscObject *) &fe);CHKERRQ(ierr); 1516 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1517 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1518 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1519 /* For now, fields only interpolate themselves */ 1520 if (fieldI == fieldJ) { 1521 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1522 for (i = 0, k = 0; i < fpdim; ++i) { 1523 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1524 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1525 for (p = 0; p < Np; ++p, ++k) { 1526 for (j = 0; j < cpdim; ++j) { 1527 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1528 } 1529 } 1530 } 1531 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1532 } 1533 offsetJ += cpdim*NcJ; 1534 } 1535 offsetI += fpdim*Nc; 1536 ierr = PetscFree(points);CHKERRQ(ierr); 1537 } 1538 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1539 for (c = cStart; c < cEnd; ++c) { 1540 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1541 } 1542 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1543 ierr = PetscFree(feRef);CHKERRQ(ierr); 1544 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1545 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1546 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1547 if (mesh->printFEM) { 1548 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1549 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1550 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1551 } 1552 #if 0 1553 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1554 #endif 1555 PetscFunctionReturn(0); 1556 } 1557 1558 #undef __FUNCT__ 1559 #define __FUNCT__ "BoundaryDuplicate" 1560 static PetscErrorCode BoundaryDuplicate(DMBoundary bd, DMBoundary *boundary) 1561 { 1562 DMBoundary b = bd, b2, bold = NULL; 1563 PetscErrorCode ierr; 1564 1565 PetscFunctionBegin; 1566 *boundary = NULL; 1567 for (; b; b = b->next, bold = b2) { 1568 ierr = PetscNew(&b2);CHKERRQ(ierr); 1569 ierr = PetscStrallocpy(b->name, (char **) &b2->name);CHKERRQ(ierr); 1570 ierr = PetscStrallocpy(b->labelname, (char **) &b2->labelname);CHKERRQ(ierr); 1571 ierr = PetscMalloc1(b->numids, &b2->ids);CHKERRQ(ierr); 1572 ierr = PetscMemcpy(b2->ids, b->ids, b->numids*sizeof(PetscInt));CHKERRQ(ierr); 1573 b2->label = NULL; 1574 b2->essential = b->essential; 1575 b2->field = b->field; 1576 b2->func = b->func; 1577 b2->numids = b->numids; 1578 b2->ctx = b->ctx; 1579 b2->next = NULL; 1580 if (!*boundary) *boundary = b2; 1581 if (bold) bold->next = b2; 1582 } 1583 PetscFunctionReturn(0); 1584 } 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "DMPlexCopyBoundary" 1588 PetscErrorCode DMPlexCopyBoundary(DM dm, DM dmNew) 1589 { 1590 DM_Plex *mesh = (DM_Plex *) dm->data; 1591 DM_Plex *meshNew = (DM_Plex *) dmNew->data; 1592 DMBoundary b; 1593 PetscErrorCode ierr; 1594 1595 PetscFunctionBegin; 1596 ierr = BoundaryDuplicate(mesh->boundary, &meshNew->boundary);CHKERRQ(ierr); 1597 for (b = meshNew->boundary; b; b = b->next) { 1598 if (b->labelname) { 1599 ierr = DMPlexGetLabel(dmNew, b->labelname, &b->label);CHKERRQ(ierr); 1600 if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname); 1601 } 1602 } 1603 PetscFunctionReturn(0); 1604 } 1605 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "DMPlexAddBoundary" 1608 /* The ids can be overridden by the command line option -bc_<boundary name> */ 1609 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], const char labelname[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx) 1610 { 1611 DM_Plex *mesh = (DM_Plex *) dm->data; 1612 DMBoundary b; 1613 PetscErrorCode ierr; 1614 1615 PetscFunctionBegin; 1616 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1617 ierr = PetscNew(&b);CHKERRQ(ierr); 1618 ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr); 1619 ierr = PetscStrallocpy(labelname, (char **) &b->labelname);CHKERRQ(ierr); 1620 ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr); 1621 ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr); 1622 if (b->labelname) { 1623 ierr = DMPlexGetLabel(dm, b->labelname, &b->label);CHKERRQ(ierr); 1624 if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname); 1625 } 1626 b->essential = isEssential; 1627 b->field = field; 1628 b->func = bcFunc; 1629 b->numids = numids; 1630 b->ctx = ctx; 1631 b->next = mesh->boundary; 1632 mesh->boundary = b; 1633 PetscFunctionReturn(0); 1634 } 1635 1636 #undef __FUNCT__ 1637 #define __FUNCT__ "DMPlexGetNumBoundary" 1638 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd) 1639 { 1640 DM_Plex *mesh = (DM_Plex *) dm->data; 1641 DMBoundary b = mesh->boundary; 1642 1643 PetscFunctionBegin; 1644 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1645 PetscValidPointer(numBd, 2); 1646 *numBd = 0; 1647 while (b) {++(*numBd); b = b->next;} 1648 PetscFunctionReturn(0); 1649 } 1650 1651 #undef __FUNCT__ 1652 #define __FUNCT__ "DMPlexGetBoundary" 1653 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, const char **labelname, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx) 1654 { 1655 DM_Plex *mesh = (DM_Plex *) dm->data; 1656 DMBoundary b = mesh->boundary; 1657 PetscInt n = 0; 1658 1659 PetscFunctionBegin; 1660 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1661 while (b) { 1662 if (n == bd) break; 1663 b = b->next; 1664 ++n; 1665 } 1666 if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n); 1667 if (isEssential) { 1668 PetscValidPointer(isEssential, 3); 1669 *isEssential = b->essential; 1670 } 1671 if (name) { 1672 PetscValidPointer(name, 4); 1673 *name = b->name; 1674 } 1675 if (labelname) { 1676 PetscValidPointer(labelname, 5); 1677 *labelname = b->labelname; 1678 } 1679 if (field) { 1680 PetscValidPointer(field, 6); 1681 *field = b->field; 1682 } 1683 if (func) { 1684 PetscValidPointer(func, 7); 1685 *func = b->func; 1686 } 1687 if (numids) { 1688 PetscValidPointer(numids, 8); 1689 *numids = b->numids; 1690 } 1691 if (ids) { 1692 PetscValidPointer(ids, 9); 1693 *ids = b->ids; 1694 } 1695 if (ctx) { 1696 PetscValidPointer(ctx, 10); 1697 *ctx = b->ctx; 1698 } 1699 PetscFunctionReturn(0); 1700 } 1701 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "DMPlexIsBoundaryPoint" 1704 PetscErrorCode DMPlexIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd) 1705 { 1706 DM_Plex *mesh = (DM_Plex *) dm->data; 1707 DMBoundary b = mesh->boundary; 1708 PetscErrorCode ierr; 1709 1710 PetscFunctionBegin; 1711 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1712 PetscValidPointer(isBd, 3); 1713 *isBd = PETSC_FALSE; 1714 while (b && !(*isBd)) { 1715 if (b->label) { 1716 PetscInt i; 1717 1718 for (i = 0; i < b->numids && !(*isBd); ++i) { 1719 ierr = DMLabelStratumHasPoint(b->label, b->ids[i], point, isBd);CHKERRQ(ierr); 1720 } 1721 } 1722 b = b->next; 1723 } 1724 PetscFunctionReturn(0); 1725 } 1726