1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscblaslapack.h> 3 4 PetscErrorCode PetscFEDestroy_Basic(PetscFE fem) 5 { 6 PetscFE_Basic *b = (PetscFE_Basic *) fem->data; 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = PetscFree(b);CHKERRQ(ierr); 11 PetscFunctionReturn(0); 12 } 13 14 PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer) 15 { 16 PetscSpace basis; 17 PetscDualSpace dual; 18 PetscQuadrature q = NULL; 19 PetscInt dim, Nc, Nq; 20 PetscViewerFormat format; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 ierr = PetscFEGetBasisSpace(fe, &basis);CHKERRQ(ierr); 25 ierr = PetscFEGetDualSpace(fe, &dual);CHKERRQ(ierr); 26 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 27 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 28 if (q) {ierr = PetscQuadratureGetData(q, &dim, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);} 29 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 30 ierr = PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");CHKERRQ(ierr); 31 ierr = PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);CHKERRQ(ierr); 32 ierr = PetscViewerASCIIPrintf(viewer, " components: %d\n", Nc);CHKERRQ(ierr); 33 ierr = PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);CHKERRQ(ierr); 34 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 35 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 36 if (q) {ierr = PetscQuadratureView(q, viewer);CHKERRQ(ierr);} 37 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 38 } 39 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 40 ierr = PetscSpaceView(basis, viewer);CHKERRQ(ierr); 41 ierr = PetscDualSpaceView(dual, viewer);CHKERRQ(ierr); 42 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer) 47 { 48 PetscBool iascii; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1); 53 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 54 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 55 if (iascii) {ierr = PetscFEView_Basic_Ascii(fe, viewer);CHKERRQ(ierr);} 56 PetscFunctionReturn(0); 57 } 58 59 /* Construct the change of basis from prime basis to nodal basis */ 60 PetscErrorCode PetscFESetUp_Basic(PetscFE fem) 61 { 62 PetscScalar *work, *invVscalar; 63 PetscBLASInt *pivots; 64 PetscBLASInt n, info; 65 PetscInt pdim, j; 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 70 ierr = PetscMalloc1(pdim*pdim,&fem->invV);CHKERRQ(ierr); 71 #if defined(PETSC_USE_COMPLEX) 72 ierr = PetscMalloc1(pdim*pdim,&invVscalar);CHKERRQ(ierr); 73 #else 74 invVscalar = fem->invV; 75 #endif 76 for (j = 0; j < pdim; ++j) { 77 PetscReal *Bf; 78 PetscQuadrature f; 79 const PetscReal *points, *weights; 80 PetscInt Nc, Nq, q, k, c; 81 82 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr); 83 ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr); 84 ierr = PetscMalloc1(Nc*Nq*pdim,&Bf);CHKERRQ(ierr); 85 ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr); 86 for (k = 0; k < pdim; ++k) { 87 /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */ 88 invVscalar[j*pdim+k] = 0.0; 89 90 for (q = 0; q < Nq; ++q) { 91 for (c = 0; c < Nc; ++c) invVscalar[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c]; 92 } 93 } 94 ierr = PetscFree(Bf);CHKERRQ(ierr); 95 } 96 ierr = PetscMalloc2(pdim,&pivots,pdim,&work);CHKERRQ(ierr); 97 n = pdim; 98 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info)); 99 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info)); 100 #if defined(PETSC_USE_COMPLEX) 101 for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]); 102 ierr = PetscFree(invVscalar);CHKERRQ(ierr); 103 #endif 104 ierr = PetscFree2(pivots,work);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim) 109 { 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 ierr = PetscDualSpaceGetDimension(fem->dualSpace, dim);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H) 118 { 119 DM dm; 120 PetscInt pdim; /* Dimension of FE space P */ 121 PetscInt dim; /* Spatial dimension */ 122 PetscInt Nc; /* Field components */ 123 PetscReal *tmpB, *tmpD, *tmpH; 124 PetscInt p, d, j, k, c; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 129 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 130 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 131 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 132 /* Evaluate the prime basis functions at all points */ 133 if (B) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 134 if (D) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 135 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 136 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);CHKERRQ(ierr); 137 /* Translate to the nodal basis */ 138 for (p = 0; p < npoints; ++p) { 139 if (B) { 140 /* Multiply by V^{-1} (pdim x pdim) */ 141 for (j = 0; j < pdim; ++j) { 142 const PetscInt i = (p*pdim + j)*Nc; 143 144 for (c = 0; c < Nc; ++c) { 145 B[i+c] = 0.0; 146 for (k = 0; k < pdim; ++k) { 147 B[i+c] += fem->invV[k*pdim+j] * tmpB[(p*pdim + k)*Nc+c]; 148 } 149 } 150 } 151 } 152 if (D) { 153 /* Multiply by V^{-1} (pdim x pdim) */ 154 for (j = 0; j < pdim; ++j) { 155 for (c = 0; c < Nc; ++c) { 156 for (d = 0; d < dim; ++d) { 157 const PetscInt i = ((p*pdim + j)*Nc + c)*dim + d; 158 159 D[i] = 0.0; 160 for (k = 0; k < pdim; ++k) { 161 D[i] += fem->invV[k*pdim+j] * tmpD[((p*pdim + k)*Nc + c)*dim + d]; 162 } 163 } 164 } 165 } 166 } 167 if (H) { 168 /* Multiply by V^{-1} (pdim x pdim) */ 169 for (j = 0; j < pdim; ++j) { 170 for (c = 0; c < Nc; ++c) { 171 for (d = 0; d < dim*dim; ++d) { 172 const PetscInt i = ((p*pdim + j)*Nc + c)*dim*dim + d; 173 174 H[i] = 0.0; 175 for (k = 0; k < pdim; ++k) { 176 H[i] += fem->invV[k*pdim+j] * tmpH[((p*pdim + k)*Nc + c)*dim*dim + d]; 177 } 178 } 179 } 180 } 181 } 182 } 183 if (B) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 184 if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 185 if (H) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 186 PetscFunctionReturn(0); 187 } 188 189 PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 190 const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 191 { 192 const PetscInt debug = 0; 193 PetscPointFunc obj_func; 194 PetscQuadrature quad; 195 PetscScalar *u, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 196 const PetscScalar *constants; 197 PetscReal *x; 198 PetscReal **B, **D, **BAux = NULL, **DAux = NULL; 199 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 200 PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 201 PetscBool isAffine; 202 const PetscReal *quadPoints, *quadWeights; 203 PetscInt qNc, Nq, q; 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 ierr = PetscDSGetObjective(prob, field, &obj_func);CHKERRQ(ierr); 208 if (!obj_func) PetscFunctionReturn(0); 209 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 210 ierr = PetscFEGetQuadrature(fem, &quad);CHKERRQ(ierr); 211 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 212 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 213 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 214 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 215 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 216 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 217 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 218 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 219 ierr = PetscDSGetTabulation(prob, &B, &D);CHKERRQ(ierr); 220 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 221 if (probAux) { 222 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 223 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 224 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 225 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 226 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 227 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 228 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 229 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 230 ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 231 } 232 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 233 Np = cgeom->numPoints; 234 dE = cgeom->dimEmbed; 235 isAffine = cgeom->isAffine; 236 for (e = 0; e < Ne; ++e) { 237 const PetscReal *v0 = &cgeom->v[e*Np*dE]; 238 const PetscReal *J = &cgeom->J[e*Np*dE*dE]; 239 240 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 241 for (q = 0; q < Nq; ++q) { 242 PetscScalar integrand; 243 const PetscReal *v; 244 const PetscReal *invJ; 245 PetscReal detJ; 246 247 if (isAffine) { 248 CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x); 249 v = x; 250 invJ = &cgeom->invJ[e*dE*dE]; 251 detJ = cgeom->detJ[e]; 252 } else { 253 v = &v0[q*dE]; 254 invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 255 detJ = cgeom->detJ[e*Np + q]; 256 } 257 if (debug > 1 && q < Np) { 258 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 259 #if !defined(PETSC_USE_COMPLEX) 260 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 261 #endif 262 } 263 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 264 EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], NULL, u, u_x, NULL); 265 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 266 obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, v, numConstants, constants, &integrand); 267 integrand *= detJ*quadWeights[q]; 268 integral[e*Nf+field] += integrand; 269 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));CHKERRQ(ierr);} 270 } 271 cOffset += totDim; 272 cOffsetAux += totDimAux; 273 } 274 PetscFunctionReturn(0); 275 } 276 277 PetscErrorCode PetscFEIntegrateBd_Basic(PetscFE fem, PetscDS prob, PetscInt field, 278 PetscBdPointFunc obj_func, 279 PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 280 { 281 const PetscInt debug = 0; 282 PetscQuadrature quad; 283 PetscScalar *u, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 284 const PetscScalar *constants; 285 PetscReal *x; 286 PetscReal **B, **D, **BAux = NULL, **DAux = NULL; 287 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 288 PetscBool isAffine, auxOnBd; 289 const PetscReal *quadPoints, *quadWeights; 290 PetscInt qNc, Nq, q, Np, dE; 291 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 292 PetscErrorCode ierr; 293 294 PetscFunctionBegin; 295 if (!obj_func) PetscFunctionReturn(0); 296 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 297 ierr = PetscFEGetFaceQuadrature(fem, &quad);CHKERRQ(ierr); 298 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 299 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 300 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 301 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 302 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 303 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 304 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 305 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 306 ierr = PetscDSGetFaceTabulation(prob, &B, &D);CHKERRQ(ierr); 307 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 308 if (probAux) { 309 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 310 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 311 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 312 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 313 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 314 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 315 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 316 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 317 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 318 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 319 if (auxOnBd) {ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr);} 320 else {ierr = PetscDSGetFaceTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr);} 321 } 322 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 323 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 324 Np = fgeom->numPoints; 325 dE = fgeom->dimEmbed; 326 isAffine = fgeom->isAffine; 327 for (e = 0; e < Ne; ++e) { 328 const PetscReal *v0 = &fgeom->v[e*Np*dE]; 329 const PetscReal *J = &fgeom->J[e*Np*dE*dE]; 330 const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */ 331 332 for (q = 0; q < Nq; ++q) { 333 const PetscReal *v; 334 const PetscReal *invJ; 335 const PetscReal *n; 336 PetscReal detJ; 337 PetscScalar integrand; 338 339 if (isAffine) { 340 CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x); 341 v = x; 342 invJ = &fgeom->suppInvJ[0][e*dE*dE]; 343 detJ = fgeom->detJ[e]; 344 n = &fgeom->n[e*dE]; 345 } else { 346 v = &v0[q*dE]; 347 invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 348 detJ = fgeom->detJ[e*Np + q]; 349 n = &fgeom->n[(e*Np+q)*dE]; 350 } 351 if (debug > 1 && q < Np) { 352 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 353 #ifndef PETSC_USE_COMPLEX 354 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 355 #endif 356 } 357 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 358 EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], NULL, u, u_x, NULL); 359 if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 360 obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, v, n, numConstants, constants, &integrand); 361 integrand *= detJ*quadWeights[q]; 362 integral[e*Nf+field] += integrand; 363 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);} 364 } 365 cOffset += totDim; 366 cOffsetAux += totDimAux; 367 } 368 PetscFunctionReturn(0); 369 } 370 371 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 372 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 373 { 374 const PetscInt debug = 0; 375 PetscPointFunc f0_func; 376 PetscPointFunc f1_func; 377 PetscQuadrature quad; 378 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 379 const PetscScalar *constants; 380 PetscReal *x; 381 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 382 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 383 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 384 PetscInt dE, Np; 385 PetscBool isAffine; 386 const PetscReal *quadPoints, *quadWeights; 387 PetscInt qNc, Nq, q; 388 PetscErrorCode ierr; 389 390 PetscFunctionBegin; 391 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 392 ierr = PetscFEGetQuadrature(fem, &quad);CHKERRQ(ierr); 393 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 394 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 395 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 396 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 397 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 398 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 399 ierr = PetscDSGetFieldOffset(prob, field, &fOffset);CHKERRQ(ierr); 400 ierr = PetscDSGetResidual(prob, field, &f0_func, &f1_func);CHKERRQ(ierr); 401 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 402 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 403 ierr = PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 404 ierr = PetscDSGetTabulation(prob, &B, &D);CHKERRQ(ierr); 405 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 406 if (probAux) { 407 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 408 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 409 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 410 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 411 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 412 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 413 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 414 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 415 ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 416 } 417 NbI = Nb[field]; 418 NcI = Nc[field]; 419 BI = B[field]; 420 DI = D[field]; 421 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 422 Np = cgeom->numPoints; 423 dE = cgeom->dimEmbed; 424 isAffine = cgeom->isAffine; 425 for (e = 0; e < Ne; ++e) { 426 const PetscReal *v0 = &cgeom->v[e*Np*dE]; 427 const PetscReal *J = &cgeom->J[e*Np*dE*dE]; 428 429 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 430 ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr); 431 ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 432 for (q = 0; q < Nq; ++q) { 433 const PetscReal *v; 434 const PetscReal *invJ; 435 PetscReal detJ; 436 437 if (isAffine) { 438 CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x); 439 v = x; 440 invJ = &cgeom->invJ[e*dE*dE]; 441 detJ = cgeom->detJ[e]; 442 } else { 443 v = &v0[q*dE]; 444 invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 445 detJ = cgeom->detJ[e*Np + q]; 446 } 447 if (debug > 1 && q < Np) { 448 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 449 #if !defined(PETSC_USE_COMPLEX) 450 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 451 #endif 452 } 453 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 454 EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 455 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 456 if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, numConstants, constants, &f0[q*NcI]); 457 if (f1_func) { 458 ierr = PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 459 f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, numConstants, constants, refSpaceDer); 460 } 461 TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL); 462 } 463 UpdateElementVec(dim, Nq, NbI, NcI, BI, DI, f0, f1, &elemVec[cOffset+fOffset]); 464 cOffset += totDim; 465 cOffsetAux += totDimAux; 466 } 467 PetscFunctionReturn(0); 468 } 469 470 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 471 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 472 { 473 const PetscInt debug = 0; 474 PetscBdPointFunc f0_func; 475 PetscBdPointFunc f1_func; 476 PetscQuadrature quad; 477 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 478 const PetscScalar *constants; 479 PetscReal *x; 480 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 481 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 482 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 483 PetscBool isAffine, auxOnBd; 484 const PetscReal *quadPoints, *quadWeights; 485 PetscInt qNc, Nq, q, Np, dE; 486 PetscErrorCode ierr; 487 488 PetscFunctionBegin; 489 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 490 ierr = PetscFEGetFaceQuadrature(fem, &quad);CHKERRQ(ierr); 491 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 492 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 493 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 494 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 495 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 496 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 497 ierr = PetscDSGetFieldOffset(prob, field, &fOffset);CHKERRQ(ierr); 498 ierr = PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);CHKERRQ(ierr); 499 if (!f0_func && !f1_func) PetscFunctionReturn(0); 500 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 501 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 502 ierr = PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 503 ierr = PetscDSGetFaceTabulation(prob, &B, &D);CHKERRQ(ierr); 504 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 505 if (probAux) { 506 ierr = PetscDSGetSpatialDimension(probAux, &dimAux);CHKERRQ(ierr); 507 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 508 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 509 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 510 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 511 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 512 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 513 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 514 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 515 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 516 if (auxOnBd) {ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr);} 517 else {ierr = PetscDSGetFaceTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr);} 518 } 519 NbI = Nb[field]; 520 NcI = Nc[field]; 521 BI = B[field]; 522 DI = D[field]; 523 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 524 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 525 Np = fgeom->numPoints; 526 dE = fgeom->dimEmbed; 527 isAffine = fgeom->isAffine; 528 for (e = 0; e < Ne; ++e) { 529 const PetscReal *v0 = &fgeom->v[e*Np*dE]; 530 const PetscReal *J = &fgeom->J[e*Np*dE*dE]; 531 const PetscInt face = fgeom->face[e][0]; 532 533 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 534 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 535 ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr); 536 ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 537 for (q = 0; q < Nq; ++q) { 538 const PetscReal *v; 539 const PetscReal *invJ; 540 const PetscReal *n; 541 PetscReal detJ; 542 if (isAffine) { 543 CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x); 544 v = x; 545 invJ = &fgeom->suppInvJ[0][e*dE*dE]; 546 detJ = fgeom->detJ[e]; 547 n = &fgeom->n[e*dE]; 548 } else { 549 v = &v0[q*dE]; 550 invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 551 detJ = fgeom->detJ[e*Np + q]; 552 n = &fgeom->n[(e*Np+q)*dE]; 553 } 554 if (debug > 1 && q < Np) { 555 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 556 #if !defined(PETSC_USE_COMPLEX) 557 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 558 #endif 559 } 560 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 561 EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 562 if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 563 if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, n, numConstants, constants, &f0[q*NcI]); 564 if (f1_func) { 565 ierr = PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 566 f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, n, numConstants, constants, refSpaceDer); 567 } 568 TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL); 569 } 570 UpdateElementVec(dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], f0, f1, &elemVec[cOffset+fOffset]); 571 cOffset += totDim; 572 cOffsetAux += totDimAux; 573 } 574 PetscFunctionReturn(0); 575 } 576 577 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *geom, 578 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 579 { 580 const PetscInt debug = 0; 581 PetscPointJac g0_func; 582 PetscPointJac g1_func; 583 PetscPointJac g2_func; 584 PetscPointJac g3_func; 585 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 586 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 587 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 588 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 589 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 590 PetscQuadrature quad; 591 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 592 const PetscScalar *constants; 593 PetscReal *x; 594 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 595 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 596 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 597 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 598 PetscInt dE, Np; 599 PetscBool isAffine; 600 const PetscReal *quadPoints, *quadWeights; 601 PetscInt qNc, Nq, q; 602 PetscErrorCode ierr; 603 604 PetscFunctionBegin; 605 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 606 ierr = PetscFEGetQuadrature(fem, &quad);CHKERRQ(ierr); 607 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 608 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 609 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 610 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 611 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 612 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 613 switch(jtype) { 614 case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 615 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 616 case PETSCFE_JACOBIAN: ierr = PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 617 } 618 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0); 619 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 620 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 621 ierr = PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 622 ierr = PetscDSGetTabulation(prob, &B, &D);CHKERRQ(ierr); 623 ierr = PetscDSGetFieldOffset(prob, fieldI, &offsetI);CHKERRQ(ierr); 624 ierr = PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);CHKERRQ(ierr); 625 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 626 if (probAux) { 627 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 628 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 629 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 630 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 631 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 632 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 633 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 634 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 635 ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 636 } 637 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 638 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 639 BI = B[fieldI], BJ = B[fieldJ]; 640 DI = D[fieldI], DJ = D[fieldJ]; 641 /* Initialize here in case the function is not defined */ 642 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 643 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 644 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 645 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 646 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 647 Np = geom->numPoints; 648 dE = geom->dimEmbed; 649 isAffine = geom->isAffine; 650 for (e = 0; e < Ne; ++e) { 651 const PetscReal *v0 = &geom->v[e*Np*dE]; 652 const PetscReal *J = &geom->J[e*Np*dE*dE]; 653 654 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 655 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 656 for (q = 0; q < Nq; ++q) { 657 const PetscReal *v; 658 const PetscReal *invJ; 659 PetscReal detJ; 660 const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ]; 661 const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim]; 662 PetscReal w; 663 PetscInt f, g, fc, gc, c; 664 665 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 666 if (isAffine) { 667 CoordinatesRefToReal(dE, dim, geom->xi, v0, J, &quadPoints[q*dim], x); 668 v = x; 669 invJ = &geom->invJ[e*dE*dE]; 670 detJ = geom->detJ[e]; 671 } else { 672 v = &v0[q*dE]; 673 invJ = &geom->invJ[(e*Np+q)*dE*dE]; 674 detJ = geom->detJ[e*Np + q]; 675 } 676 w = detJ*quadWeights[q]; 677 if (coefficients) EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 678 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 679 if (g0_func) { 680 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 681 g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, g0); 682 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 683 } 684 if (g1_func) { 685 PetscInt d, d2; 686 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 687 g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer); 688 for (fc = 0; fc < NcI; ++fc) { 689 for (gc = 0; gc < NcJ; ++gc) { 690 for (d = 0; d < dim; ++d) { 691 g1[(fc*NcJ+gc)*dim+d] = 0.0; 692 for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 693 g1[(fc*NcJ+gc)*dim+d] *= w; 694 } 695 } 696 } 697 } 698 if (g2_func) { 699 PetscInt d, d2; 700 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 701 g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer); 702 for (fc = 0; fc < NcI; ++fc) { 703 for (gc = 0; gc < NcJ; ++gc) { 704 for (d = 0; d < dim; ++d) { 705 g2[(fc*NcJ+gc)*dim+d] = 0.0; 706 for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 707 g2[(fc*NcJ+gc)*dim+d] *= w; 708 } 709 } 710 } 711 } 712 if (g3_func) { 713 PetscInt d, d2, dp, d3; 714 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 715 g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer); 716 for (fc = 0; fc < NcI; ++fc) { 717 for (gc = 0; gc < NcJ; ++gc) { 718 for (d = 0; d < dim; ++d) { 719 for (dp = 0; dp < dim; ++dp) { 720 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0; 721 for (d2 = 0; d2 < dim; ++d2) { 722 for (d3 = 0; d3 < dim; ++d3) { 723 g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3]; 724 } 725 } 726 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w; 727 } 728 } 729 } 730 } 731 } 732 733 for (f = 0; f < NbI; ++f) { 734 for (fc = 0; fc < NcI; ++fc) { 735 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 736 const PetscInt i = offsetI+f; /* Element matrix row */ 737 for (g = 0; g < NbJ; ++g) { 738 for (gc = 0; gc < NcJ; ++gc) { 739 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 740 const PetscInt j = offsetJ+g; /* Element matrix column */ 741 const PetscInt fOff = eOffset+i*totDim+j; 742 PetscInt d, d2; 743 744 elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx]; 745 for (d = 0; d < dim; ++d) { 746 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d]; 747 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx]; 748 for (d2 = 0; d2 < dim; ++d2) { 749 elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2]; 750 } 751 } 752 } 753 } 754 } 755 } 756 } 757 if (debug > 1) { 758 PetscInt fc, f, gc, g; 759 760 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 761 for (fc = 0; fc < NcI; ++fc) { 762 for (f = 0; f < NbI; ++f) { 763 const PetscInt i = offsetI + f*NcI+fc; 764 for (gc = 0; gc < NcJ; ++gc) { 765 for (g = 0; g < NbJ; ++g) { 766 const PetscInt j = offsetJ + g*NcJ+gc; 767 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 768 } 769 } 770 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 771 } 772 } 773 } 774 cOffset += totDim; 775 cOffsetAux += totDimAux; 776 eOffset += PetscSqr(totDim); 777 } 778 PetscFunctionReturn(0); 779 } 780 781 PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 782 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 783 { 784 const PetscInt debug = 0; 785 PetscBdPointJac g0_func; 786 PetscBdPointJac g1_func; 787 PetscBdPointJac g2_func; 788 PetscBdPointJac g3_func; 789 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 790 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 791 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 792 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 793 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 794 PetscQuadrature quad; 795 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 796 const PetscScalar *constants; 797 PetscReal *x; 798 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 799 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 800 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 801 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 802 PetscBool isAffine; 803 const PetscReal *quadPoints, *quadWeights; 804 PetscInt qNc, Nq, q, Np, dE; 805 PetscErrorCode ierr; 806 807 PetscFunctionBegin; 808 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 809 ierr = PetscFEGetFaceQuadrature(fem, &quad);CHKERRQ(ierr); 810 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 811 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 812 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 813 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 814 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 815 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 816 ierr = PetscDSGetFieldOffset(prob, fieldI, &offsetI);CHKERRQ(ierr); 817 ierr = PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);CHKERRQ(ierr); 818 ierr = PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr); 819 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 820 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 821 ierr = PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 822 ierr = PetscDSGetFaceTabulation(prob, &B, &D);CHKERRQ(ierr); 823 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 824 if (probAux) { 825 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 826 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 827 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 828 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 829 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 830 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 831 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 832 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 833 ierr = PetscDSGetFaceTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 834 } 835 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 836 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 837 BI = B[fieldI], BJ = B[fieldJ]; 838 DI = D[fieldI], DJ = D[fieldJ]; 839 /* Initialize here in case the function is not defined */ 840 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 841 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 842 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 843 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 844 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 845 Np = fgeom->numPoints; 846 dE = fgeom->dimEmbed; 847 isAffine = fgeom->isAffine; 848 for (e = 0; e < Ne; ++e) { 849 const PetscReal *v0 = &fgeom->v[e*Np*dE]; 850 const PetscReal *J = &fgeom->J[e*Np*dE*dE]; 851 const PetscInt face = fgeom->face[e][0]; 852 853 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 854 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 855 for (q = 0; q < Nq; ++q) { 856 const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ]; 857 const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim]; 858 PetscReal w; 859 PetscInt f, g, fc, gc, c; 860 const PetscReal *v; 861 const PetscReal *invJ; 862 const PetscReal *n; 863 PetscReal detJ; 864 865 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 866 if (isAffine) { 867 CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x); 868 v = x; 869 invJ = &fgeom->suppInvJ[0][e*dE*dE]; 870 detJ = fgeom->detJ[e]; 871 n = &fgeom->n[e*dE]; 872 } else { 873 v = &v0[q*dE]; 874 invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 875 detJ = fgeom->detJ[e*Np + q]; 876 n = &fgeom->n[(e*Np+q)*dE]; 877 } 878 w = detJ*quadWeights[q]; 879 880 if (coefficients) EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 881 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 882 if (g0_func) { 883 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 884 g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, g0); 885 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 886 } 887 if (g1_func) { 888 PetscInt d, d2; 889 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 890 g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer); 891 for (fc = 0; fc < NcI; ++fc) { 892 for (gc = 0; gc < NcJ; ++gc) { 893 for (d = 0; d < dim; ++d) { 894 g1[(fc*NcJ+gc)*dim+d] = 0.0; 895 for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 896 g1[(fc*NcJ+gc)*dim+d] *= w; 897 } 898 } 899 } 900 } 901 if (g2_func) { 902 PetscInt d, d2; 903 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 904 g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer); 905 for (fc = 0; fc < NcI; ++fc) { 906 for (gc = 0; gc < NcJ; ++gc) { 907 for (d = 0; d < dim; ++d) { 908 g2[(fc*NcJ+gc)*dim+d] = 0.0; 909 for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 910 g2[(fc*NcJ+gc)*dim+d] *= w; 911 } 912 } 913 } 914 } 915 if (g3_func) { 916 PetscInt d, d2, dp, d3; 917 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 918 g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer); 919 for (fc = 0; fc < NcI; ++fc) { 920 for (gc = 0; gc < NcJ; ++gc) { 921 for (d = 0; d < dim; ++d) { 922 for (dp = 0; dp < dim; ++dp) { 923 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0; 924 for (d2 = 0; d2 < dim; ++d2) { 925 for (d3 = 0; d3 < dim; ++d3) { 926 g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3]; 927 } 928 } 929 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w; 930 } 931 } 932 } 933 } 934 } 935 936 for (f = 0; f < NbI; ++f) { 937 for (fc = 0; fc < NcI; ++fc) { 938 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 939 const PetscInt i = offsetI+f; /* Element matrix row */ 940 for (g = 0; g < NbJ; ++g) { 941 for (gc = 0; gc < NcJ; ++gc) { 942 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 943 const PetscInt j = offsetJ+g; /* Element matrix column */ 944 const PetscInt fOff = eOffset+i*totDim+j; 945 PetscInt d, d2; 946 947 elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx]; 948 for (d = 0; d < dim; ++d) { 949 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d]; 950 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx]; 951 for (d2 = 0; d2 < dim; ++d2) { 952 elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2]; 953 } 954 } 955 } 956 } 957 } 958 } 959 } 960 if (debug > 1) { 961 PetscInt fc, f, gc, g; 962 963 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 964 for (fc = 0; fc < NcI; ++fc) { 965 for (f = 0; f < NbI; ++f) { 966 const PetscInt i = offsetI + f*NcI+fc; 967 for (gc = 0; gc < NcJ; ++gc) { 968 for (g = 0; g < NbJ; ++g) { 969 const PetscInt j = offsetJ + g*NcJ+gc; 970 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 971 } 972 } 973 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 974 } 975 } 976 } 977 cOffset += totDim; 978 cOffsetAux += totDimAux; 979 eOffset += PetscSqr(totDim); 980 } 981 PetscFunctionReturn(0); 982 } 983 984 PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 985 { 986 PetscFunctionBegin; 987 fem->ops->setfromoptions = NULL; 988 fem->ops->setup = PetscFESetUp_Basic; 989 fem->ops->view = PetscFEView_Basic; 990 fem->ops->destroy = PetscFEDestroy_Basic; 991 fem->ops->getdimension = PetscFEGetDimension_Basic; 992 fem->ops->gettabulation = PetscFEGetTabulation_Basic; 993 fem->ops->integrate = PetscFEIntegrate_Basic; 994 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 995 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 996 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 997 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 998 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 999 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 1000 PetscFunctionReturn(0); 1001 } 1002 1003 /*MC 1004 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization 1005 1006 Level: intermediate 1007 1008 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 1009 M*/ 1010 1011 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 1012 { 1013 PetscFE_Basic *b; 1014 PetscErrorCode ierr; 1015 1016 PetscFunctionBegin; 1017 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1018 ierr = PetscNewLog(fem,&b);CHKERRQ(ierr); 1019 fem->data = b; 1020 1021 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr); 1022 PetscFunctionReturn(0); 1023 } 1024 1025 1026