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