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 PetscQuadrature *quadBd = fem->quadBd; 466 PetscSection section; 467 PetscReal *v0, *J, *invJ, *detJ; 468 PetscScalar *elemVec, *u; 469 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 470 PetscInt cellDof, numComponents; 471 PetscBool has; 472 PetscErrorCode ierr; 473 474 PetscFunctionBegin; 475 /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 476 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 477 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 478 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 479 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 480 numCells = cEnd - cStart; 481 for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) { 482 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 483 numComponents += quad[field].numComponents; 484 } 485 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 486 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 487 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); 488 for (c = cStart; c < cEnd; ++c) { 489 PetscScalar *x; 490 PetscInt i; 491 492 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 493 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 494 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 495 496 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 497 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 498 } 499 for (field = 0; field < numFields; ++field) { 500 const PetscInt numQuadPoints = quad[field].numQuadPoints; 501 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 502 void (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0Funcs[field]; 503 void (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1Funcs[field]; 504 /* Conforming batches */ 505 PetscInt blockSize = numBasisFuncs*numQuadPoints; 506 PetscInt numBlocks = 1; 507 PetscInt batchSize = numBlocks * blockSize; 508 PetscInt numBatches = numBatchesTmp; 509 PetscInt numChunks = numCells / (numBatches*batchSize); 510 /* Remainder */ 511 PetscInt numRemainder = numCells % (numBatches * batchSize); 512 PetscInt offset = numCells - numRemainder; 513 514 ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr); 515 ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 516 f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 517 } 518 for (c = cStart; c < cEnd; ++c) { 519 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 520 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 521 } 522 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 523 /* Integration over the boundary: 524 - This can probably be generalized to integration over a set of labels, however 525 the idea here is to do integration where we need the cell normal 526 - We can replace hardcoding with a registration process, and this is how we hook 527 up the system to something like FEniCS 528 */ 529 ierr = DMPlexHasLabel(dm, "boundary", &has);CHKERRQ(ierr); 530 if (has && quadBd) { 531 DMLabel label; 532 IS pointIS; 533 const PetscInt *points; 534 PetscInt numPoints, p; 535 536 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 537 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 538 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 539 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 540 ierr = DMPlexGetLabelValue(dm, "depth", points[0], &dim);CHKERRQ(ierr); 541 for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) { 542 cellDof += quadBd[field].numBasisFuncs*quadBd[field].numComponents; 543 numComponents += quadBd[field].numComponents; 544 } 545 ierr = PetscMalloc6(numPoints*cellDof,PetscScalar,&u,numPoints*dim,PetscReal,&v0,numPoints*dim*dim,PetscReal,&J,numPoints*dim*dim,PetscReal,&invJ,numPoints,PetscReal,&detJ,numPoints*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr); 546 for (p = 0; p < numPoints; ++p) { 547 const PetscInt point = points[p]; 548 PetscScalar *x; 549 PetscInt i; 550 551 ierr = DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);CHKERRQ(ierr); 552 if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[p], p); 553 ierr = DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr); 554 555 for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i]; 556 ierr = DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr); 557 } 558 for (field = 0; field < numFields; ++field) { 559 const PetscInt numQuadPoints = quadBd[field].numQuadPoints; 560 const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs; 561 void (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0BdFuncs[field]; 562 void (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1BdFuncs[field]; 563 /* Conforming batches */ 564 PetscInt blockSize = numBasisFuncs*numQuadPoints; 565 PetscInt numBlocks = 1; 566 PetscInt batchSize = numBlocks * blockSize; 567 PetscInt numBatches = numBatchesTmp; 568 PetscInt numChunks = numPoints / (numBatches*batchSize); 569 /* Remainder */ 570 PetscInt numRemainder = numPoints % (numBatches * batchSize); 571 PetscInt offset = numPoints - numRemainder; 572 573 ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr); 574 ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 575 f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 576 } 577 for (p = 0; p < numPoints; ++p) { 578 const PetscInt point = points[p]; 579 580 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);CHKERRQ(ierr);} 581 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);CHKERRQ(ierr); 582 } 583 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 584 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 585 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 586 } 587 if (mesh->printFEM) { 588 PetscMPIInt rank, numProcs; 589 PetscInt p; 590 591 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 592 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 593 ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr); 594 for (p = 0; p < numProcs; ++p) { 595 if (p == rank) { 596 Vec f; 597 598 ierr = VecDuplicate(F, &f);CHKERRQ(ierr); 599 ierr = VecCopy(F, f);CHKERRQ(ierr); 600 ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr); 601 ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 602 ierr = VecDestroy(&f);CHKERRQ(ierr); 603 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 604 } 605 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 606 } 607 } 608 /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 609 PetscFunctionReturn(0); 610 } 611 612 #undef __FUNCT__ 613 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 614 /*@C 615 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 616 617 Input Parameters: 618 + dm - The mesh 619 . J - The Jacobian shell matrix 620 . X - Local input vector 621 - user - The user context 622 623 Output Parameter: 624 . F - Local output vector 625 626 Note: 627 The second member of the user context must be an FEMContext. 628 629 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 630 like a GPU, or vectorize on a multicore machine. 631 632 Level: developer 633 634 .seealso: DMPlexComputeResidualFEM() 635 @*/ 636 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 637 { 638 DM_Plex *mesh = (DM_Plex*) dm->data; 639 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 640 PetscQuadrature *quad = fem->quad; 641 PetscSection section; 642 JacActionCtx *jctx; 643 PetscReal *v0, *J, *invJ, *detJ; 644 PetscScalar *elemVec, *u, *a; 645 PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c; 646 PetscInt cellDof = 0; 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 651 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 652 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 653 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 654 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 655 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 656 numCells = cEnd - cStart; 657 for (field = 0; field < numFields; ++field) { 658 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 659 } 660 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 661 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); 662 for (c = cStart; c < cEnd; ++c) { 663 PetscScalar *x; 664 PetscInt i; 665 666 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 667 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 668 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 669 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 670 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 671 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 672 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 673 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 674 } 675 for (field = 0; field < numFields; ++field) { 676 const PetscInt numQuadPoints = quad[field].numQuadPoints; 677 const PetscInt numBasisFuncs = quad[field].numBasisFuncs; 678 /* Conforming batches */ 679 PetscInt blockSize = numBasisFuncs*numQuadPoints; 680 PetscInt numBlocks = 1; 681 PetscInt batchSize = numBlocks * blockSize; 682 PetscInt numBatches = numBatchesTmp; 683 PetscInt numChunks = numCells / (numBatches*batchSize); 684 /* Remainder */ 685 PetscInt numRemainder = numCells % (numBatches * batchSize); 686 PetscInt offset = numCells - numRemainder; 687 688 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); 689 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], 690 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 691 } 692 for (c = cStart; c < cEnd; ++c) { 693 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 694 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 695 } 696 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 697 if (mesh->printFEM) { 698 PetscMPIInt rank, numProcs; 699 PetscInt p; 700 701 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 702 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 703 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action:\n");CHKERRQ(ierr); 704 for (p = 0; p < numProcs; ++p) { 705 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 706 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 707 } 708 } 709 /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "DMPlexComputeJacobianFEM" 715 /*@ 716 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 717 718 Input Parameters: 719 + dm - The mesh 720 . X - Local input vector 721 - user - The user context 722 723 Output Parameter: 724 . Jac - Jacobian matrix 725 726 Note: 727 The second member of the user context must be an FEMContext. 728 729 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 730 like a GPU, or vectorize on a multicore machine. 731 732 Level: developer 733 734 .seealso: FormFunctionLocal() 735 @*/ 736 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user) 737 { 738 DM_Plex *mesh = (DM_Plex*) dm->data; 739 PetscFEM *fem = (PetscFEM*) &((DM*) user)[1]; 740 PetscQuadrature *quad = fem->quad; 741 PetscSection section; 742 PetscReal *v0, *J, *invJ, *detJ; 743 PetscScalar *elemMat, *u; 744 PetscInt dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c; 745 PetscInt cellDof = 0, numComponents = 0; 746 PetscBool isShell; 747 PetscErrorCode ierr; 748 749 PetscFunctionBegin; 750 /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 751 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 752 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 753 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 754 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 755 numCells = cEnd - cStart; 756 for (field = 0; field < numFields; ++field) { 757 cellDof += quad[field].numBasisFuncs*quad[field].numComponents; 758 numComponents += quad[field].numComponents; 759 } 760 ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr); 761 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 762 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); 763 for (c = cStart; c < cEnd; ++c) { 764 PetscScalar *x; 765 PetscInt i; 766 767 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 768 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 769 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 770 771 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 772 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 773 } 774 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 775 for (fieldI = 0; fieldI < numFields; ++fieldI) { 776 const PetscInt numQuadPoints = quad[fieldI].numQuadPoints; 777 const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs; 778 PetscInt fieldJ; 779 780 for (fieldJ = 0; fieldJ < numFields; ++fieldJ) { 781 void (*g0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]) = fem->g0Funcs[fieldI*numFields+fieldJ]; 782 void (*g1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]) = fem->g1Funcs[fieldI*numFields+fieldJ]; 783 void (*g2)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]) = fem->g2Funcs[fieldI*numFields+fieldJ]; 784 void (*g3)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]) = fem->g3Funcs[fieldI*numFields+fieldJ]; 785 /* Conforming batches */ 786 PetscInt blockSize = numBasisFuncs*numQuadPoints; 787 PetscInt numBlocks = 1; 788 PetscInt batchSize = numBlocks * blockSize; 789 PetscInt numBatches = numBatchesTmp; 790 PetscInt numChunks = numCells / (numBatches*batchSize); 791 /* Remainder */ 792 PetscInt numRemainder = numCells % (numBatches * batchSize); 793 PetscInt offset = numCells - numRemainder; 794 795 ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 796 ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset], 797 g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr); 798 } 799 } 800 for (c = cStart; c < cEnd; ++c) { 801 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 802 ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 803 } 804 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 805 806 /* Assemble matrix, using the 2-step process: 807 MatAssemblyBegin(), MatAssemblyEnd(). */ 808 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 809 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 810 811 if (mesh->printFEM) { 812 ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr); 813 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 814 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 815 } 816 /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */ 817 ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr); 818 if (isShell) { 819 JacActionCtx *jctx; 820 821 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 822 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 823 } 824 *str = SAME_NONZERO_PATTERN; 825 PetscFunctionReturn(0); 826 } 827