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