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