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 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 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 #ifndef 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 PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 278 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 279 { 280 const PetscInt debug = 0; 281 PetscPointFunc f0_func; 282 PetscPointFunc f1_func; 283 PetscQuadrature quad; 284 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 285 const PetscScalar *constants; 286 PetscReal *x; 287 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 288 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 289 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 290 PetscInt dE, Np; 291 PetscBool isAffine; 292 const PetscReal *quadPoints, *quadWeights; 293 PetscInt qNc, Nq, q; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 298 ierr = PetscFEGetQuadrature(fem, &quad);CHKERRQ(ierr); 299 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 300 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 301 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 302 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 303 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 304 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 305 ierr = PetscDSGetFieldOffset(prob, field, &fOffset);CHKERRQ(ierr); 306 ierr = PetscDSGetResidual(prob, field, &f0_func, &f1_func);CHKERRQ(ierr); 307 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 308 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 309 ierr = PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 310 ierr = PetscDSGetTabulation(prob, &B, &D);CHKERRQ(ierr); 311 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 312 if (probAux) { 313 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 314 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 315 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 316 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 317 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 318 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 319 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 320 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 321 ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 322 } 323 NbI = Nb[field]; 324 NcI = Nc[field]; 325 BI = B[field]; 326 DI = D[field]; 327 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 328 Np = cgeom->numPoints; 329 dE = cgeom->dimEmbed; 330 isAffine = cgeom->isAffine; 331 for (e = 0; e < Ne; ++e) { 332 const PetscReal *v0 = &cgeom->v[e*Np*dE]; 333 const PetscReal *J = &cgeom->J[e*Np*dE*dE]; 334 335 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 336 ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr); 337 ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 338 for (q = 0; q < Nq; ++q) { 339 const PetscReal *v; 340 const PetscReal *invJ; 341 PetscReal detJ; 342 343 if (isAffine) { 344 CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x); 345 v = x; 346 invJ = &cgeom->invJ[e*dE*dE]; 347 detJ = cgeom->detJ[e]; 348 } else { 349 v = &v0[q*dE]; 350 invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 351 detJ = cgeom->detJ[e*Np + q]; 352 } 353 if (debug > 1 && q < Np) { 354 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 355 #ifndef PETSC_USE_COMPLEX 356 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 357 #endif 358 } 359 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 360 EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 361 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 362 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]); 363 if (f1_func) { 364 ierr = PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 365 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); 366 } 367 TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL); 368 } 369 UpdateElementVec(dim, Nq, NbI, NcI, BI, DI, f0, f1, &elemVec[cOffset+fOffset]); 370 cOffset += totDim; 371 cOffsetAux += totDimAux; 372 } 373 PetscFunctionReturn(0); 374 } 375 376 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 377 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 378 { 379 const PetscInt debug = 0; 380 PetscBdPointFunc f0_func; 381 PetscBdPointFunc f1_func; 382 PetscQuadrature quad; 383 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 384 const PetscScalar *constants; 385 PetscReal *x; 386 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 387 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 388 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 389 PetscBool isAffine; 390 const PetscReal *quadPoints, *quadWeights; 391 PetscInt qNc, Nq, q, Np, dE; 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 396 ierr = PetscFEGetFaceQuadrature(fem, &quad);CHKERRQ(ierr); 397 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 398 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 399 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 400 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 401 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 402 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 403 ierr = PetscDSGetFieldOffset(prob, field, &fOffset);CHKERRQ(ierr); 404 ierr = PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);CHKERRQ(ierr); 405 if (!f0_func && !f1_func) PetscFunctionReturn(0); 406 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 407 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 408 ierr = PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 409 ierr = PetscDSGetFaceTabulation(prob, &B, &D);CHKERRQ(ierr); 410 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 411 if (probAux) { 412 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 413 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 414 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 415 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 416 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 417 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 418 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 419 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 420 ierr = PetscDSGetFaceTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 421 } 422 NbI = Nb[field]; 423 NcI = Nc[field]; 424 BI = B[field]; 425 DI = D[field]; 426 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 427 Np = fgeom->numPoints; 428 dE = fgeom->dimEmbed; 429 isAffine = fgeom->isAffine; 430 for (e = 0; e < Ne; ++e) { 431 const PetscReal *v0 = &fgeom->v[e*Np*dE]; 432 const PetscReal *J = &fgeom->J[e*Np*dE*dE]; 433 const PetscInt face = fgeom->face[e][0]; 434 435 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 436 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 437 ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr); 438 ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 439 for (q = 0; q < Nq; ++q) { 440 const PetscReal *v; 441 const PetscReal *invJ; 442 const PetscReal *n; 443 PetscReal detJ; 444 if (isAffine) { 445 CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x); 446 v = x; 447 invJ = &fgeom->suppInvJ[0][e*dE*dE]; 448 detJ = fgeom->detJ[e]; 449 n = &fgeom->n[e*dE]; 450 } else { 451 v = &v0[q*dE]; 452 invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 453 detJ = fgeom->detJ[e*Np + q]; 454 n = &fgeom->n[(e*Np+q)*dE]; 455 } 456 if (debug > 1 && q < Np) { 457 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);CHKERRQ(ierr); 458 #ifndef PETSC_USE_COMPLEX 459 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr); 460 #endif 461 } 462 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 463 EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 464 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 465 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]); 466 if (f1_func) { 467 ierr = PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr); 468 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); 469 } 470 TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL); 471 } 472 UpdateElementVec(dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], f0, f1, &elemVec[cOffset+fOffset]); 473 cOffset += totDim; 474 cOffsetAux += totDimAux; 475 } 476 PetscFunctionReturn(0); 477 } 478 479 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *geom, 480 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 481 { 482 const PetscInt debug = 0; 483 PetscPointJac g0_func; 484 PetscPointJac g1_func; 485 PetscPointJac g2_func; 486 PetscPointJac g3_func; 487 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 488 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 489 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 490 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 491 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 492 PetscQuadrature quad; 493 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 494 const PetscScalar *constants; 495 PetscReal *x; 496 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 497 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 498 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 499 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 500 PetscInt dE, Np; 501 PetscBool isAffine; 502 const PetscReal *quadPoints, *quadWeights; 503 PetscInt qNc, Nq, q; 504 PetscErrorCode ierr; 505 506 PetscFunctionBegin; 507 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 508 ierr = PetscFEGetQuadrature(fem, &quad);CHKERRQ(ierr); 509 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 510 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 511 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 512 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 513 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 514 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 515 switch(jtype) { 516 case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 517 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 518 case PETSCFE_JACOBIAN: ierr = PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 519 } 520 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0); 521 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 522 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 523 ierr = PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 524 ierr = PetscDSGetTabulation(prob, &B, &D);CHKERRQ(ierr); 525 ierr = PetscDSGetFieldOffset(prob, fieldI, &offsetI);CHKERRQ(ierr); 526 ierr = PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);CHKERRQ(ierr); 527 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 528 if (probAux) { 529 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 530 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 531 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 532 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 533 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 534 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 535 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 536 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 537 ierr = PetscDSGetTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 538 } 539 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 540 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 541 BI = B[fieldI], BJ = B[fieldJ]; 542 DI = D[fieldI], DJ = D[fieldJ]; 543 /* Initialize here in case the function is not defined */ 544 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 545 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 546 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 547 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 548 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 549 Np = geom->numPoints; 550 dE = geom->dimEmbed; 551 isAffine = geom->isAffine; 552 for (e = 0; e < Ne; ++e) { 553 const PetscReal *v0 = &geom->v[e*Np*dE]; 554 const PetscReal *J = &geom->J[e*Np*dE*dE]; 555 556 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 557 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 558 for (q = 0; q < Nq; ++q) { 559 const PetscReal *v; 560 const PetscReal *invJ; 561 PetscReal detJ; 562 const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ]; 563 const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim]; 564 PetscReal w; 565 PetscInt f, g, fc, gc, c; 566 567 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 568 if (isAffine) { 569 CoordinatesRefToReal(dE, dim, geom->xi, v0, J, &quadPoints[q*dim], x); 570 v = x; 571 invJ = &geom->invJ[e*dE*dE]; 572 detJ = geom->detJ[e]; 573 } else { 574 v = &v0[q*dE]; 575 invJ = &geom->invJ[(e*Np+q)*dE*dE]; 576 detJ = geom->detJ[e*Np + q]; 577 } 578 w = detJ*quadWeights[q]; 579 EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 580 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 581 if (g0_func) { 582 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 583 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); 584 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 585 } 586 if (g1_func) { 587 PetscInt d, d2; 588 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 589 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); 590 for (fc = 0; fc < NcI; ++fc) { 591 for (gc = 0; gc < NcJ; ++gc) { 592 for (d = 0; d < dim; ++d) { 593 g1[(fc*NcJ+gc)*dim+d] = 0.0; 594 for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 595 g1[(fc*NcJ+gc)*dim+d] *= w; 596 } 597 } 598 } 599 } 600 if (g2_func) { 601 PetscInt d, d2; 602 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 603 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); 604 for (fc = 0; fc < NcI; ++fc) { 605 for (gc = 0; gc < NcJ; ++gc) { 606 for (d = 0; d < dim; ++d) { 607 g2[(fc*NcJ+gc)*dim+d] = 0.0; 608 for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 609 g2[(fc*NcJ+gc)*dim+d] *= w; 610 } 611 } 612 } 613 } 614 if (g3_func) { 615 PetscInt d, d2, dp, d3; 616 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 617 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); 618 for (fc = 0; fc < NcI; ++fc) { 619 for (gc = 0; gc < NcJ; ++gc) { 620 for (d = 0; d < dim; ++d) { 621 for (dp = 0; dp < dim; ++dp) { 622 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0; 623 for (d2 = 0; d2 < dim; ++d2) { 624 for (d3 = 0; d3 < dim; ++d3) { 625 g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3]; 626 } 627 } 628 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w; 629 } 630 } 631 } 632 } 633 } 634 635 for (f = 0; f < NbI; ++f) { 636 for (fc = 0; fc < NcI; ++fc) { 637 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 638 const PetscInt i = offsetI+f; /* Element matrix row */ 639 for (g = 0; g < NbJ; ++g) { 640 for (gc = 0; gc < NcJ; ++gc) { 641 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 642 const PetscInt j = offsetJ+g; /* Element matrix column */ 643 const PetscInt fOff = eOffset+i*totDim+j; 644 PetscInt d, d2; 645 646 elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx]; 647 for (d = 0; d < dim; ++d) { 648 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d]; 649 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx]; 650 for (d2 = 0; d2 < dim; ++d2) { 651 elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2]; 652 } 653 } 654 } 655 } 656 } 657 } 658 } 659 if (debug > 1) { 660 PetscInt fc, f, gc, g; 661 662 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 663 for (fc = 0; fc < NcI; ++fc) { 664 for (f = 0; f < NbI; ++f) { 665 const PetscInt i = offsetI + f*NcI+fc; 666 for (gc = 0; gc < NcJ; ++gc) { 667 for (g = 0; g < NbJ; ++g) { 668 const PetscInt j = offsetJ + g*NcJ+gc; 669 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 670 } 671 } 672 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 673 } 674 } 675 } 676 cOffset += totDim; 677 cOffsetAux += totDimAux; 678 eOffset += PetscSqr(totDim); 679 } 680 PetscFunctionReturn(0); 681 } 682 683 PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 684 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 685 { 686 const PetscInt debug = 0; 687 PetscBdPointJac g0_func; 688 PetscBdPointJac g1_func; 689 PetscBdPointJac g2_func; 690 PetscBdPointJac g3_func; 691 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 692 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 693 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 694 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 695 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 696 PetscQuadrature quad; 697 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux; 698 const PetscScalar *constants; 699 PetscReal *x; 700 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 701 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 702 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 703 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 704 PetscBool isAffine; 705 const PetscReal *quadPoints, *quadWeights; 706 PetscInt qNc, Nq, q, Np, dE; 707 PetscErrorCode ierr; 708 709 PetscFunctionBegin; 710 ierr = PetscFEGetSpatialDimension(fem, &dim);CHKERRQ(ierr); 711 ierr = PetscFEGetFaceQuadrature(fem, &quad);CHKERRQ(ierr); 712 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 713 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 714 ierr = PetscDSGetDimensions(prob, &Nb);CHKERRQ(ierr); 715 ierr = PetscDSGetComponents(prob, &Nc);CHKERRQ(ierr); 716 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 717 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 718 ierr = PetscDSGetFieldOffset(prob, fieldI, &offsetI);CHKERRQ(ierr); 719 ierr = PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);CHKERRQ(ierr); 720 ierr = PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr); 721 ierr = PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 722 ierr = PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);CHKERRQ(ierr); 723 ierr = PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 724 ierr = PetscDSGetFaceTabulation(prob, &B, &D);CHKERRQ(ierr); 725 ierr = PetscDSGetConstants(prob, &numConstants, &constants);CHKERRQ(ierr); 726 if (probAux) { 727 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 728 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 729 ierr = PetscDSGetDimensions(probAux, &NbAux);CHKERRQ(ierr); 730 ierr = PetscDSGetComponents(probAux, &NcAux);CHKERRQ(ierr); 731 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 732 ierr = PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);CHKERRQ(ierr); 733 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 734 ierr = PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);CHKERRQ(ierr); 735 ierr = PetscDSGetFaceTabulation(probAux, &BAux, &DAux);CHKERRQ(ierr); 736 } 737 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 738 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 739 BI = B[fieldI], BJ = B[fieldJ]; 740 DI = D[fieldI], DJ = D[fieldJ]; 741 /* Initialize here in case the function is not defined */ 742 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 743 ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 744 ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 745 ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 746 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 747 Np = fgeom->numPoints; 748 dE = fgeom->dimEmbed; 749 isAffine = fgeom->isAffine; 750 for (e = 0; e < Ne; ++e) { 751 const PetscReal *v0 = &fgeom->v[e*Np*dE]; 752 const PetscReal *J = &fgeom->J[e*Np*dE*dE]; 753 const PetscInt face = fgeom->face[e][0]; 754 755 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 756 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 757 for (q = 0; q < Nq; ++q) { 758 const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ]; 759 const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim]; 760 PetscReal w; 761 PetscInt f, g, fc, gc, c; 762 const PetscReal *v; 763 const PetscReal *invJ; 764 const PetscReal *n; 765 PetscReal detJ; 766 767 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 768 if (isAffine) { 769 CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x); 770 v = x; 771 invJ = &fgeom->suppInvJ[0][e*dE*dE]; 772 detJ = fgeom->detJ[e]; 773 n = &fgeom->n[e*dE]; 774 } else { 775 v = &v0[q*dE]; 776 invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 777 detJ = fgeom->detJ[e*Np + q]; 778 n = &fgeom->n[(e*Np+q)*dE]; 779 } 780 w = detJ*quadWeights[q]; 781 782 EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t); 783 if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL); 784 if (g0_func) { 785 ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr); 786 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); 787 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 788 } 789 if (g1_func) { 790 PetscInt d, d2; 791 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 792 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); 793 for (fc = 0; fc < NcI; ++fc) { 794 for (gc = 0; gc < NcJ; ++gc) { 795 for (d = 0; d < dim; ++d) { 796 g1[(fc*NcJ+gc)*dim+d] = 0.0; 797 for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 798 g1[(fc*NcJ+gc)*dim+d] *= w; 799 } 800 } 801 } 802 } 803 if (g2_func) { 804 PetscInt d, d2; 805 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr); 806 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); 807 for (fc = 0; fc < NcI; ++fc) { 808 for (gc = 0; gc < NcJ; ++gc) { 809 for (d = 0; d < dim; ++d) { 810 g2[(fc*NcJ+gc)*dim+d] = 0.0; 811 for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2]; 812 g2[(fc*NcJ+gc)*dim+d] *= w; 813 } 814 } 815 } 816 } 817 if (g3_func) { 818 PetscInt d, d2, dp, d3; 819 ierr = PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr); 820 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); 821 for (fc = 0; fc < NcI; ++fc) { 822 for (gc = 0; gc < NcJ; ++gc) { 823 for (d = 0; d < dim; ++d) { 824 for (dp = 0; dp < dim; ++dp) { 825 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0; 826 for (d2 = 0; d2 < dim; ++d2) { 827 for (d3 = 0; d3 < dim; ++d3) { 828 g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3]; 829 } 830 } 831 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w; 832 } 833 } 834 } 835 } 836 } 837 838 for (f = 0; f < NbI; ++f) { 839 for (fc = 0; fc < NcI; ++fc) { 840 const PetscInt fidx = f*NcI+fc; /* Test function basis index */ 841 const PetscInt i = offsetI+f; /* Element matrix row */ 842 for (g = 0; g < NbJ; ++g) { 843 for (gc = 0; gc < NcJ; ++gc) { 844 const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */ 845 const PetscInt j = offsetJ+g; /* Element matrix column */ 846 const PetscInt fOff = eOffset+i*totDim+j; 847 PetscInt d, d2; 848 849 elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx]; 850 for (d = 0; d < dim; ++d) { 851 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d]; 852 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx]; 853 for (d2 = 0; d2 < dim; ++d2) { 854 elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2]; 855 } 856 } 857 } 858 } 859 } 860 } 861 } 862 if (debug > 1) { 863 PetscInt fc, f, gc, g; 864 865 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 866 for (fc = 0; fc < NcI; ++fc) { 867 for (f = 0; f < NbI; ++f) { 868 const PetscInt i = offsetI + f*NcI+fc; 869 for (gc = 0; gc < NcJ; ++gc) { 870 for (g = 0; g < NbJ; ++g) { 871 const PetscInt j = offsetJ + g*NcJ+gc; 872 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 873 } 874 } 875 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 876 } 877 } 878 } 879 cOffset += totDim; 880 cOffsetAux += totDimAux; 881 eOffset += PetscSqr(totDim); 882 } 883 PetscFunctionReturn(0); 884 } 885 886 PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 887 { 888 PetscFunctionBegin; 889 fem->ops->setfromoptions = NULL; 890 fem->ops->setup = PetscFESetUp_Basic; 891 fem->ops->view = PetscFEView_Basic; 892 fem->ops->destroy = PetscFEDestroy_Basic; 893 fem->ops->getdimension = PetscFEGetDimension_Basic; 894 fem->ops->gettabulation = PetscFEGetTabulation_Basic; 895 fem->ops->integrate = PetscFEIntegrate_Basic; 896 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 897 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 898 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 899 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 900 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 901 PetscFunctionReturn(0); 902 } 903 904 /*MC 905 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization 906 907 Level: intermediate 908 909 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 910 M*/ 911 912 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 913 { 914 PetscFE_Basic *b; 915 PetscErrorCode ierr; 916 917 PetscFunctionBegin; 918 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 919 ierr = PetscNewLog(fem,&b);CHKERRQ(ierr); 920 fem->data = b; 921 922 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr); 923 PetscFunctionReturn(0); 924 } 925 926 927