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 (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], 192 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 193 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), 194 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]), 195 PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[], 196 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 197 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 198 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 199 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 200 void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]), 201 PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], 202 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], 203 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 204 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 205 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), 206 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[])) 207 { 208 DM_Plex *mesh = (DM_Plex*) dm->data; 209 210 PetscFunctionBegin; 211 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 212 mesh->integrateResidualFEM = integrateResidualFEM; 213 mesh->integrateBdResidualFEM = integrateBdResidualFEM; 214 mesh->integrateJacobianActionFEM = integrateJacobianActionFEM; 215 mesh->integrateJacobianFEM = integrateJacobianFEM; 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "DMPlexProjectFunctionLocal" 221 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec localX) 222 { 223 Vec coordinates; 224 PetscSection section, cSection; 225 PetscInt dim, vStart, vEnd, v, c, d; 226 PetscScalar *values, *cArray; 227 PetscReal *coords; 228 PetscErrorCode ierr; 229 230 PetscFunctionBegin; 231 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 232 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 233 ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr); 234 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 235 ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr); 236 ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr); 237 ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr); 238 ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr); 239 for (v = vStart; v < vEnd; ++v) { 240 PetscInt dof, off; 241 242 ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr); 243 ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr); 244 if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim); 245 for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]); 246 for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords); 247 ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr); 248 } 249 ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr); 250 /* Temporary, must be replaced by a projection on the finite element basis */ 251 { 252 PetscInt eStart = 0, eEnd = 0, e, depth; 253 254 ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr); 255 --depth; 256 if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);} 257 for (e = eStart; e < eEnd; ++e) { 258 const PetscInt *cone = NULL; 259 PetscInt coneSize, d; 260 PetscScalar *coordsA, *coordsB; 261 262 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 263 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 264 if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e); 265 ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr); 266 ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr); 267 for (d = 0; d < dim; ++d) { 268 coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d])); 269 } 270 for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords); 271 ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr); 272 } 273 } 274 275 ierr = PetscFree(coords);CHKERRQ(ierr); 276 ierr = PetscFree(values);CHKERRQ(ierr); 277 #if 0 278 const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin()); 279 PetscReal detJ; 280 281 ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr); 282 ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr); 283 ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true); 284 285 for (PetscInt c = cStart; c < cEnd; ++c) { 286 ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV); 287 const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints(); 288 const int oSize = pV.getSize(); 289 int v = 0; 290 291 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 292 for (PetscInt cl = 0; cl < oSize; ++cl) { 293 const PetscInt fDim; 294 295 ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr); 296 if (pointDim) { 297 for (PetscInt d = 0; d < fDim; ++d, ++v) { 298 values[v] = (*this->_options.integrate)(v0, J, v, initFunc); 299 } 300 } 301 } 302 ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr); 303 pV.clear(); 304 } 305 ierr = PetscFree2(v0,J);CHKERRQ(ierr); 306 ierr = PetscFree(values);CHKERRQ(ierr); 307 #endif 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "DMPlexProjectFunction" 313 /*@C 314 DMPlexProjectFunction - This projects the given function into the function space provided. 315 316 Input Parameters: 317 + dm - The DM 318 . numComp - The number of components (functions) 319 . funcs - The coordinate functions to evaluate 320 - mode - The insertion mode for values 321 322 Output Parameter: 323 . X - vector 324 325 Level: developer 326 327 Note: 328 This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector. 329 We will eventually fix it. 330 331 .seealso: DMPlexComputeL2Diff() 332 @*/ 333 PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec X) 334 { 335 Vec localX; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 340 ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr); 341 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 342 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 343 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "DMPlexComputeL2Diff" 349 /*@C 350 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 351 352 Input Parameters: 353 + dm - The DM 354 . quad - The PetscQuadrature object for each field 355 . funcs - The functions to evaluate for each field component 356 - X - The coefficient vector u_h 357 358 Output Parameter: 359 . diff - The diff ||u - u_h||_2 360 361 Level: developer 362 363 .seealso: DMPlexProjectFunction() 364 @*/ 365 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], PetscScalar (**funcs)(const PetscReal []), Vec X, PetscReal *diff) 366 { 367 const PetscInt debug = 0; 368 PetscSection section; 369 Vec localX; 370 PetscReal *coords, *v0, *J, *invJ, detJ; 371 PetscReal localDiff = 0.0; 372 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 373 PetscErrorCode ierr; 374 375 PetscFunctionBegin; 376 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 377 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 378 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 379 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 380 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 381 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 382 for (field = 0; field < numFields; ++field) { 383 numComponents += quad[field].numComponents; 384 } 385 ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 386 ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr); 387 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 388 for (c = cStart; c < cEnd; ++c) { 389 PetscScalar *x; 390 PetscReal elemDiff = 0.0; 391 392 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 393 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 394 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 395 396 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 397 const PetscInt numQuadPoints = quad[field].numQuadPoints; 398 const PetscReal *quadPoints = quad[field].quadPoints; 399 const PetscReal *quadWeights = quad[field].quadWeights; 400 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 401 const PetscInt numBasisComps = quad[field].numComponents; 402 const PetscReal *basis = quad[field].basis; 403 PetscInt q, d, e, fc, f; 404 405 if (debug) { 406 char title[1024]; 407 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 408 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 409 } 410 for (q = 0; q < numQuadPoints; ++q) { 411 for (d = 0; d < dim; d++) { 412 coords[d] = v0[d]; 413 for (e = 0; e < dim; e++) { 414 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 415 } 416 } 417 for (fc = 0; fc < numBasisComps; ++fc) { 418 const PetscReal funcVal = PetscRealPart((*funcs[comp+fc])(coords)); 419 PetscReal interpolant = 0.0; 420 for (f = 0; f < numBasisFuncs; ++f) { 421 const PetscInt fidx = f*numBasisComps+fc; 422 interpolant += PetscRealPart(x[fieldOffset+fidx])*basis[q*numBasisFuncs*numBasisComps+fidx]; 423 } 424 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ);CHKERRQ(ierr);} 425 elemDiff += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ; 426 } 427 } 428 comp += numBasisComps; 429 fieldOffset += numBasisFuncs*numBasisComps; 430 } 431 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 432 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 433 localDiff += elemDiff; 434 } 435 ierr = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr); 436 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 437 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);CHKERRQ(ierr); 438 *diff = PetscSqrtReal(*diff); 439 PetscFunctionReturn(0); 440 } 441 442 #undef __FUNCT__ 443 #define __FUNCT__ "DMPlexComputeResidualFEM" 444 /*@ 445 DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 446 447 Input Parameters: 448 + dm - The mesh 449 . X - Local input vector 450 - user - The user context 451 452 Output Parameter: 453 . F - Local output vector 454 455 Note: 456 The second member of the user context must be an FEMContext. 457 458 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 459 like a GPU, or vectorize on a multicore machine. 460 461 Level: developer 462 463 .seealso: DMPlexComputeJacobianActionFEM() 464 @*/ 465 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 466 { 467 DM_Plex *mesh = (DM_Plex*) dm->data; 468 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 469 PetscQuadrature *quad = fem->quad; 470 PetscQuadrature *quadBd = fem->quadBd; 471 PetscSection section; 472 PetscReal *v0, *n, *J, *invJ, *detJ; 473 PetscScalar *elemVec, *u; 474 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 475 PetscInt cellDof, numComponents; 476 PetscBool has; 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 481 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 482 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 483 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 484 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 485 numCells = cEnd - cStart; 486 for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) { 487 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 488 numComponents += quad[field].numComponents; 489 } 490 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 491 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 492 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); 493 for (c = cStart; c < cEnd; ++c) { 494 PetscScalar *x; 495 PetscInt i; 496 497 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 498 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 499 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 500 501 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 502 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 503 } 504 for (field = 0; field < numFields; ++field) { 505 const PetscInt numQuadPoints = quad[field].numQuadPoints; 506 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 507 void (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0Funcs[field]; 508 void (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1Funcs[field]; 509 /* Conforming batches */ 510 PetscInt blockSize = numBasisFuncs*numQuadPoints; 511 PetscInt numBlocks = 1; 512 PetscInt batchSize = numBlocks * blockSize; 513 PetscInt numBatches = numBatchesTmp; 514 PetscInt numChunks = numCells / (numBatches*batchSize); 515 /* Remainder */ 516 PetscInt numRemainder = numCells % (numBatches * batchSize); 517 PetscInt offset = numCells - numRemainder; 518 519 ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr); 520 ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 521 f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 522 } 523 for (c = cStart; c < cEnd; ++c) { 524 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 525 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 526 } 527 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 528 /* Integration over the boundary: 529 - This can probably be generalized to integration over a set of labels, however 530 the idea here is to do integration where we need the cell normal 531 - We can replace hardcoding with a registration process, and this is how we hook 532 up the system to something like FEniCS 533 */ 534 ierr = DMPlexHasLabel(dm, "boundary", &has);CHKERRQ(ierr); 535 if (has && quadBd) { 536 DMLabel label; 537 IS pointIS; 538 const PetscInt *points; 539 PetscInt numPoints, p; 540 541 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 542 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 543 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 544 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 545 for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) { 546 cellDof += quadBd[field].numBasisFuncs*quadBd[field].numComponents; 547 numComponents += quadBd[field].numComponents; 548 } 549 ierr = PetscMalloc7(numPoints*cellDof,PetscScalar,&u,numPoints*dim,PetscReal,&v0,numPoints*dim,PetscReal,&n,numPoints*dim*dim,PetscReal,&J,numPoints*dim*dim,PetscReal,&invJ,numPoints,PetscReal,&detJ,numPoints*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr); 550 for (p = 0; p < numPoints; ++p) { 551 const PetscInt point = points[p]; 552 PetscScalar *x; 553 PetscInt i; 554 555 /* TODO: Add normal determination here */ 556 ierr = DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);CHKERRQ(ierr); 557 if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[p], p); 558 ierr = DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr); 559 560 for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i]; 561 ierr = DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr); 562 } 563 for (field = 0; field < numFields; ++field) { 564 const PetscInt numQuadPoints = quadBd[field].numQuadPoints; 565 const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs; 566 void (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]) = fem->f0BdFuncs[field]; 567 void (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]) = fem->f1BdFuncs[field]; 568 /* Conforming batches */ 569 PetscInt blockSize = numBasisFuncs*numQuadPoints; 570 PetscInt numBlocks = 1; 571 PetscInt batchSize = numBlocks * blockSize; 572 PetscInt numBatches = numBatchesTmp; 573 PetscInt numChunks = numPoints / (numBatches*batchSize); 574 /* Remainder */ 575 PetscInt numRemainder = numPoints % (numBatches * batchSize); 576 PetscInt offset = numPoints - numRemainder; 577 578 ierr = (*mesh->integrateBdResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quadBd, u, v0, n, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr); 579 ierr = (*mesh->integrateBdResidualFEM)(numRemainder, numFields, field, quadBd, &u[offset*cellDof], &v0[offset*dim], &n[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 580 f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 581 } 582 for (p = 0; p < numPoints; ++p) { 583 const PetscInt point = points[p]; 584 585 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);CHKERRQ(ierr);} 586 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);CHKERRQ(ierr); 587 } 588 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 589 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 590 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 591 } 592 if (mesh->printFEM) { 593 PetscMPIInt rank, numProcs; 594 PetscInt p; 595 596 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 597 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 598 ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr); 599 for (p = 0; p < numProcs; ++p) { 600 if (p == rank) { 601 Vec f; 602 603 ierr = VecDuplicate(F, &f);CHKERRQ(ierr); 604 ierr = VecCopy(F, f);CHKERRQ(ierr); 605 ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr); 606 ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 607 ierr = VecDestroy(&f);CHKERRQ(ierr); 608 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 609 } 610 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 611 } 612 } 613 /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 614 PetscFunctionReturn(0); 615 } 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 619 /*@C 620 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 621 622 Input Parameters: 623 + dm - The mesh 624 . J - The Jacobian shell matrix 625 . X - Local input vector 626 - user - The user context 627 628 Output Parameter: 629 . F - Local output vector 630 631 Note: 632 The second member of the user context must be an FEMContext. 633 634 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 635 like a GPU, or vectorize on a multicore machine. 636 637 Level: developer 638 639 .seealso: DMPlexComputeResidualFEM() 640 @*/ 641 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 642 { 643 DM_Plex *mesh = (DM_Plex*) dm->data; 644 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 645 PetscQuadrature *quad = fem->quad; 646 PetscSection section; 647 JacActionCtx *jctx; 648 PetscReal *v0, *J, *invJ, *detJ; 649 PetscScalar *elemVec, *u, *a; 650 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 651 PetscInt cellDof = 0; 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 656 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 657 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 658 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 659 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 660 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 661 numCells = cEnd - cStart; 662 for (field = 0; field < numFields; ++field) { 663 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 664 } 665 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 666 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); 667 for (c = cStart; c < cEnd; ++c) { 668 PetscScalar *x; 669 PetscInt i; 670 671 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 672 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 673 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 674 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 675 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 676 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 677 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 678 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 679 } 680 for (field = 0; field < numFields; ++field) { 681 const PetscInt numQuadPoints = quad[field].numQuadPoints; 682 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 683 /* Conforming batches */ 684 PetscInt blockSize = numBasisFuncs*numQuadPoints; 685 PetscInt numBlocks = 1; 686 PetscInt batchSize = numBlocks * blockSize; 687 PetscInt numBatches = numBatchesTmp; 688 PetscInt numChunks = numCells / (numBatches*batchSize); 689 /* Remainder */ 690 PetscInt numRemainder = numCells % (numBatches * batchSize); 691 PetscInt offset = numCells - numRemainder; 692 693 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); 694 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], 695 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 696 } 697 for (c = cStart; c < cEnd; ++c) { 698 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 699 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 700 } 701 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 702 if (mesh->printFEM) { 703 PetscMPIInt rank, numProcs; 704 PetscInt p; 705 706 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 707 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 708 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action:\n");CHKERRQ(ierr); 709 for (p = 0; p < numProcs; ++p) { 710 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 711 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 712 } 713 } 714 /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 715 PetscFunctionReturn(0); 716 } 717 718 #undef __FUNCT__ 719 #define __FUNCT__ "DMPlexComputeJacobianFEM" 720 /*@ 721 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 722 723 Input Parameters: 724 + dm - The mesh 725 . X - Local input vector 726 - user - The user context 727 728 Output Parameter: 729 . Jac - Jacobian matrix 730 731 Note: 732 The second member of the user context must be an FEMContext. 733 734 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 735 like a GPU, or vectorize on a multicore machine. 736 737 Level: developer 738 739 .seealso: FormFunctionLocal() 740 @*/ 741 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user) 742 { 743 DM_Plex *mesh = (DM_Plex*) dm->data; 744 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 745 PetscQuadrature *quad = fem->quad; 746 PetscSection section; 747 PetscReal *v0, *J, *invJ, *detJ; 748 PetscScalar *elemMat, *u; 749 PetscInt dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c; 750 PetscInt cellDof = 0, numComponents = 0; 751 PetscBool isShell; 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 756 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 757 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 758 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 759 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 760 numCells = cEnd - cStart; 761 for (field = 0; field < numFields; ++field) { 762 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 763 numComponents += quad[field].numComponents; 764 } 765 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 766 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 767 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); 768 for (c = cStart; c < cEnd; ++c) { 769 PetscScalar *x; 770 PetscInt i; 771 772 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 773 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 774 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 775 776 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 777 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 778 } 779 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 780 for (fieldI = 0; fieldI < numFields; ++fieldI) { 781 const PetscInt numQuadPoints = quad[fieldI].numQuadPoints; 782 const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs; 783 PetscInt fieldJ; 784 785 for (fieldJ = 0; fieldJ < numFields; ++fieldJ) { 786 void (*g0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]) = fem->g0Funcs[fieldI*numFields+fieldJ]; 787 void (*g1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]) = fem->g1Funcs[fieldI*numFields+fieldJ]; 788 void (*g2)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]) = fem->g2Funcs[fieldI*numFields+fieldJ]; 789 void (*g3)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]) = fem->g3Funcs[fieldI*numFields+fieldJ]; 790 /* Conforming batches */ 791 PetscInt blockSize = numBasisFuncs*numQuadPoints; 792 PetscInt numBlocks = 1; 793 PetscInt batchSize = numBlocks * blockSize; 794 PetscInt numBatches = numBatchesTmp; 795 PetscInt numChunks = numCells / (numBatches*batchSize); 796 /* Remainder */ 797 PetscInt numRemainder = numCells % (numBatches * batchSize); 798 PetscInt offset = numCells - numRemainder; 799 800 ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 801 ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 802 g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr); 803 } 804 } 805 for (c = cStart; c < cEnd; ++c) { 806 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 807 ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 808 } 809 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 810 811 /* Assemble matrix, using the 2-step process: 812 MatAssemblyBegin(), MatAssemblyEnd(). */ 813 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 814 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 815 816 if (mesh->printFEM) { 817 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr); 818 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 819 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 820 } 821 /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 822 ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr); 823 if (isShell) { 824 JacActionCtx *jctx; 825 826 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 827 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 828 } 829 *str = SAME_NONZERO_PATTERN; 830 PetscFunctionReturn(0); 831 } 832