1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexGetScale" 5 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 6 { 7 DM_Plex *mesh = (DM_Plex*) dm->data; 8 9 PetscFunctionBegin; 10 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 11 PetscValidPointer(scale, 3); 12 *scale = mesh->scale[unit]; 13 PetscFunctionReturn(0); 14 } 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "DMPlexSetScale" 18 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 19 { 20 DM_Plex *mesh = (DM_Plex*) dm->data; 21 22 PetscFunctionBegin; 23 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 24 mesh->scale[unit] = scale; 25 PetscFunctionReturn(0); 26 } 27 28 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 29 { 30 switch (i) { 31 case 0: 32 switch (j) { 33 case 0: return 0; 34 case 1: 35 switch (k) { 36 case 0: return 0; 37 case 1: return 0; 38 case 2: return 1; 39 } 40 case 2: 41 switch (k) { 42 case 0: return 0; 43 case 1: return -1; 44 case 2: return 0; 45 } 46 } 47 case 1: 48 switch (j) { 49 case 0: 50 switch (k) { 51 case 0: return 0; 52 case 1: return 0; 53 case 2: return -1; 54 } 55 case 1: return 0; 56 case 2: 57 switch (k) { 58 case 0: return 1; 59 case 1: return 0; 60 case 2: return 0; 61 } 62 } 63 case 2: 64 switch (j) { 65 case 0: 66 switch (k) { 67 case 0: return 0; 68 case 1: return 1; 69 case 2: return 0; 70 } 71 case 1: 72 switch (k) { 73 case 0: return -1; 74 case 1: return 0; 75 case 2: return 0; 76 } 77 case 2: return 0; 78 } 79 } 80 return 0; 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "DMPlexCreateRigidBody" 85 /*@C 86 DMPlexCreateRigidBody - create rigid body modes from coordinates 87 88 Collective on DM 89 90 Input Arguments: 91 + dm - the DM 92 . section - the local section associated with the rigid field, or NULL for the default section 93 - globalSection - the global section associated with the rigid field, or NULL for the default section 94 95 Output Argument: 96 . sp - the null space 97 98 Note: This is necessary to take account of Dirichlet conditions on the displacements 99 100 Level: advanced 101 102 .seealso: MatNullSpaceCreate() 103 @*/ 104 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp) 105 { 106 MPI_Comm comm; 107 Vec coordinates, localMode, mode[6]; 108 PetscSection coordSection; 109 PetscScalar *coords; 110 PetscInt dim, vStart, vEnd, v, n, m, d, i, j; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 115 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 116 if (dim == 1) { 117 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 121 if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);} 122 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 123 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 124 ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 125 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 126 m = (dim*(dim+1))/2; 127 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 128 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 129 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 130 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 131 /* Assume P1 */ 132 ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr); 133 for (d = 0; d < dim; ++d) { 134 PetscScalar values[3] = {0.0, 0.0, 0.0}; 135 136 values[d] = 1.0; 137 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 138 for (v = vStart; v < vEnd; ++v) { 139 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 140 } 141 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 142 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 143 } 144 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 145 for (d = dim; d < dim*(dim+1)/2; ++d) { 146 PetscInt i, j, k = dim > 2 ? d - dim : d; 147 148 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 149 for (v = vStart; v < vEnd; ++v) { 150 PetscScalar values[3] = {0.0, 0.0, 0.0}; 151 PetscInt off; 152 153 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 154 for (i = 0; i < dim; ++i) { 155 for (j = 0; j < dim; ++j) { 156 values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]); 157 } 158 } 159 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 160 } 161 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 162 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 163 } 164 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 165 ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr); 166 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 167 /* Orthonormalize system */ 168 for (i = dim; i < m; ++i) { 169 PetscScalar dots[6]; 170 171 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 172 for (j = 0; j < i; ++j) dots[j] *= -1.0; 173 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 174 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 175 } 176 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 177 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 178 PetscFunctionReturn(0); 179 } 180 /******************************************************************************* 181 This should be in a separate Discretization object, but I am not sure how to lay 182 it out yet, so I am stuffing things here while I experiment. 183 *******************************************************************************/ 184 #undef __FUNCT__ 185 #define __FUNCT__ "DMPlexSetFEMIntegration" 186 PetscErrorCode DMPlexSetFEMIntegration(DM dm, 187 PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], 188 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 189 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 190 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]), 191 PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[], 192 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 193 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 194 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 195 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 196 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]), 197 PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], 198 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 199 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 200 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 201 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 202 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[])) 203 { 204 DM_Plex *mesh = (DM_Plex*) dm->data; 205 206 PetscFunctionBegin; 207 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 208 mesh->integrateResidualFEM = integrateResidualFEM; 209 mesh->integrateJacobianActionFEM = integrateJacobianActionFEM; 210 mesh->integrateJacobianFEM = integrateJacobianFEM; 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "DMPlexProjectFunctionLocal" 216 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec localX) 217 { 218 Vec coordinates; 219 PetscSection section, cSection; 220 PetscInt dim, vStart, vEnd, v, c, d; 221 PetscScalar *values, *cArray; 222 PetscReal *coords; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 227 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 228 ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr); 229 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 230 ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr); 231 ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr); 232 ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr); 233 ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr); 234 for (v = vStart; v < vEnd; ++v) { 235 PetscInt dof, off; 236 237 ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr); 238 ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr); 239 if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim); 240 for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]); 241 for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords); 242 ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr); 243 } 244 ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr); 245 /* Temporary, must be replaced by a projection on the finite element basis */ 246 { 247 PetscInt eStart = 0, eEnd = 0, e, depth; 248 249 ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr); 250 --depth; 251 if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);} 252 for (e = eStart; e < eEnd; ++e) { 253 const PetscInt *cone = NULL; 254 PetscInt coneSize, d; 255 PetscScalar *coordsA, *coordsB; 256 257 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 258 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 259 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e); 260 ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr); 261 ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr); 262 for (d = 0; d < dim; ++d) { 263 coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d])); 264 } 265 for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords); 266 ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr); 267 } 268 } 269 270 ierr = PetscFree(coords);CHKERRQ(ierr); 271 ierr = PetscFree(values);CHKERRQ(ierr); 272 #if 0 273 const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin()); 274 PetscReal detJ; 275 276 ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr); 277 ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr); 278 ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true); 279 280 for (PetscInt c = cStart; c < cEnd; ++c) { 281 ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV); 282 const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints(); 283 const int oSize = pV.getSize(); 284 int v = 0; 285 286 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 287 for (PetscInt cl = 0; cl < oSize; ++cl) { 288 const PetscInt fDim; 289 290 ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr); 291 if (pointDim) { 292 for (PetscInt d = 0; d < fDim; ++d, ++v) { 293 values[v] = (*this->_options.integrate)(v0, J, v, initFunc); 294 } 295 } 296 } 297 ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr); 298 pV.clear(); 299 } 300 ierr = PetscFree2(v0,J);CHKERRQ(ierr); 301 ierr = PetscFree(values);CHKERRQ(ierr); 302 #endif 303 PetscFunctionReturn(0); 304 } 305 306 #undef __FUNCT__ 307 #define __FUNCT__ "DMPlexProjectFunction" 308 /*@C 309 DMPlexProjectFunction - This projects the given function into the function space provided. 310 311 Input Parameters: 312 + dm - The DM 313 . numComp - The number of components (functions) 314 . funcs - The coordinate functions to evaluate 315 - mode - The insertion mode for values 316 317 Output Parameter: 318 . X - vector 319 320 Level: developer 321 322 Note: 323 This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector. 324 We will eventually fix it. 325 326 .seealso: DMPlexComputeL2Diff() 327 @*/ 328 PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec X) 329 { 330 Vec localX; 331 PetscErrorCode ierr; 332 333 PetscFunctionBegin; 334 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 335 ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr); 336 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 337 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 338 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "DMPlexComputeL2Diff" 344 /*@C 345 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 346 347 Input Parameters: 348 + dm - The DM 349 . quad - The PetscQuadrature object for each field 350 . funcs - The functions to evaluate for each field component 351 - X - The coefficient vector u_h 352 353 Output Parameter: 354 . diff - The diff ||u - u_h||_2 355 356 Level: developer 357 358 .seealso: DMPlexProjectFunction() 359 @*/ 360 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], PetscScalar (**funcs)(const PetscReal []), Vec X, PetscReal *diff) 361 { 362 const PetscInt debug = 0; 363 PetscSection section; 364 Vec localX; 365 PetscReal *coords, *v0, *J, *invJ, detJ; 366 PetscReal localDiff = 0.0; 367 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 368 PetscErrorCode ierr; 369 370 PetscFunctionBegin; 371 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 372 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 373 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 374 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 375 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 376 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 377 for (field = 0; field < numFields; ++field) { 378 numComponents += quad[field].numComponents; 379 } 380 ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 381 ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr); 382 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 383 for (c = cStart; c < cEnd; ++c) { 384 PetscScalar *x; 385 PetscReal elemDiff = 0.0; 386 387 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 388 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 389 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 390 391 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 392 const PetscInt numQuadPoints = quad[field].numQuadPoints; 393 const PetscReal *quadPoints = quad[field].quadPoints; 394 const PetscReal *quadWeights = quad[field].quadWeights; 395 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 396 const PetscInt numBasisComps = quad[field].numComponents; 397 const PetscReal *basis = quad[field].basis; 398 PetscInt q, d, e, fc, f; 399 400 if (debug) { 401 char title[1024]; 402 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 403 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 404 } 405 for (q = 0; q < numQuadPoints; ++q) { 406 for (d = 0; d < dim; d++) { 407 coords[d] = v0[d]; 408 for (e = 0; e < dim; e++) { 409 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 410 } 411 } 412 for (fc = 0; fc < numBasisComps; ++fc) { 413 const PetscReal funcVal = PetscRealPart((*funcs[comp+fc])(coords)); 414 PetscReal interpolant = 0.0; 415 for (f = 0; f < numBasisFuncs; ++f) { 416 const PetscInt fidx = f*numBasisComps+fc; 417 interpolant += PetscRealPart(x[fieldOffset+fidx])*basis[q*numBasisFuncs*numBasisComps+fidx]; 418 } 419 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ);CHKERRQ(ierr);} 420 elemDiff += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ; 421 } 422 } 423 comp += numBasisComps; 424 fieldOffset += numBasisFuncs*numBasisComps; 425 } 426 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 427 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 428 localDiff += elemDiff; 429 } 430 ierr = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr); 431 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 432 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);CHKERRQ(ierr); 433 *diff = PetscSqrtReal(*diff); 434 PetscFunctionReturn(0); 435 } 436 437 #undef __FUNCT__ 438 #define __FUNCT__ "DMPlexComputeResidualFEM" 439 /*@ 440 DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 441 442 Input Parameters: 443 + dm - The mesh 444 . X - Local input vector 445 - user - The user context 446 447 Output Parameter: 448 . F - Local output vector 449 450 Note: 451 The second member of the user context must be an FEMContext. 452 453 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 454 like a GPU, or vectorize on a multicore machine. 455 456 Level: developer 457 458 .seealso: DMPlexComputeJacobianActionFEM() 459 @*/ 460 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 461 { 462 DM_Plex *mesh = (DM_Plex*) dm->data; 463 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 464 PetscQuadrature *quad = fem->quad; 465 PetscSection section; 466 PetscReal *v0, *J, *invJ, *detJ; 467 PetscScalar *elemVec, *u; 468 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 469 PetscInt cellDof = 0, numComponents = 0; 470 PetscErrorCode ierr; 471 472 PetscFunctionBegin; 473 /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 474 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 475 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 476 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 477 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 478 numCells = cEnd - cStart; 479 for (field = 0; field < numFields; ++field) { 480 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 481 numComponents += quad[field].numComponents; 482 } 483 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 484 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 485 ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr); 486 for (c = cStart; c < cEnd; ++c) { 487 PetscScalar *x; 488 PetscInt i; 489 490 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 491 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 492 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 493 494 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 495 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 496 } 497 for (field = 0; field < numFields; ++field) { 498 const PetscInt numQuadPoints = quad[field].numQuadPoints; 499 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 500 void (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0Funcs[field]; 501 void (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1Funcs[field]; 502 /* Conforming batches */ 503 PetscInt blockSize = numBasisFuncs*numQuadPoints; 504 PetscInt numBlocks = 1; 505 PetscInt batchSize = numBlocks * blockSize; 506 PetscInt numBatches = numBatchesTmp; 507 PetscInt numChunks = numCells / (numBatches*batchSize); 508 /* Remainder */ 509 PetscInt numRemainder = numCells % (numBatches * batchSize); 510 PetscInt offset = numCells - numRemainder; 511 512 ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr); 513 ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 514 f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 515 } 516 for (c = cStart; c < cEnd; ++c) { 517 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 518 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 519 } 520 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 521 if (mesh->printFEM) { 522 PetscMPIInt rank, numProcs; 523 PetscInt p; 524 525 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 526 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 527 ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr); 528 for (p = 0; p < numProcs; ++p) { 529 if (p == rank) { 530 Vec f; 531 532 ierr = VecDuplicate(F, &f);CHKERRQ(ierr); 533 ierr = VecCopy(F, f);CHKERRQ(ierr); 534 ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr); 535 ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 536 ierr = VecDestroy(&f);CHKERRQ(ierr); 537 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 538 } 539 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 540 } 541 } 542 /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 543 PetscFunctionReturn(0); 544 } 545 546 #undef __FUNCT__ 547 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 548 /*@C 549 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 550 551 Input Parameters: 552 + dm - The mesh 553 . J - The Jacobian shell matrix 554 . X - Local input vector 555 - user - The user context 556 557 Output Parameter: 558 . F - Local output vector 559 560 Note: 561 The second member of the user context must be an FEMContext. 562 563 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 564 like a GPU, or vectorize on a multicore machine. 565 566 Level: developer 567 568 .seealso: DMPlexComputeResidualFEM() 569 @*/ 570 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 571 { 572 DM_Plex *mesh = (DM_Plex*) dm->data; 573 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 574 PetscQuadrature *quad = fem->quad; 575 PetscSection section; 576 JacActionCtx *jctx; 577 PetscReal *v0, *J, *invJ, *detJ; 578 PetscScalar *elemVec, *u, *a; 579 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 580 PetscInt cellDof = 0; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 585 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 586 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 587 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 588 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 589 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 590 numCells = cEnd - cStart; 591 for (field = 0; field < numFields; ++field) { 592 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 593 } 594 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 595 ierr = PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr); 596 for (c = cStart; c < cEnd; ++c) { 597 PetscScalar *x; 598 PetscInt i; 599 600 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 601 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 602 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 603 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 604 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 605 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 606 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 607 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 608 } 609 for (field = 0; field < numFields; ++field) { 610 const PetscInt numQuadPoints = quad[field].numQuadPoints; 611 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 612 /* Conforming batches */ 613 PetscInt blockSize = numBasisFuncs*numQuadPoints; 614 PetscInt numBlocks = 1; 615 PetscInt batchSize = numBlocks * blockSize; 616 PetscInt numBatches = numBatchesTmp; 617 PetscInt numChunks = numCells / (numBatches*batchSize); 618 /* Remainder */ 619 PetscInt numRemainder = numCells % (numBatches * batchSize); 620 PetscInt offset = numCells - numRemainder; 621 622 ierr = (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr); 623 ierr = (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 624 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 625 } 626 for (c = cStart; c < cEnd; ++c) { 627 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 628 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 629 } 630 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 631 if (mesh->printFEM) { 632 PetscMPIInt rank, numProcs; 633 PetscInt p; 634 635 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 636 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 637 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action:\n");CHKERRQ(ierr); 638 for (p = 0; p < numProcs; ++p) { 639 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 640 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 641 } 642 } 643 /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 644 PetscFunctionReturn(0); 645 } 646 647 #undef __FUNCT__ 648 #define __FUNCT__ "DMPlexComputeJacobianFEM" 649 /*@ 650 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 651 652 Input Parameters: 653 + dm - The mesh 654 . X - Local input vector 655 - user - The user context 656 657 Output Parameter: 658 . Jac - Jacobian matrix 659 660 Note: 661 The second member of the user context must be an FEMContext. 662 663 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 664 like a GPU, or vectorize on a multicore machine. 665 666 Level: developer 667 668 .seealso: FormFunctionLocal() 669 @*/ 670 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user) 671 { 672 DM_Plex *mesh = (DM_Plex*) dm->data; 673 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 674 PetscQuadrature *quad = fem->quad; 675 PetscSection section; 676 PetscReal *v0, *J, *invJ, *detJ; 677 PetscScalar *elemMat, *u; 678 PetscInt dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c; 679 PetscInt cellDof = 0, numComponents = 0; 680 PetscBool isShell; 681 PetscErrorCode ierr; 682 683 PetscFunctionBegin; 684 /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 685 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 686 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 687 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 688 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 689 numCells = cEnd - cStart; 690 for (field = 0; field < numFields; ++field) { 691 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 692 numComponents += quad[field].numComponents; 693 } 694 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 695 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 696 ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof*cellDof,PetscScalar,&elemMat);CHKERRQ(ierr); 697 for (c = cStart; c < cEnd; ++c) { 698 PetscScalar *x; 699 PetscInt i; 700 701 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 702 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 703 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 704 705 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 706 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 707 } 708 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 709 for (fieldI = 0; fieldI < numFields; ++fieldI) { 710 const PetscInt numQuadPoints = quad[fieldI].numQuadPoints; 711 const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs; 712 PetscInt fieldJ; 713 714 for (fieldJ = 0; fieldJ < numFields; ++fieldJ) { 715 void (*g0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]) = fem->g0Funcs[fieldI*numFields+fieldJ]; 716 void (*g1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]) = fem->g1Funcs[fieldI*numFields+fieldJ]; 717 void (*g2)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]) = fem->g2Funcs[fieldI*numFields+fieldJ]; 718 void (*g3)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]) = fem->g3Funcs[fieldI*numFields+fieldJ]; 719 /* Conforming batches */ 720 PetscInt blockSize = numBasisFuncs*numQuadPoints; 721 PetscInt numBlocks = 1; 722 PetscInt batchSize = numBlocks * blockSize; 723 PetscInt numBatches = numBatchesTmp; 724 PetscInt numChunks = numCells / (numBatches*batchSize); 725 /* Remainder */ 726 PetscInt numRemainder = numCells % (numBatches * batchSize); 727 PetscInt offset = numCells - numRemainder; 728 729 ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 730 ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 731 g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr); 732 } 733 } 734 for (c = cStart; c < cEnd; ++c) { 735 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 736 ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 737 } 738 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 739 740 /* Assemble matrix, using the 2-step process: 741 MatAssemblyBegin(), MatAssemblyEnd(). */ 742 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 743 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 744 745 if (mesh->printFEM) { 746 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr); 747 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 748 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 749 } 750 /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 751 ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr); 752 if (isShell) { 753 JacActionCtx *jctx; 754 755 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 756 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 757 } 758 *str = SAME_NONZERO_PATTERN; 759 PetscFunctionReturn(0); 760 } 761