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 PETSC_INTERN 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 PETSC_INTERN 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 // LCOV_EXCL_START 432 PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " wt %g x:", q, (double)quadWeights[q])); 433 for (c = 0; c < dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)fegeom.v[c])); 434 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 435 if (debug > 2) { 436 PetscCall(PetscPrintf(PETSC_COMM_SELF, " field %" PetscInt_FMT ":", field)); 437 for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c]))); 438 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 439 PetscCall(PetscPrintf(PETSC_COMM_SELF, " field der %" PetscInt_FMT ":", field)); 440 for (c = 0; c < T[field]->Nc * dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u_x[uOff[field] + c]))); 441 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 442 PetscCall(PetscPrintf(PETSC_COMM_SELF, " resid %" PetscInt_FMT ":", field)); 443 for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c]))); 444 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 445 PetscCall(PetscPrintf(PETSC_COMM_SELF, " res der %" PetscInt_FMT ":", field)); 446 for (c = 0; c < T[field]->Nc; ++c) { 447 for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f1[(q * T[field]->Nc + c) * dim + d]))); 448 } 449 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 450 } 451 // LCOV_EXCL_STOP 452 } 453 } 454 PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset])); 455 cOffset += totDim; 456 cOffsetAux += totDimAux; 457 } 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460 461 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[]) 462 { 463 const PetscInt debug = 0; 464 const PetscInt field = key.field; 465 PetscFE fe; 466 PetscInt n0, n1, i; 467 PetscBdPointFunc *f0_func, *f1_func; 468 PetscQuadrature quad; 469 PetscTabulation *Tf, *TfAux = NULL; 470 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 471 const PetscScalar *constants; 472 PetscReal *x; 473 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 474 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI; 475 PetscBool auxOnBd = PETSC_FALSE; 476 const PetscReal *quadPoints, *quadWeights; 477 PetscInt qdim, qNc, Nq, q, dE; 478 479 PetscFunctionBegin; 480 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 481 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 482 PetscCall(PetscFEGetFaceQuadrature(fe, &quad)); 483 PetscCall(PetscDSGetNumFields(ds, &Nf)); 484 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 485 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 486 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 487 PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset)); 488 PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func)); 489 if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS); 490 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 491 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 492 PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL)); 493 PetscCall(PetscDSGetFaceTabulation(ds, &Tf)); 494 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 495 if (dsAux) { 496 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 497 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 498 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 499 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 500 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 501 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 502 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 503 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux)); 504 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux)); 505 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); 506 } 507 NcI = Tf[field]->Nc; 508 PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights)); 509 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 510 dE = fgeom->dimEmbed; 511 /* TODO FIX THIS */ 512 fgeom->dim = dim - 1; 513 PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim); 514 for (e = 0; e < Ne; ++e) { 515 PetscFEGeom fegeom, cgeom; 516 const PetscInt face = fgeom->face[e][0]; 517 518 fegeom.v = x; /* Workspace */ 519 PetscCall(PetscArrayzero(f0, Nq * NcI)); 520 PetscCall(PetscArrayzero(f1, Nq * NcI * dE)); 521 for (q = 0; q < Nq; ++q) { 522 PetscReal w; 523 PetscInt c, d; 524 525 PetscCall(PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom)); 526 PetscCall(PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom)); 527 w = fegeom.detJ[0] * quadWeights[q]; 528 if (debug > 1) { 529 if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) { 530 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 531 #if !defined(PETSC_USE_COMPLEX) 532 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 533 PetscCall(DMPrintCellVector(e, "n", dim, fegeom.n)); 534 #endif 535 } 536 } 537 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 538 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 539 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]); 540 for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w; 541 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]); 542 for (c = 0; c < NcI; ++c) 543 for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w; 544 if (debug) { 545 PetscCall(PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q)); 546 for (c = 0; c < NcI; ++c) { 547 if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c]))); 548 if (n1) { 549 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]))); 550 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 551 } 552 } 553 } 554 } 555 PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 556 cOffset += totDim; 557 cOffsetAux += totDimAux; 558 } 559 PetscFunctionReturn(PETSC_SUCCESS); 560 } 561 562 /* 563 BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but 564 all transforms operate in the full space and are square. 565 566 HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square. 567 1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces 568 2) We need to assume that the orientation is 0 for both 569 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() 570 */ 571 PETSC_INTERN 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[]) 572 { 573 const PetscInt debug = 0; 574 const PetscInt field = key.field; 575 PetscFE fe; 576 PetscWeakForm wf; 577 PetscInt n0, n1, i; 578 PetscBdPointFunc *f0_func, *f1_func; 579 PetscQuadrature quad; 580 DMPolytopeType ct; 581 PetscTabulation *Tf, *TfIn, *TfAux = NULL; 582 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 583 const PetscScalar *constants; 584 PetscReal *x; 585 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 586 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS; 587 PetscBool isCohesiveField, auxOnBd = PETSC_FALSE; 588 const PetscReal *quadPoints, *quadWeights; 589 PetscInt qdim, qNc, Nq, q, dE; 590 591 PetscFunctionBegin; 592 /* Hybrid discretization is posed directly on faces */ 593 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe)); 594 PetscCall(PetscFEGetSpatialDimension(fe, &dim)); 595 PetscCall(PetscFEGetQuadrature(fe, &quad)); 596 PetscCall(PetscDSGetNumFields(ds, &Nf)); 597 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 598 PetscCall(PetscDSGetTotalDimension(dsIn, &totDimIn)); 599 PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 0, &uOff)); // Change 0 to s for one-sided offsets 600 PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(dsIn, s, &uOff_x)); 601 PetscCall(PetscDSGetFieldOffsetCohesive(ds, field, &fOffset)); 602 PetscCall(PetscDSGetWeakForm(ds, &wf)); 603 PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func)); 604 if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS); 605 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 606 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL)); 607 PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL)); 608 /* NOTE This is a bulk tabulation because the DS is a face discretization */ 609 PetscCall(PetscDSGetTabulation(ds, &Tf)); 610 PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn)); 611 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 612 if (dsAux) { 613 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 614 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 615 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 616 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 617 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 618 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 619 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE; 620 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux)); 621 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux)); 622 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); 623 } 624 PetscCall(PetscDSGetCohesive(ds, field, &isCohesiveField)); 625 NcI = Tf[field]->Nc; 626 NcS = NcI; 627 if (!isCohesiveField && s == 2) { 628 // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides 629 NcS *= 2; 630 } 631 PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights)); 632 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 633 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 634 dE = fgeom->dimEmbed; 635 PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim); 636 for (e = 0; e < Ne; ++e) { 637 PetscFEGeom fegeom; 638 const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]}; 639 const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]}; 640 const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]}; 641 642 fegeom.v = x; /* Workspace */ 643 PetscCall(PetscArrayzero(f0, Nq * NcS)); 644 PetscCall(PetscArrayzero(f1, Nq * NcS * dE)); 645 for (q = 0; q < Nq; ++q) { 646 PetscInt qpt[2]; 647 PetscReal w; 648 PetscInt c, d; 649 650 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0])); 651 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), field, q, &qpt[1])); 652 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 653 w = fegeom.detJ[0] * quadWeights[q]; 654 if (debug > 1 && q < fgeom->numPoints) { 655 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 656 #if !defined(PETSC_USE_COMPLEX) 657 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ)); 658 #endif 659 } 660 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 661 /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */ 662 PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, &coefficients[cOffsetIn], &coefficients_t[cOffsetIn], u, u_x, u_t)); 663 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 664 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]); 665 for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w; 666 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]); 667 for (c = 0; c < NcS; ++c) 668 for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w; 669 } 670 if (isCohesiveField) { 671 PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 672 } else { 673 PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset])); 674 } 675 cOffset += totDim; 676 cOffsetIn += totDimIn; 677 cOffsetAux += totDimAux; 678 } 679 PetscFunctionReturn(PETSC_SUCCESS); 680 } 681 682 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[]) 683 { 684 const PetscInt debug = 0; 685 PetscFE feI, feJ; 686 PetscWeakForm wf; 687 PetscPointJac *g0_func, *g1_func, *g2_func, *g3_func; 688 PetscInt n0, n1, n2, n3, i; 689 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 690 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 691 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 692 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 693 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 694 PetscQuadrature quad; 695 PetscTabulation *T, *TAux = NULL; 696 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 697 const PetscScalar *constants; 698 PetscReal *x; 699 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 700 PetscInt NcI = 0, NcJ = 0; 701 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 702 PetscInt dE, Np; 703 PetscBool isAffine; 704 const PetscReal *quadPoints, *quadWeights; 705 PetscInt qNc, Nq, q; 706 707 PetscFunctionBegin; 708 PetscCall(PetscDSGetNumFields(ds, &Nf)); 709 fieldI = key.field / Nf; 710 fieldJ = key.field % Nf; 711 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 712 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 713 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 714 PetscCall(PetscFEGetQuadrature(feI, &quad)); 715 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 716 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 717 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 718 PetscCall(PetscDSGetWeakForm(ds, &wf)); 719 switch (jtype) { 720 case PETSCFE_JACOBIAN_DYN: 721 PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 722 break; 723 case PETSCFE_JACOBIAN_PRE: 724 PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 725 break; 726 case PETSCFE_JACOBIAN: 727 PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 728 break; 729 } 730 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 731 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 732 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 733 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 734 PetscCall(PetscDSGetTabulation(ds, &T)); 735 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 736 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 737 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 738 if (dsAux) { 739 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 740 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 741 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 742 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 743 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 744 PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 745 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); 746 } 747 NcI = T[fieldI]->Nc; 748 NcJ = T[fieldJ]->Nc; 749 Np = cgeom->numPoints; 750 dE = cgeom->dimEmbed; 751 isAffine = cgeom->isAffine; 752 /* Initialize here in case the function is not defined */ 753 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 754 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 755 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 756 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 757 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 758 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 759 for (e = 0; e < Ne; ++e) { 760 PetscFEGeom fegeom; 761 762 fegeom.dim = cgeom->dim; 763 fegeom.dimEmbed = cgeom->dimEmbed; 764 if (isAffine) { 765 fegeom.v = x; 766 fegeom.xi = cgeom->xi; 767 fegeom.J = &cgeom->J[e * Np * dE * dE]; 768 fegeom.invJ = &cgeom->invJ[e * Np * dE * dE]; 769 fegeom.detJ = &cgeom->detJ[e * Np]; 770 } 771 for (q = 0; q < Nq; ++q) { 772 PetscReal w; 773 PetscInt c; 774 775 if (isAffine) { 776 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x); 777 } else { 778 fegeom.v = &cgeom->v[(e * Np + q) * dE]; 779 fegeom.J = &cgeom->J[(e * Np + q) * dE * dE]; 780 fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE]; 781 fegeom.detJ = &cgeom->detJ[e * Np + q]; 782 } 783 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0])); 784 w = fegeom.detJ[0] * quadWeights[q]; 785 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 786 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 787 if (n0) { 788 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 789 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); 790 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 791 } 792 if (n1) { 793 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 794 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); 795 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 796 } 797 if (n2) { 798 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 799 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); 800 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 801 } 802 if (n3) { 803 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 804 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); 805 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 806 } 807 808 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)); 809 } 810 if (debug > 1) { 811 PetscInt fc, f, gc, g; 812 813 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 814 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 815 for (f = 0; f < T[fieldI]->Nb; ++f) { 816 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 817 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 818 for (g = 0; g < T[fieldJ]->Nb; ++g) { 819 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 820 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]))); 821 } 822 } 823 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 824 } 825 } 826 } 827 cOffset += totDim; 828 cOffsetAux += totDimAux; 829 eOffset += PetscSqr(totDim); 830 } 831 PetscFunctionReturn(PETSC_SUCCESS); 832 } 833 834 PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFEJacobianType jtype, 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[]) 835 { 836 const PetscInt debug = 0; 837 PetscFE feI, feJ; 838 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 839 PetscInt n0, n1, n2, n3, i; 840 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 841 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 842 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 843 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 844 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 845 PetscQuadrature quad; 846 PetscTabulation *T, *TAux = NULL; 847 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 848 const PetscScalar *constants; 849 PetscReal *x; 850 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 851 PetscInt NcI = 0, NcJ = 0; 852 PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 853 PetscBool isAffine; 854 const PetscReal *quadPoints, *quadWeights; 855 PetscInt qNc, Nq, q, Np, dE; 856 857 PetscFunctionBegin; 858 PetscCall(PetscDSGetNumFields(ds, &Nf)); 859 fieldI = key.field / Nf; 860 fieldJ = key.field % Nf; 861 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 862 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 863 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 864 PetscCall(PetscFEGetFaceQuadrature(feI, &quad)); 865 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 866 PetscCall(PetscDSGetComponentOffsets(ds, &uOff)); 867 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x)); 868 PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI)); 869 PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ)); 870 switch (jtype) { 871 case PETSCFE_JACOBIAN_PRE: 872 PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 873 break; 874 case PETSCFE_JACOBIAN: 875 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 876 break; 877 case PETSCFE_JACOBIAN_DYN: 878 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSCFE_JACOBIAN_DYN is not supported for PetscFEIntegrateBdJacobian()"); 879 } 880 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 881 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 882 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 883 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 884 PetscCall(PetscDSGetFaceTabulation(ds, &T)); 885 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 886 if (dsAux) { 887 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 888 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 889 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 890 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 891 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 892 PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 893 } 894 NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc; 895 Np = fgeom->numPoints; 896 dE = fgeom->dimEmbed; 897 isAffine = fgeom->isAffine; 898 /* Initialize here in case the function is not defined */ 899 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 900 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 901 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 902 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 903 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 904 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 905 for (e = 0; e < Ne; ++e) { 906 PetscFEGeom fegeom, cgeom; 907 const PetscInt face = fgeom->face[e][0]; 908 fegeom.n = NULL; 909 fegeom.v = NULL; 910 fegeom.J = NULL; 911 fegeom.detJ = NULL; 912 fegeom.dim = fgeom->dim; 913 fegeom.dimEmbed = fgeom->dimEmbed; 914 cgeom.dim = fgeom->dim; 915 cgeom.dimEmbed = fgeom->dimEmbed; 916 if (isAffine) { 917 fegeom.v = x; 918 fegeom.xi = fgeom->xi; 919 fegeom.J = &fgeom->J[e * Np * dE * dE]; 920 fegeom.invJ = &fgeom->invJ[e * Np * dE * dE]; 921 fegeom.detJ = &fgeom->detJ[e * Np]; 922 fegeom.n = &fgeom->n[e * Np * dE]; 923 924 cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE]; 925 cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE]; 926 cgeom.detJ = &fgeom->suppDetJ[0][e * Np]; 927 } 928 for (q = 0; q < Nq; ++q) { 929 PetscReal w; 930 PetscInt c; 931 932 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 933 if (isAffine) { 934 CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x); 935 } else { 936 fegeom.v = &fgeom->v[(e * Np + q) * dE]; 937 fegeom.J = &fgeom->J[(e * Np + q) * dE * dE]; 938 fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE]; 939 fegeom.detJ = &fgeom->detJ[e * Np + q]; 940 fegeom.n = &fgeom->n[(e * Np + q) * dE]; 941 942 cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE]; 943 cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE]; 944 cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q]; 945 } 946 w = fegeom.detJ[0] * quadWeights[q]; 947 if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t)); 948 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 949 if (n0) { 950 PetscCall(PetscArrayzero(g0, NcI * NcJ)); 951 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); 952 for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w; 953 } 954 if (n1) { 955 PetscCall(PetscArrayzero(g1, NcI * NcJ * dE)); 956 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); 957 for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w; 958 } 959 if (n2) { 960 PetscCall(PetscArrayzero(g2, NcI * NcJ * dE)); 961 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); 962 for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w; 963 } 964 if (n3) { 965 PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE)); 966 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); 967 for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w; 968 } 969 970 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)); 971 } 972 if (debug > 1) { 973 PetscInt fc, f, gc, g; 974 975 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ)); 976 for (fc = 0; fc < T[fieldI]->Nc; ++fc) { 977 for (f = 0; f < T[fieldI]->Nb; ++f) { 978 const PetscInt i = offsetI + f * T[fieldI]->Nc + fc; 979 for (gc = 0; gc < T[fieldJ]->Nc; ++gc) { 980 for (g = 0; g < T[fieldJ]->Nb; ++g) { 981 const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc; 982 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]))); 983 } 984 } 985 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 986 } 987 } 988 } 989 cOffset += totDim; 990 cOffsetAux += totDimAux; 991 eOffset += PetscSqr(totDim); 992 } 993 PetscFunctionReturn(PETSC_SUCCESS); 994 } 995 996 PETSC_INTERN 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[]) 997 { 998 const PetscInt debug = 0; 999 PetscFE feI, feJ; 1000 PetscWeakForm wf; 1001 PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func; 1002 PetscInt n0, n1, n2, n3, i; 1003 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 1004 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 1005 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 1006 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 1007 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 1008 PetscQuadrature quad; 1009 DMPolytopeType ct; 1010 PetscTabulation *T, *TfIn, *TAux = NULL; 1011 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 1012 const PetscScalar *constants; 1013 PetscReal *x; 1014 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL; 1015 PetscInt NcI = 0, NcJ = 0, NcS, NcT; 1016 PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e; 1017 PetscBool isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE; 1018 const PetscReal *quadPoints, *quadWeights; 1019 PetscInt qNc, Nq, q; 1020 1021 PetscFunctionBegin; 1022 PetscCall(PetscDSGetNumFields(ds, &Nf)); 1023 fieldI = key.field / Nf; 1024 fieldJ = key.field % Nf; 1025 /* Hybrid discretization is posed directly on faces */ 1026 PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI)); 1027 PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ)); 1028 PetscCall(PetscFEGetSpatialDimension(feI, &dim)); 1029 PetscCall(PetscFEGetQuadrature(feI, &quad)); 1030 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 1031 PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets 1032 PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x)); 1033 PetscCall(PetscDSGetWeakForm(ds, &wf)); 1034 switch (jtype) { 1035 case PETSCFE_JACOBIAN_PRE: 1036 PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1037 break; 1038 case PETSCFE_JACOBIAN: 1039 PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func)); 1040 break; 1041 case PETSCFE_JACOBIAN_DYN: 1042 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)"); 1043 } 1044 if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS); 1045 PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x)); 1046 PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal)); 1047 PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3)); 1048 PetscCall(PetscDSGetTabulation(ds, &T)); 1049 PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn)); 1050 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI)); 1051 PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ)); 1052 PetscCall(PetscDSGetConstants(ds, &numConstants, &constants)); 1053 if (dsAux) { 1054 PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux)); 1055 PetscCall(PetscDSGetNumFields(dsAux, &NfAux)); 1056 PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux)); 1057 PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff)); 1058 PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x)); 1059 PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x)); 1060 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE; 1061 if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux)); 1062 else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux)); 1063 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); 1064 } 1065 PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI)); 1066 PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ)); 1067 NcI = T[fieldI]->Nc; 1068 NcJ = T[fieldJ]->Nc; 1069 NcS = isCohesiveFieldI ? NcI : 2 * NcI; 1070 NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ; 1071 if (!isCohesiveFieldI && s == 2) { 1072 // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides 1073 NcS *= 2; 1074 } 1075 if (!isCohesiveFieldJ && s == 2) { 1076 // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides 1077 NcT *= 2; 1078 } 1079 // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though 1080 // the coordinates are in dE dimensions 1081 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1082 PetscCall(PetscArrayzero(g1, NcS * NcT * dim)); 1083 PetscCall(PetscArrayzero(g2, NcS * NcT * dim)); 1084 PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim)); 1085 PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights)); 1086 PetscCall(PetscQuadratureGetCellType(quad, &ct)); 1087 PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc); 1088 for (e = 0; e < Ne; ++e) { 1089 PetscFEGeom fegeom; 1090 const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]}; 1091 const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]}; 1092 const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]}; 1093 1094 fegeom.v = x; /* Workspace */ 1095 for (q = 0; q < Nq; ++q) { 1096 PetscInt qpt[2]; 1097 PetscReal w; 1098 PetscInt c; 1099 1100 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0])); 1101 PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1])); 1102 PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom)); 1103 w = fegeom.detJ[0] * quadWeights[q]; 1104 if (debug > 1 && q < fgeom->numPoints) { 1105 PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0])); 1106 #if !defined(PETSC_USE_COMPLEX) 1107 PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ)); 1108 #endif 1109 } 1110 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q)); 1111 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)); 1112 if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL)); 1113 if (n0) { 1114 PetscCall(PetscArrayzero(g0, NcS * NcT)); 1115 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); 1116 for (c = 0; c < NcS * NcT; ++c) g0[c] *= w; 1117 } 1118 if (n1) { 1119 PetscCall(PetscArrayzero(g1, NcS * NcT * dim)); 1120 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); 1121 for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w; 1122 } 1123 if (n2) { 1124 PetscCall(PetscArrayzero(g2, NcS * NcT * dim)); 1125 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); 1126 for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w; 1127 } 1128 if (n3) { 1129 PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim)); 1130 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); 1131 for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w; 1132 } 1133 1134 if (isCohesiveFieldI) { 1135 if (isCohesiveFieldJ) { 1136 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)); 1137 } else { 1138 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1139 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 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)); 1140 } 1141 } else { 1142 if (s == 2) { 1143 if (isCohesiveFieldJ) { 1144 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1145 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 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)); 1146 } else { 1147 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1148 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 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)); 1149 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ * 2], &g1[NcI * NcJ * dim * 2], &g2[NcI * NcJ * dim * 2], &g3[NcI * NcJ * dim * dim * 2], eOffset, totDim, offsetI, offsetJ, elemMat)); 1150 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ * 3], &g1[NcI * NcJ * dim * 3], &g2[NcI * NcJ * dim * 3], &g3[NcI * NcJ * dim * dim * 3], eOffset, totDim, offsetI, offsetJ, elemMat)); 1151 } 1152 } else 1153 PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, s, s, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat)); 1154 } 1155 } 1156 if (debug > 1) { 1157 const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb)); 1158 const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb)); 1159 const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb)); 1160 const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb)); 1161 PetscInt f, g; 1162 1163 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)); 1164 for (f = fS; f < fE; ++f) { 1165 const PetscInt i = offsetI + f; 1166 for (g = gS; g < gE; ++g) { 1167 const PetscInt j = offsetJ + g; 1168 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); 1169 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]))); 1170 } 1171 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n")); 1172 } 1173 } 1174 cOffset += totDim; 1175 cOffsetAux += totDimAux; 1176 eOffset += PetscSqr(totDim); 1177 } 1178 PetscFunctionReturn(PETSC_SUCCESS); 1179 } 1180 1181 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 1182 { 1183 PetscFunctionBegin; 1184 fem->ops->setfromoptions = NULL; 1185 fem->ops->setup = PetscFESetUp_Basic; 1186 fem->ops->view = PetscFEView_Basic; 1187 fem->ops->destroy = PetscFEDestroy_Basic; 1188 fem->ops->getdimension = PetscFEGetDimension_Basic; 1189 fem->ops->createtabulation = PetscFECreateTabulation_Basic; 1190 fem->ops->integrate = PetscFEIntegrate_Basic; 1191 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 1192 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 1193 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 1194 fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic; 1195 fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */; 1196 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 1197 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 1198 fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic; 1199 PetscFunctionReturn(PETSC_SUCCESS); 1200 } 1201 1202 /*MC 1203 PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization 1204 1205 Level: intermediate 1206 1207 .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()` 1208 M*/ 1209 1210 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 1211 { 1212 PetscFE_Basic *b; 1213 1214 PetscFunctionBegin; 1215 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 1216 PetscCall(PetscNew(&b)); 1217 fem->data = b; 1218 1219 PetscCall(PetscFEInitialize_Basic(fem)); 1220 PetscFunctionReturn(PETSC_SUCCESS); 1221 } 1222