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