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, 0, &uOff)); // Change 0 to s for one-sided offsets 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 const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]}; 625 626 fegeom.v = x; /* Workspace */ 627 PetscCall(PetscArrayzero(f0, Nq * NcS)); 628 PetscCall(PetscArrayzero(f1, Nq * NcS * dE)); 629 for (q = 0; q < Nq; ++q) { 630 PetscInt qpt[2]; 631 PetscReal w; 632 PetscInt c, d; 633 634 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0])); 635 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), field, q, &qpt[1])); 636 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 637 w = fegeom.detJ[0] * quadWeights[q]; 638 if (debug > 1 && q < fgeom->numPoints) { 639 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 640 #if !defined(PETSC_USE_COMPLEX) 641 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ)); 642 #endif 643 } 644 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 645 /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */ 646 PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, &coefficients[cOffsetIn], &coefficients_t[cOffsetIn], u, u_x, u_t)); 647 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 648 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]); 649 for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w; 650 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]); 651 for (c = 0; c < NcS; ++c) 652 for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w; 653 } 654 if (isCohesiveField) { 655 PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 656 } else { 657 PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 658 } 659 cOffset += totDim; 660 cOffsetIn += totDimIn; 661 cOffsetAux += totDimAux; 662 } 663 PetscFunctionReturn(PETSC_SUCCESS); 664 } 665 666 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[]) 667 { 668 const PetscInt debug = 0; 669 PetscFE feI, feJ; 670 PetscWeakForm wf; 671 PetscPointJac *g0_func, *g1_func, *g2_func, *g3_func; 672 PetscInt n0, n1, n2, n3, i; 673 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 674 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 675 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 676 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 677 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 678 PetscQuadrature quad; 679 PetscTabulation *T, *TAux = NULL; 680 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 681 const PetscScalar *constants; 682 PetscReal *x; 683 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 684 PetscInt NcI = 0, NcJ = 0; 685 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 686 PetscInt dE, Np; 687 PetscBool isAffine; 688 const PetscReal *quadPoints, *quadWeights; 689 PetscInt qNc, Nq, q; 690 691 PetscFunctionBegin; 692 PetscCall(PetscDSGetNumFields(ds, &Nf)); 693 fieldI = key.field / Nf; 694 fieldJ = key.field % Nf; 695 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 696 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 697 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 698 PetscCall(PetscFEGetQuadrature(feI, &quad)); 699 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 700 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 701 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 702 PetscCall(PetscDSGetWeakForm(ds, &wf)); 703 switch (jtype) { 704 case PETSCFE_JACOBIAN_DYN: 705 PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 706 break; 707 case PETSCFE_JACOBIAN_PRE: 708 PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 709 break; 710 case PETSCFE_JACOBIAN: 711 PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 712 break; 713 } 714 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 715 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 716 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 717 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 718 PetscCall(PetscDSGetTabulation(ds, &T)); 719 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 720 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 721 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 722 if (dsAux) { 723 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 724 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 725 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 726 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 727 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 728 PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 729 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); 730 } 731 NcI = T[fieldI]->Nc; 732 NcJ = T[fieldJ]->Nc; 733 Np = cgeom->numPoints; 734 dE = cgeom->dimEmbed; 735 isAffine = cgeom->isAffine; 736 /* Initialize here in case the function is not defined */ 737 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 738 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 739 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 740 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 741 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 742 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 743 for (e = 0; e < Ne; ++e) { 744 PetscFEGeom fegeom; 745 746 fegeom.dim = cgeom->dim; 747 fegeom.dimEmbed = cgeom->dimEmbed; 748 if (isAffine) { 749 fegeom.v = x; 750 fegeom.xi = cgeom->xi; 751 fegeom.J = &cgeom->J[e * Np * dE * dE]; 752 fegeom.invJ = &cgeom->invJ[e * Np * dE * dE]; 753 fegeom.detJ = &cgeom->detJ[e * Np]; 754 } 755 for (q = 0; q < Nq; ++q) { 756 PetscReal w; 757 PetscInt c; 758 759 if (isAffine) { 760 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x); 761 } else { 762 fegeom.v = &cgeom->v[(e * Np + q) * dE]; 763 fegeom.J = &cgeom->J[(e * Np + q) * dE * dE]; 764 fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE]; 765 fegeom.detJ = &cgeom->detJ[e * Np + q]; 766 } 767 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 768 w = fegeom.detJ[0] * quadWeights[q]; 769 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 770 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 771 if (n0) { 772 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 773 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); 774 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 775 } 776 if (n1) { 777 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 778 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); 779 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 780 } 781 if (n2) { 782 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 783 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); 784 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 785 } 786 if (n3) { 787 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 788 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); 789 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 790 } 791 792 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)); 793 } 794 if (debug > 1) { 795 PetscInt fc, f, gc, g; 796 797 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 798 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 799 for (f = 0; f < T[fieldI]->Nb; ++f) { 800 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 801 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 802 for (g = 0; g < T[fieldJ]->Nb; ++g) { 803 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 804 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]))); 805 } 806 } 807 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 808 } 809 } 810 } 811 cOffset += totDim; 812 cOffsetAux += totDimAux; 813 eOffset += PetscSqr(totDim); 814 } 815 PetscFunctionReturn(PETSC_SUCCESS); 816 } 817 818 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[]) 819 { 820 const PetscInt debug = 0; 821 PetscFE feI, feJ; 822 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 823 PetscInt n0, n1, n2, n3, i; 824 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 825 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 826 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 827 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 828 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 829 PetscQuadrature quad; 830 PetscTabulation *T, *TAux = NULL; 831 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 832 const PetscScalar *constants; 833 PetscReal *x; 834 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 835 PetscInt NcI = 0, NcJ = 0; 836 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 837 PetscBool isAffine; 838 const PetscReal *quadPoints, *quadWeights; 839 PetscInt qNc, Nq, q, Np, dE; 840 841 PetscFunctionBegin; 842 PetscCall(PetscDSGetNumFields(ds, &Nf)); 843 fieldI = key.field / Nf; 844 fieldJ = key.field % Nf; 845 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 846 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 847 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 848 PetscCall(PetscFEGetFaceQuadrature(feI, &quad)); 849 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 850 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 851 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 852 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 853 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 854 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 855 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 856 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 857 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 858 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 859 PetscCall(PetscDSGetFaceTabulation(ds, &T)); 860 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 861 if (dsAux) { 862 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 863 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 864 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 865 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 866 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 867 PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 868 } 869 NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc; 870 Np = fgeom->numPoints; 871 dE = fgeom->dimEmbed; 872 isAffine = fgeom->isAffine; 873 /* Initialize here in case the function is not defined */ 874 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 875 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 876 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 877 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 878 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 879 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 880 for (e = 0; e < Ne; ++e) { 881 PetscFEGeom fegeom, cgeom; 882 const PetscInt face = fgeom->face[e][0]; 883 fegeom.n = NULL; 884 fegeom.v = NULL; 885 fegeom.J = NULL; 886 fegeom.detJ = NULL; 887 fegeom.dim = fgeom->dim; 888 fegeom.dimEmbed = fgeom->dimEmbed; 889 cgeom.dim = fgeom->dim; 890 cgeom.dimEmbed = fgeom->dimEmbed; 891 if (isAffine) { 892 fegeom.v = x; 893 fegeom.xi = fgeom->xi; 894 fegeom.J = &fgeom->J[e * Np * dE * dE]; 895 fegeom.invJ = &fgeom->invJ[e * Np * dE * dE]; 896 fegeom.detJ = &fgeom->detJ[e * Np]; 897 fegeom.n = &fgeom->n[e * Np * dE]; 898 899 cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE]; 900 cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE]; 901 cgeom.detJ = &fgeom->suppDetJ[0][e * Np]; 902 } 903 for (q = 0; q < Nq; ++q) { 904 PetscReal w; 905 PetscInt c; 906 907 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 908 if (isAffine) { 909 CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x); 910 } else { 911 fegeom.v = &fgeom->v[(e * Np + q) * dE]; 912 fegeom.J = &fgeom->J[(e * Np + q) * dE * dE]; 913 fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE]; 914 fegeom.detJ = &fgeom->detJ[e * Np + q]; 915 fegeom.n = &fgeom->n[(e * Np + q) * dE]; 916 917 cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE]; 918 cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE]; 919 cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q]; 920 } 921 w = fegeom.detJ[0] * quadWeights[q]; 922 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 923 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 924 if (n0) { 925 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 926 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); 927 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 928 } 929 if (n1) { 930 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 931 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); 932 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 933 } 934 if (n2) { 935 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 936 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); 937 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 938 } 939 if (n3) { 940 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 941 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); 942 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 943 } 944 945 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)); 946 } 947 if (debug > 1) { 948 PetscInt fc, f, gc, g; 949 950 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 951 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 952 for (f = 0; f < T[fieldI]->Nb; ++f) { 953 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 954 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 955 for (g = 0; g < T[fieldJ]->Nb; ++g) { 956 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 957 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]))); 958 } 959 } 960 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 961 } 962 } 963 } 964 cOffset += totDim; 965 cOffsetAux += totDimAux; 966 eOffset += PetscSqr(totDim); 967 } 968 PetscFunctionReturn(PETSC_SUCCESS); 969 } 970 971 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[]) 972 { 973 const PetscInt debug = 0; 974 PetscFE feI, feJ; 975 PetscWeakForm wf; 976 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 977 PetscInt n0, n1, n2, n3, i; 978 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 979 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 980 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 981 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 982 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 983 PetscQuadrature quad; 984 DMPolytopeType ct; 985 PetscTabulation *T, *TfIn, *TAux = NULL; 986 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 987 const PetscScalar *constants; 988 PetscReal *x; 989 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 990 PetscInt NcI = 0, NcJ = 0, NcS, NcT; 991 PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 992 PetscBool isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE; 993 const PetscReal *quadPoints, *quadWeights; 994 PetscInt qNc, Nq, q; 995 996 PetscFunctionBegin; 997 PetscCall(PetscDSGetNumFields(ds, &Nf)); 998 fieldI = key.field / Nf; 999 fieldJ = key.field % Nf; 1000 /* Hybrid discretization is posed directly on faces */ 1001 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 1002 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 1003 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 1004 PetscCall(PetscFEGetQuadrature(feI, &quad)); 1005 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 1006 PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets 1007 PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x)); 1008 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1009 switch (jtype) { 1010 case PETSCFE_JACOBIAN_PRE: 1011 PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1012 break; 1013 case PETSCFE_JACOBIAN: 1014 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1015 break; 1016 case PETSCFE_JACOBIAN_DYN: 1017 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)"); 1018 } 1019 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 1020 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 1021 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 1022 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 1023 PetscCall(PetscDSGetTabulation(ds, &T)); 1024 PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn)); 1025 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI)); 1026 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ)); 1027 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 1028 if (dsAux) { 1029 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 1030 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 1031 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 1032 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 1033 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 1034 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 1035 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE; 1036 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 1037 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 1038 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); 1039 } 1040 PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI)); 1041 PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ)); 1042 NcI = T[fieldI]->Nc; 1043 NcJ = T[fieldJ]->Nc; 1044 NcS = isCohesiveFieldI ? NcI : 2 * NcI; 1045 NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ; 1046 // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though 1047 // the coordinates are in dE dimensions 1048 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1049 PetscCall(PetscArrayzero(g1, NcS * NcT * dim)); 1050 PetscCall(PetscArrayzero(g2, NcS * NcT * dim)); 1051 PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim)); 1052 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 1053 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 1054 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 1055 for (e = 0; e < Ne; ++e) { 1056 PetscFEGeom fegeom; 1057 const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]}; 1058 const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]}; 1059 const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]}; 1060 1061 fegeom.v = x; /* Workspace */ 1062 for (q = 0; q < Nq; ++q) { 1063 PetscInt qpt[2]; 1064 PetscReal w; 1065 PetscInt c; 1066 1067 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0])); 1068 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1])); 1069 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 1070 w = fegeom.detJ[0] * quadWeights[q]; 1071 if (debug > 1 && q < fgeom->numPoints) { 1072 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 1073 #if !defined(PETSC_USE_COMPLEX) 1074 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 1075 #endif 1076 } 1077 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 1078 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)); 1079 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 1080 if (n0) { 1081 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1082 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); 1083 for (c = 0; c < NcS * NcT; ++c) g0[c] *= w; 1084 } 1085 if (n1) { 1086 PetscCall(PetscArrayzero(g1, NcS * NcT * dim)); 1087 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); 1088 for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w; 1089 } 1090 if (n2) { 1091 PetscCall(PetscArrayzero(g2, NcS * NcT * dim)); 1092 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); 1093 for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w; 1094 } 1095 if (n3) { 1096 PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim)); 1097 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); 1098 for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w; 1099 } 1100 1101 if (isCohesiveFieldI) { 1102 if (isCohesiveFieldJ) { 1103 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)); 1104 } else { 1105 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)); 1106 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)); 1107 } 1108 } else 1109 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)); 1110 } 1111 if (debug > 1) { 1112 const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb)); 1113 const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb)); 1114 const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb)); 1115 const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb)); 1116 PetscInt f, g; 1117 1118 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT " s %s totDim %" PetscInt_FMT " offsets (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", fieldI, fieldJ, s ? (s > 1 ? "Coh" : "Pos") : "Neg", totDim, eOffset, offsetI, offsetJ)); 1119 for (f = fS; f < fE; ++f) { 1120 const PetscInt i = offsetI + f; 1121 for (g = gS; g < gE; ++g) { 1122 const PetscInt j = offsetJ + g; 1123 PetscCheck(i < totDim && j < totDim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fuck up %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", f, i, g, j); 1124 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f / NcI, f % NcI, g / NcJ, g % NcJ, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]))); 1125 } 1126 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1127 } 1128 } 1129 cOffset += totDim; 1130 cOffsetAux += totDimAux; 1131 eOffset += PetscSqr(totDim); 1132 } 1133 PetscFunctionReturn(PETSC_SUCCESS); 1134 } 1135 1136 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 1137 { 1138 PetscFunctionBegin; 1139 fem->ops->setfromoptions = NULL; 1140 fem->ops->setup = PetscFESetUp_Basic; 1141 fem->ops->view = PetscFEView_Basic; 1142 fem->ops->destroy = PetscFEDestroy_Basic; 1143 fem->ops->getdimension = PetscFEGetDimension_Basic; 1144 fem->ops->createtabulation = PetscFECreateTabulation_Basic; 1145 fem->ops->integrate = PetscFEIntegrate_Basic; 1146 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 1147 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 1148 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 1149 fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic; 1150 fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */; 1151 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 1152 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 1153 fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic; 1154 PetscFunctionReturn(PETSC_SUCCESS); 1155 } 1156 1157 /*MC 1158 PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization 1159 1160 Level: intermediate 1161 1162 .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()` 1163 M*/ 1164 1165 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 1166 { 1167 PetscFE_Basic *b; 1168 1169 PetscFunctionBegin; 1170 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1171 PetscCall(PetscNew(&b)); 1172 fem->data = b; 1173 1174 PetscCall(PetscFEInitialize_Basic(fem)); 1175 PetscFunctionReturn(PETSC_SUCCESS); 1176 } 1177