1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscblaslapack.h> 3 4 static PetscErrorCode PetscFEDestroy_Basic(PetscFE fem) 5 { 6 PetscFE_Basic *b = (PetscFE_Basic *)fem->data; 7 8 PetscFunctionBegin; 9 PetscCall(PetscFree(b)); 10 PetscFunctionReturn(PETSC_SUCCESS); 11 } 12 13 static PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer v) 14 { 15 PetscInt dim, Nc; 16 PetscSpace basis = NULL; 17 PetscDualSpace dual = NULL; 18 PetscQuadrature quad = NULL; 19 20 PetscFunctionBegin; 21 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 22 PetscCall(PetscFEGetNumComponents(fe, &Nc)); 23 PetscCall(PetscFEGetBasisSpace(fe, &basis)); 24 PetscCall(PetscFEGetDualSpace(fe, &dual)); 25 PetscCall(PetscFEGetQuadrature(fe, &quad)); 26 PetscCall(PetscViewerASCIIPushTab(v)); 27 PetscCall(PetscViewerASCIIPrintf(v, "Basic Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components\n", dim, Nc)); 28 if (basis) PetscCall(PetscSpaceView(basis, v)); 29 if (dual) PetscCall(PetscDualSpaceView(dual, v)); 30 if (quad) PetscCall(PetscQuadratureView(quad, v)); 31 PetscCall(PetscViewerASCIIPopTab(v)); 32 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 static PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v) 36 { 37 PetscBool iascii; 38 39 PetscFunctionBegin; 40 PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii)); 41 if (iascii) PetscCall(PetscFEView_Basic_Ascii(fe, v)); 42 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 /* Construct the change of basis from prime basis to nodal basis */ 46 PETSC_INTERN PetscErrorCode PetscFESetUp_Basic(PetscFE fem) 47 { 48 PetscReal *work; 49 PetscBLASInt *pivots; 50 PetscBLASInt n, info; 51 PetscInt pdim, j; 52 53 PetscFunctionBegin; 54 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim)); 55 PetscCall(PetscMalloc1(pdim * pdim, &fem->invV)); 56 for (j = 0; j < pdim; ++j) { 57 PetscReal *Bf; 58 PetscQuadrature f; 59 const PetscReal *points, *weights; 60 PetscInt Nc, Nq, q, k, c; 61 62 PetscCall(PetscDualSpaceGetFunctional(fem->dualSpace, j, &f)); 63 PetscCall(PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights)); 64 PetscCall(PetscMalloc1(Nc * Nq * pdim, &Bf)); 65 PetscCall(PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL)); 66 for (k = 0; k < pdim; ++k) { 67 /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */ 68 fem->invV[j * pdim + k] = 0.0; 69 70 for (q = 0; q < Nq; ++q) { 71 for (c = 0; c < Nc; ++c) fem->invV[j * pdim + k] += Bf[(q * pdim + k) * Nc + c] * weights[q * Nc + c]; 72 } 73 } 74 PetscCall(PetscFree(Bf)); 75 } 76 77 PetscCall(PetscMalloc2(pdim, &pivots, pdim, &work)); 78 n = pdim; 79 PetscCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info)); 80 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info); 81 PetscCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info)); 82 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info); 83 PetscCall(PetscFree2(pivots, work)); 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86 87 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim) 88 { 89 PetscFunctionBegin; 90 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, dim)); 91 PetscFunctionReturn(PETSC_SUCCESS); 92 } 93 94 /* Tensor contraction on the middle index, 95 * C[m,n,p] = A[m,k,p] * B[k,n] 96 * where all matrices use C-style ordering. 97 */ 98 static PetscErrorCode TensorContract_Private(PetscInt m, PetscInt n, PetscInt p, PetscInt k, const PetscReal *A, const PetscReal *B, PetscReal *C) 99 { 100 PetscInt i; 101 102 PetscFunctionBegin; 103 PetscCheck(n && p, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Empty tensor is not allowed %" PetscInt_FMT " %" PetscInt_FMT, n, p); 104 for (i = 0; i < m; i++) { 105 PetscBLASInt n_, p_, k_, lda, ldb, ldc; 106 PetscReal one = 1, zero = 0; 107 /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n] 108 * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k) 109 */ 110 PetscCall(PetscBLASIntCast(n, &n_)); 111 PetscCall(PetscBLASIntCast(p, &p_)); 112 PetscCall(PetscBLASIntCast(k, &k_)); 113 lda = p_; 114 ldb = n_; 115 ldc = p_; 116 PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc)); 117 } 118 PetscCall(PetscLogFlops(2. * m * n * p * k)); 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 PETSC_INTERN PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 123 { 124 DM dm; 125 PetscInt pdim; /* Dimension of FE space P */ 126 PetscInt dim; /* Spatial dimension */ 127 PetscInt Nc; /* Field components */ 128 PetscReal *B = K >= 0 ? T->T[0] : NULL; 129 PetscReal *D = K >= 1 ? T->T[1] : NULL; 130 PetscReal *H = K >= 2 ? T->T[2] : NULL; 131 PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL; 132 133 PetscFunctionBegin; 134 PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm)); 135 PetscCall(DMGetDimension(dm, &dim)); 136 PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim)); 137 PetscCall(PetscFEGetNumComponents(fem, &Nc)); 138 /* Evaluate the prime basis functions at all points */ 139 if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB)); 140 if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD)); 141 if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH)); 142 PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH)); 143 /* Translate from prime to nodal basis */ 144 if (B) { 145 /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */ 146 PetscCall(TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B)); 147 } 148 if (D && dim) { 149 /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */ 150 PetscCall(TensorContract_Private(npoints, pdim, Nc * dim, pdim, tmpD, fem->invV, D)); 151 } 152 if (H && dim) { 153 /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */ 154 PetscCall(TensorContract_Private(npoints, pdim, Nc * dim * dim, pdim, tmpH, fem->invV, H)); 155 } 156 if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB)); 157 if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD)); 158 if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH)); 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 162 static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 163 { 164 const PetscInt debug = 0; 165 PetscFE fe; 166 PetscPointFunc obj_func; 167 PetscQuadrature quad; 168 PetscTabulation *T, *TAux = NULL; 169 PetscScalar *u, *u_x, *a, *a_x; 170 const PetscScalar *constants; 171 PetscReal *x; 172 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 173 PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 174 PetscBool isAffine; 175 const PetscReal *quadPoints, *quadWeights; 176 PetscInt qNc, Nq, q; 177 178 PetscFunctionBegin; 179 PetscCall(PetscDSGetObjective(ds, field, &obj_func)); 180 if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS); 181 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 182 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 183 PetscCall(PetscFEGetQuadrature(fe, &quad)); 184 PetscCall(PetscDSGetNumFields(ds, &Nf)); 185 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 186 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 187 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 188 PetscCall(PetscDSGetTabulation(ds, &T)); 189 PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x)); 190 PetscCall(PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL)); 191 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 192 if (dsAux) { 193 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 194 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 195 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 196 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 197 PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 198 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 199 PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np); 200 } 201 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 202 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 203 Np = cgeom->numPoints; 204 dE = cgeom->dimEmbed; 205 isAffine = cgeom->isAffine; 206 for (e = 0; e < Ne; ++e) { 207 PetscFEGeom fegeom; 208 209 fegeom.dim = cgeom->dim; 210 fegeom.dimEmbed = cgeom->dimEmbed; 211 if (isAffine) { 212 fegeom.v = x; 213 fegeom.xi = cgeom->xi; 214 fegeom.J = &cgeom->J[e * Np * dE * dE]; 215 fegeom.invJ = &cgeom->invJ[e * Np * dE * dE]; 216 fegeom.detJ = &cgeom->detJ[e * Np]; 217 } 218 for (q = 0; q < Nq; ++q) { 219 PetscScalar integrand; 220 PetscReal w; 221 222 if (isAffine) { 223 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x); 224 } else { 225 fegeom.v = &cgeom->v[(e * Np + q) * dE]; 226 fegeom.J = &cgeom->J[(e * Np + q) * dE * dE]; 227 fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE]; 228 fegeom.detJ = &cgeom->detJ[e * Np + q]; 229 } 230 w = fegeom.detJ[0] * quadWeights[q]; 231 if (debug > 1 && q < Np) { 232 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 233 #if !defined(PETSC_USE_COMPLEX) 234 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 235 #endif 236 } 237 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 238 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL)); 239 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 240 obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, fegeom.v, numConstants, constants, &integrand); 241 integrand *= w; 242 integral[e * Nf + field] += integrand; 243 if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[field]))); 244 } 245 cOffset += totDim; 246 cOffsetAux += totDimAux; 247 } 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 251 static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, PetscBdPointFunc obj_func, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 252 { 253 const PetscInt debug = 0; 254 PetscFE fe; 255 PetscQuadrature quad; 256 PetscTabulation *Tf, *TfAux = NULL; 257 PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal; 258 const PetscScalar *constants; 259 PetscReal *x; 260 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 261 PetscBool isAffine, auxOnBd; 262 const PetscReal *quadPoints, *quadWeights; 263 PetscInt qNc, Nq, q, Np, dE; 264 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 265 266 PetscFunctionBegin; 267 if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS); 268 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 269 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 270 PetscCall(PetscFEGetFaceQuadrature(fe, &quad)); 271 PetscCall(PetscDSGetNumFields(ds, &Nf)); 272 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 273 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 274 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 275 PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x)); 276 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 277 PetscCall(PetscDSGetFaceTabulation(ds, &Tf)); 278 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 279 if (dsAux) { 280 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 281 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 282 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 283 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 284 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 285 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 286 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 287 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux)); 288 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux)); 289 PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np); 290 } 291 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 292 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 293 Np = fgeom->numPoints; 294 dE = fgeom->dimEmbed; 295 isAffine = fgeom->isAffine; 296 for (e = 0; e < Ne; ++e) { 297 PetscFEGeom fegeom, cgeom; 298 const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */ 299 fegeom.n = NULL; 300 fegeom.v = NULL; 301 fegeom.J = NULL; 302 fegeom.detJ = NULL; 303 fegeom.dim = fgeom->dim; 304 fegeom.dimEmbed = fgeom->dimEmbed; 305 cgeom.dim = fgeom->dim; 306 cgeom.dimEmbed = fgeom->dimEmbed; 307 if (isAffine) { 308 fegeom.v = x; 309 fegeom.xi = fgeom->xi; 310 fegeom.J = &fgeom->J[e * Np * dE * dE]; 311 fegeom.invJ = &fgeom->invJ[e * Np * dE * dE]; 312 fegeom.detJ = &fgeom->detJ[e * Np]; 313 fegeom.n = &fgeom->n[e * Np * dE]; 314 315 cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE]; 316 cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE]; 317 cgeom.detJ = &fgeom->suppDetJ[0][e * Np]; 318 } 319 for (q = 0; q < Nq; ++q) { 320 PetscScalar integrand; 321 PetscReal w; 322 323 if (isAffine) { 324 CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x); 325 } else { 326 fegeom.v = &fgeom->v[(e * Np + q) * dE]; 327 fegeom.J = &fgeom->J[(e * Np + q) * dE * dE]; 328 fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE]; 329 fegeom.detJ = &fgeom->detJ[e * Np + q]; 330 fegeom.n = &fgeom->n[(e * Np + q) * dE]; 331 332 cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE]; 333 cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE]; 334 cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q]; 335 } 336 w = fegeom.detJ[0] * quadWeights[q]; 337 if (debug > 1 && q < Np) { 338 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 339 #ifndef PETSC_USE_COMPLEX 340 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 341 #endif 342 } 343 if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 344 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL)); 345 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 346 obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, fegeom.v, fegeom.n, numConstants, constants, &integrand); 347 integrand *= w; 348 integral[e * Nf + field] += integrand; 349 if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field]))); 350 } 351 cOffset += totDim; 352 cOffsetAux += totDimAux; 353 } 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 358 { 359 const PetscInt debug = 0; 360 const PetscInt field = key.field; 361 PetscFE fe; 362 PetscWeakForm wf; 363 PetscInt n0, n1, i; 364 PetscPointFunc *f0_func, *f1_func; 365 PetscQuadrature quad; 366 PetscTabulation *T, *TAux = NULL; 367 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 368 const PetscScalar *constants; 369 PetscReal *x; 370 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 371 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e; 372 const PetscReal *quadPoints, *quadWeights; 373 PetscInt qdim, qNc, Nq, q, dE; 374 375 PetscFunctionBegin; 376 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 377 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 378 PetscCall(PetscFEGetQuadrature(fe, &quad)); 379 PetscCall(PetscDSGetNumFields(ds, &Nf)); 380 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 381 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 382 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 383 PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset)); 384 PetscCall(PetscDSGetWeakForm(ds, &wf)); 385 PetscCall(PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func)); 386 if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS); 387 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 388 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 389 PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL)); 390 PetscCall(PetscDSGetTabulation(ds, &T)); 391 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 392 if (dsAux) { 393 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 394 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 395 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 396 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 397 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 398 PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 399 PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np); 400 } 401 PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights)); 402 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 403 dE = cgeom->dimEmbed; 404 PetscCheck(cgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", cgeom->dim, qdim); 405 for (e = 0; e < Ne; ++e) { 406 PetscFEGeom fegeom; 407 408 fegeom.v = x; /* workspace */ 409 PetscCall(PetscArrayzero(f0, Nq * T[field]->Nc)); 410 PetscCall(PetscArrayzero(f1, Nq * T[field]->Nc * dE)); 411 for (q = 0; q < Nq; ++q) { 412 PetscReal w; 413 PetscInt c, d; 414 415 PetscCall(PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom)); 416 w = fegeom.detJ[0] * quadWeights[q]; 417 if (debug > 1 && q < cgeom->numPoints) { 418 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 419 #if !defined(PETSC_USE_COMPLEX) 420 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 421 #endif 422 } 423 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 424 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 425 for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q * T[field]->Nc]); 426 for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w; 427 for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q * T[field]->Nc * dim]); 428 for (c = 0; c < T[field]->Nc; ++c) 429 for (d = 0; d < dim; ++d) f1[(q * T[field]->Nc + c) * dim + d] *= w; 430 if (debug) { 431 PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " wt %g\n", q, (double)quadWeights[q])); 432 if (debug > 2) { 433 PetscCall(PetscPrintf(PETSC_COMM_SELF, " field %" PetscInt_FMT ":", field)); 434 for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c]))); 435 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 436 PetscCall(PetscPrintf(PETSC_COMM_SELF, " resid %" PetscInt_FMT ":", field)); 437 for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c]))); 438 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 439 } 440 } 441 } 442 PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset])); 443 cOffset += totDim; 444 cOffsetAux += totDimAux; 445 } 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 450 { 451 const PetscInt debug = 0; 452 const PetscInt field = key.field; 453 PetscFE fe; 454 PetscInt n0, n1, i; 455 PetscBdPointFunc *f0_func, *f1_func; 456 PetscQuadrature quad; 457 PetscTabulation *Tf, *TfAux = NULL; 458 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 459 const PetscScalar *constants; 460 PetscReal *x; 461 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 462 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI; 463 PetscBool auxOnBd = PETSC_FALSE; 464 const PetscReal *quadPoints, *quadWeights; 465 PetscInt qdim, qNc, Nq, q, dE; 466 467 PetscFunctionBegin; 468 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 469 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 470 PetscCall(PetscFEGetFaceQuadrature(fe, &quad)); 471 PetscCall(PetscDSGetNumFields(ds, &Nf)); 472 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 473 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 474 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 475 PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset)); 476 PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func)); 477 if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS); 478 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 479 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 480 PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL)); 481 PetscCall(PetscDSGetFaceTabulation(ds, &Tf)); 482 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 483 if (dsAux) { 484 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 485 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 486 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 487 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 488 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 489 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 490 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 491 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux)); 492 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux)); 493 PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np); 494 } 495 NcI = Tf[field]->Nc; 496 PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights)); 497 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 498 dE = fgeom->dimEmbed; 499 /* TODO FIX THIS */ 500 fgeom->dim = dim - 1; 501 PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim); 502 for (e = 0; e < Ne; ++e) { 503 PetscFEGeom fegeom, cgeom; 504 const PetscInt face = fgeom->face[e][0]; 505 506 fegeom.v = x; /* Workspace */ 507 PetscCall(PetscArrayzero(f0, Nq * NcI)); 508 PetscCall(PetscArrayzero(f1, Nq * NcI * dE)); 509 for (q = 0; q < Nq; ++q) { 510 PetscReal w; 511 PetscInt c, d; 512 513 PetscCall(PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom)); 514 PetscCall(PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom)); 515 w = fegeom.detJ[0] * quadWeights[q]; 516 if (debug > 1) { 517 if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) { 518 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 519 #if !defined(PETSC_USE_COMPLEX) 520 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 521 PetscCall(DMPrintCellVector(e, "n", dim, fegeom.n)); 522 #endif 523 } 524 } 525 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 526 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 527 for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q * NcI]); 528 for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w; 529 for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q * NcI * dim]); 530 for (c = 0; c < NcI; ++c) 531 for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w; 532 if (debug) { 533 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q)); 534 for (c = 0; c < NcI; ++c) { 535 if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c]))); 536 if (n1) { 537 for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcI + c) * dim + d]))); 538 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 539 } 540 } 541 } 542 } 543 PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 544 cOffset += totDim; 545 cOffsetAux += totDimAux; 546 } 547 PetscFunctionReturn(PETSC_SUCCESS); 548 } 549 550 /* 551 BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but 552 all transforms operate in the full space and are square. 553 554 HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square. 555 1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces 556 2) We need to assume that the orientation is 0 for both 557 3) TODO We need to use a non-square Jacobian for the derivative maps, meaning the embedding dimension has to go to EvaluateFieldJets() and UpdateElementVec() 558 */ 559 static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 560 { 561 const PetscInt debug = 0; 562 const PetscInt field = key.field; 563 PetscFE fe; 564 PetscWeakForm wf; 565 PetscInt n0, n1, i; 566 PetscBdPointFunc *f0_func, *f1_func; 567 PetscQuadrature quad; 568 DMPolytopeType ct; 569 PetscTabulation *Tf, *TfIn, *TfAux = NULL; 570 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 571 const PetscScalar *constants; 572 PetscReal *x; 573 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 574 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS; 575 PetscBool isCohesiveField, auxOnBd = PETSC_FALSE; 576 const PetscReal *quadPoints, *quadWeights; 577 PetscInt qdim, qNc, Nq, q, dE; 578 579 PetscFunctionBegin; 580 /* Hybrid discretization is posed directly on faces */ 581 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 582 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 583 PetscCall(PetscFEGetQuadrature(fe, &quad)); 584 PetscCall(PetscDSGetNumFields(ds, &Nf)); 585 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 586 PetscCall(PetscDSGetTotalDimension(dsIn, &totDimIn)); 587 PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, s, &uOff)); 588 PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(dsIn, s, &uOff_x)); 589 PetscCall(PetscDSGetFieldOffsetCohesive(ds, field, &fOffset)); 590 PetscCall(PetscDSGetWeakForm(ds, &wf)); 591 PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func)); 592 if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS); 593 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 594 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 595 PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL)); 596 /* NOTE This is a bulk tabulation because the DS is a face discretization */ 597 PetscCall(PetscDSGetTabulation(ds, &Tf)); 598 PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn)); 599 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 600 if (dsAux) { 601 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 602 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 603 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 604 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 605 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 606 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 607 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE; 608 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux)); 609 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux)); 610 PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np); 611 } 612 PetscCall(PetscDSGetCohesive(ds, field, &isCohesiveField)); 613 NcI = Tf[field]->Nc; 614 NcS = NcI; 615 PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights)); 616 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 617 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 618 dE = fgeom->dimEmbed; 619 PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim); 620 for (e = 0; e < Ne; ++e) { 621 PetscFEGeom fegeom; 622 const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]}; 623 const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]}; 624 625 fegeom.v = x; /* Workspace */ 626 PetscCall(PetscArrayzero(f0, Nq * NcS)); 627 PetscCall(PetscArrayzero(f1, Nq * NcS * dE)); 628 for (q = 0; q < Nq; ++q) { 629 PetscInt qpt[2]; 630 PetscReal w; 631 PetscInt c, d; 632 633 PetscCall(PetscDSPermuteQuadPoint(ds, ornt[0], field, q, &qpt[0])); 634 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], 0), field, q, &qpt[1])); 635 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 636 w = fegeom.detJ[0] * quadWeights[q]; 637 if (debug > 1 && q < fgeom->numPoints) { 638 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 639 #if !defined(PETSC_USE_COMPLEX) 640 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ)); 641 #endif 642 } 643 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 644 /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */ 645 PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, &coefficients[cOffsetIn], &coefficients_t[cOffsetIn], u, u_x, u_t)); 646 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 647 for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q * NcS]); 648 for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w; 649 for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q * NcS * dE]); 650 for (c = 0; c < NcS; ++c) 651 for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w; 652 } 653 if (isCohesiveField) { 654 PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 655 } else { 656 PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 657 } 658 cOffset += totDim; 659 cOffsetIn += totDimIn; 660 cOffsetAux += totDimAux; 661 } 662 PetscFunctionReturn(PETSC_SUCCESS); 663 } 664 665 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 666 { 667 const PetscInt debug = 0; 668 PetscFE feI, feJ; 669 PetscWeakForm wf; 670 PetscPointJac *g0_func, *g1_func, *g2_func, *g3_func; 671 PetscInt n0, n1, n2, n3, i; 672 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 673 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 674 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 675 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 676 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 677 PetscQuadrature quad; 678 PetscTabulation *T, *TAux = NULL; 679 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 680 const PetscScalar *constants; 681 PetscReal *x; 682 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 683 PetscInt NcI = 0, NcJ = 0; 684 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 685 PetscInt dE, Np; 686 PetscBool isAffine; 687 const PetscReal *quadPoints, *quadWeights; 688 PetscInt qNc, Nq, q; 689 690 PetscFunctionBegin; 691 PetscCall(PetscDSGetNumFields(ds, &Nf)); 692 fieldI = key.field / Nf; 693 fieldJ = key.field % Nf; 694 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 695 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 696 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 697 PetscCall(PetscFEGetQuadrature(feI, &quad)); 698 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 699 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 700 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 701 PetscCall(PetscDSGetWeakForm(ds, &wf)); 702 switch (jtype) { 703 case PETSCFE_JACOBIAN_DYN: 704 PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 705 break; 706 case PETSCFE_JACOBIAN_PRE: 707 PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 708 break; 709 case PETSCFE_JACOBIAN: 710 PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 711 break; 712 } 713 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 714 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 715 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 716 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 717 PetscCall(PetscDSGetTabulation(ds, &T)); 718 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 719 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 720 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 721 if (dsAux) { 722 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 723 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 724 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 725 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 726 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 727 PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 728 PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np); 729 } 730 NcI = T[fieldI]->Nc; 731 NcJ = T[fieldJ]->Nc; 732 Np = cgeom->numPoints; 733 dE = cgeom->dimEmbed; 734 isAffine = cgeom->isAffine; 735 /* Initialize here in case the function is not defined */ 736 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 737 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 738 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 739 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 740 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 741 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 742 for (e = 0; e < Ne; ++e) { 743 PetscFEGeom fegeom; 744 745 fegeom.dim = cgeom->dim; 746 fegeom.dimEmbed = cgeom->dimEmbed; 747 if (isAffine) { 748 fegeom.v = x; 749 fegeom.xi = cgeom->xi; 750 fegeom.J = &cgeom->J[e * Np * dE * dE]; 751 fegeom.invJ = &cgeom->invJ[e * Np * dE * dE]; 752 fegeom.detJ = &cgeom->detJ[e * Np]; 753 } 754 for (q = 0; q < Nq; ++q) { 755 PetscReal w; 756 PetscInt c; 757 758 if (isAffine) { 759 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x); 760 } else { 761 fegeom.v = &cgeom->v[(e * Np + q) * dE]; 762 fegeom.J = &cgeom->J[(e * Np + q) * dE * dE]; 763 fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE]; 764 fegeom.detJ = &cgeom->detJ[e * Np + q]; 765 } 766 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 767 w = fegeom.detJ[0] * quadWeights[q]; 768 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 769 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 770 if (n0) { 771 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 772 for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g0); 773 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 774 } 775 if (n1) { 776 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 777 for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g1); 778 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 779 } 780 if (n2) { 781 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 782 for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g2); 783 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 784 } 785 if (n3) { 786 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 787 for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g3); 788 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 789 } 790 791 PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 792 } 793 if (debug > 1) { 794 PetscInt fc, f, gc, g; 795 796 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 797 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 798 for (f = 0; f < T[fieldI]->Nb; ++f) { 799 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 800 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 801 for (g = 0; g < T[fieldJ]->Nb; ++g) { 802 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 803 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]))); 804 } 805 } 806 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 807 } 808 } 809 } 810 cOffset += totDim; 811 cOffsetAux += totDimAux; 812 eOffset += PetscSqr(totDim); 813 } 814 PetscFunctionReturn(PETSC_SUCCESS); 815 } 816 817 static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 818 { 819 const PetscInt debug = 0; 820 PetscFE feI, feJ; 821 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 822 PetscInt n0, n1, n2, n3, i; 823 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 824 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 825 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 826 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 827 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 828 PetscQuadrature quad; 829 PetscTabulation *T, *TAux = NULL; 830 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 831 const PetscScalar *constants; 832 PetscReal *x; 833 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 834 PetscInt NcI = 0, NcJ = 0; 835 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 836 PetscBool isAffine; 837 const PetscReal *quadPoints, *quadWeights; 838 PetscInt qNc, Nq, q, Np, dE; 839 840 PetscFunctionBegin; 841 PetscCall(PetscDSGetNumFields(ds, &Nf)); 842 fieldI = key.field / Nf; 843 fieldJ = key.field % Nf; 844 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 845 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 846 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 847 PetscCall(PetscFEGetFaceQuadrature(feI, &quad)); 848 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 849 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 850 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 851 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 852 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 853 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 854 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 855 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 856 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 857 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 858 PetscCall(PetscDSGetFaceTabulation(ds, &T)); 859 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 860 if (dsAux) { 861 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 862 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 863 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 864 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 865 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 866 PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 867 } 868 NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc; 869 Np = fgeom->numPoints; 870 dE = fgeom->dimEmbed; 871 isAffine = fgeom->isAffine; 872 /* Initialize here in case the function is not defined */ 873 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 874 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 875 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 876 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 877 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 878 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 879 for (e = 0; e < Ne; ++e) { 880 PetscFEGeom fegeom, cgeom; 881 const PetscInt face = fgeom->face[e][0]; 882 fegeom.n = NULL; 883 fegeom.v = NULL; 884 fegeom.J = NULL; 885 fegeom.detJ = NULL; 886 fegeom.dim = fgeom->dim; 887 fegeom.dimEmbed = fgeom->dimEmbed; 888 cgeom.dim = fgeom->dim; 889 cgeom.dimEmbed = fgeom->dimEmbed; 890 if (isAffine) { 891 fegeom.v = x; 892 fegeom.xi = fgeom->xi; 893 fegeom.J = &fgeom->J[e * Np * dE * dE]; 894 fegeom.invJ = &fgeom->invJ[e * Np * dE * dE]; 895 fegeom.detJ = &fgeom->detJ[e * Np]; 896 fegeom.n = &fgeom->n[e * Np * dE]; 897 898 cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE]; 899 cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE]; 900 cgeom.detJ = &fgeom->suppDetJ[0][e * Np]; 901 } 902 for (q = 0; q < Nq; ++q) { 903 PetscReal w; 904 PetscInt c; 905 906 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 907 if (isAffine) { 908 CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x); 909 } else { 910 fegeom.v = &fgeom->v[(e * Np + q) * dE]; 911 fegeom.J = &fgeom->J[(e * Np + q) * dE * dE]; 912 fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE]; 913 fegeom.detJ = &fgeom->detJ[e * Np + q]; 914 fegeom.n = &fgeom->n[(e * Np + q) * dE]; 915 916 cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE]; 917 cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE]; 918 cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q]; 919 } 920 w = fegeom.detJ[0] * quadWeights[q]; 921 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 922 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 923 if (n0) { 924 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 925 for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0); 926 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 927 } 928 if (n1) { 929 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 930 for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1); 931 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 932 } 933 if (n2) { 934 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 935 for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2); 936 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 937 } 938 if (n3) { 939 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 940 for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3); 941 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 942 } 943 944 PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 945 } 946 if (debug > 1) { 947 PetscInt fc, f, gc, g; 948 949 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 950 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 951 for (f = 0; f < T[fieldI]->Nb; ++f) { 952 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 953 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 954 for (g = 0; g < T[fieldJ]->Nb; ++g) { 955 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 956 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]))); 957 } 958 } 959 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 960 } 961 } 962 } 963 cOffset += totDim; 964 cOffsetAux += totDimAux; 965 eOffset += PetscSqr(totDim); 966 } 967 PetscFunctionReturn(PETSC_SUCCESS); 968 } 969 970 PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscDS dsIn, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 971 { 972 const PetscInt debug = 0; 973 PetscFE feI, feJ; 974 PetscWeakForm wf; 975 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 976 PetscInt n0, n1, n2, n3, i; 977 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 978 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 979 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 980 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 981 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 982 PetscQuadrature quad; 983 DMPolytopeType ct; 984 PetscTabulation *T, *TfIn, *TAux = NULL; 985 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 986 const PetscScalar *constants; 987 PetscReal *x; 988 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 989 PetscInt NcI = 0, NcJ = 0, NcS, NcT; 990 PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 991 PetscBool isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE; 992 const PetscReal *quadPoints, *quadWeights; 993 PetscInt qNc, Nq, q, dE; 994 995 PetscFunctionBegin; 996 PetscCall(PetscDSGetNumFields(ds, &Nf)); 997 fieldI = key.field / Nf; 998 fieldJ = key.field % Nf; 999 /* Hybrid discretization is posed directly on faces */ 1000 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 1001 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 1002 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 1003 PetscCall(PetscFEGetQuadrature(feI, &quad)); 1004 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 1005 PetscCall(PetscDSGetComponentOffsetsCohesive(ds, s, &uOff)); 1006 PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x)); 1007 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1008 switch (jtype) { 1009 case PETSCFE_JACOBIAN_PRE: 1010 PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1011 break; 1012 case PETSCFE_JACOBIAN: 1013 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1014 break; 1015 case PETSCFE_JACOBIAN_DYN: 1016 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)"); 1017 } 1018 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 1019 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 1020 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 1021 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 1022 PetscCall(PetscDSGetTabulation(ds, &T)); 1023 PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn)); 1024 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI)); 1025 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ)); 1026 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 1027 if (dsAux) { 1028 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 1029 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 1030 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 1031 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 1032 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 1033 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 1034 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE; 1035 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 1036 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 1037 PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np); 1038 } 1039 PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI)); 1040 PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ)); 1041 NcI = T[fieldI]->Nc; 1042 NcJ = T[fieldJ]->Nc; 1043 NcS = isCohesiveFieldI ? NcI : 2 * NcI; 1044 NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ; 1045 dE = fgeom->dimEmbed; 1046 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1047 PetscCall(PetscArrayzero(g1, NcS * NcT * dE)); 1048 PetscCall(PetscArrayzero(g2, NcS * NcT * dE)); 1049 PetscCall(PetscArrayzero(g3, NcS * NcT * dE * dE)); 1050 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 1051 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 1052 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 1053 for (e = 0; e < Ne; ++e) { 1054 PetscFEGeom fegeom; 1055 const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]}; 1056 const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]}; 1057 1058 fegeom.v = x; /* Workspace */ 1059 for (q = 0; q < Nq; ++q) { 1060 PetscInt qpt[2]; 1061 PetscReal w; 1062 PetscInt c; 1063 1064 PetscCall(PetscDSPermuteQuadPoint(ds, ornt[0], fieldI, q, &qpt[0])); 1065 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], 0), fieldI, q, &qpt[1])); 1066 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 1067 w = fegeom.detJ[0] * quadWeights[q]; 1068 if (debug > 1 && q < fgeom->numPoints) { 1069 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 1070 #if !defined(PETSC_USE_COMPLEX) 1071 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 1072 #endif 1073 } 1074 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 1075 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, face, qpt, TfIn, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 1076 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 1077 if (n0) { 1078 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1079 for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0); 1080 for (c = 0; c < NcS * NcT; ++c) g0[c] *= w; 1081 } 1082 if (n1) { 1083 PetscCall(PetscArrayzero(g1, NcS * NcT * dE)); 1084 for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1); 1085 for (c = 0; c < NcS * NcT * dE; ++c) g1[c] *= w; 1086 } 1087 if (n2) { 1088 PetscCall(PetscArrayzero(g2, NcS * NcT * dE)); 1089 for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2); 1090 for (c = 0; c < NcS * NcT * dE; ++c) g2[c] *= w; 1091 } 1092 if (n3) { 1093 PetscCall(PetscArrayzero(g3, NcS * NcT * dE * dE)); 1094 for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3); 1095 for (c = 0; c < NcS * NcT * dE * dE; ++c) g3[c] *= w; 1096 } 1097 1098 if (isCohesiveFieldI) { 1099 if (isCohesiveFieldJ) { 1100 PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1101 } else { 1102 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1103 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ], &g1[NcI * NcJ * dim], &g2[NcI * NcJ * dim], &g3[NcI * NcJ * dim * dim], eOffset, totDim, offsetI, offsetJ, elemMat)); 1104 } 1105 } else 1106 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, s, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1107 } 1108 if (debug > 1) { 1109 PetscInt fc, f, gc, g; 1110 1111 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 1112 for (fc = 0; fc < NcI; ++fc) { 1113 for (f = 0; f < T[fieldI]->Nb; ++f) { 1114 const PetscInt i = offsetI + f * NcI + fc; 1115 for (gc = 0; gc < NcJ; ++gc) { 1116 for (g = 0; g < T[fieldJ]->Nb; ++g) { 1117 const PetscInt j = offsetJ + g * NcJ + gc; 1118 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]))); 1119 } 1120 } 1121 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1122 } 1123 } 1124 } 1125 cOffset += totDim; 1126 cOffsetAux += totDimAux; 1127 eOffset += PetscSqr(totDim); 1128 } 1129 PetscFunctionReturn(PETSC_SUCCESS); 1130 } 1131 1132 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 1133 { 1134 PetscFunctionBegin; 1135 fem->ops->setfromoptions = NULL; 1136 fem->ops->setup = PetscFESetUp_Basic; 1137 fem->ops->view = PetscFEView_Basic; 1138 fem->ops->destroy = PetscFEDestroy_Basic; 1139 fem->ops->getdimension = PetscFEGetDimension_Basic; 1140 fem->ops->createtabulation = PetscFECreateTabulation_Basic; 1141 fem->ops->integrate = PetscFEIntegrate_Basic; 1142 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 1143 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 1144 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 1145 fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic; 1146 fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */; 1147 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 1148 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 1149 fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic; 1150 PetscFunctionReturn(PETSC_SUCCESS); 1151 } 1152 1153 /*MC 1154 PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization 1155 1156 Level: intermediate 1157 1158 .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()` 1159 M*/ 1160 1161 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 1162 { 1163 PetscFE_Basic *b; 1164 1165 PetscFunctionBegin; 1166 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1167 PetscCall(PetscNew(&b)); 1168 fem->data = b; 1169 1170 PetscCall(PetscFEInitialize_Basic(fem)); 1171 PetscFunctionReturn(PETSC_SUCCESS); 1172 } 1173