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 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = PetscFree(b);CHKERRQ(ierr); 11 PetscFunctionReturn(0); 12 } 13 14 static PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer v) 15 { 16 PetscInt dim, Nc; 17 PetscSpace basis = NULL; 18 PetscDualSpace dual = NULL; 19 PetscQuadrature quad = NULL; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 24 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 25 ierr = PetscFEGetBasisSpace(fe, &basis);CHKERRQ(ierr); 26 ierr = PetscFEGetDualSpace(fe, &dual);CHKERRQ(ierr); 27 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 28 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 29 ierr = PetscViewerASCIIPrintf(v, "Basic Finite Element in %D dimensions with %D components\n",dim,Nc);CHKERRQ(ierr); 30 if (basis) {ierr = PetscSpaceView(basis, v);CHKERRQ(ierr);} 31 if (dual) {ierr = PetscDualSpaceView(dual, v);CHKERRQ(ierr);} 32 if (quad) {ierr = PetscQuadratureView(quad, v);CHKERRQ(ierr);} 33 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 static PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v) 38 { 39 PetscBool iascii; 40 PetscErrorCode ierr; 41 42 PetscFunctionBegin; 43 ierr = PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 44 if (iascii) {ierr = PetscFEView_Basic_Ascii(fe, v);CHKERRQ(ierr);} 45 PetscFunctionReturn(0); 46 } 47 48 /* Construct the change of basis from prime basis to nodal basis */ 49 PetscErrorCode PetscFESetUp_Basic(PetscFE fem) 50 { 51 PetscScalar *work, *invVscalar; 52 PetscBLASInt *pivots; 53 PetscBLASInt n, info; 54 PetscInt pdim, j; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 59 ierr = PetscMalloc1(pdim*pdim,&fem->invV);CHKERRQ(ierr); 60 #if defined(PETSC_USE_COMPLEX) 61 ierr = PetscMalloc1(pdim*pdim,&invVscalar);CHKERRQ(ierr); 62 #else 63 invVscalar = fem->invV; 64 #endif 65 for (j = 0; j < pdim; ++j) { 66 PetscReal *Bf; 67 PetscQuadrature f; 68 const PetscReal *points, *weights; 69 PetscInt Nc, Nq, q, k, c; 70 71 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr); 72 ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr); 73 ierr = PetscMalloc1(Nc*Nq*pdim,&Bf);CHKERRQ(ierr); 74 ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr); 75 for (k = 0; k < pdim; ++k) { 76 /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */ 77 invVscalar[j*pdim+k] = 0.0; 78 79 for (q = 0; q < Nq; ++q) { 80 for (c = 0; c < Nc; ++c) invVscalar[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c]; 81 } 82 } 83 ierr = PetscFree(Bf);CHKERRQ(ierr); 84 } 85 86 ierr = PetscMalloc2(pdim,&pivots,pdim,&work);CHKERRQ(ierr); 87 n = pdim; 88 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info)); 89 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info); 90 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info)); 91 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info); 92 #if defined(PETSC_USE_COMPLEX) 93 for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]); 94 ierr = PetscFree(invVscalar);CHKERRQ(ierr); 95 #endif 96 ierr = PetscFree2(pivots,work);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim) 101 { 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 ierr = PetscDualSpaceGetDimension(fem->dualSpace, dim);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H) 110 { 111 DM dm; 112 PetscInt pdim; /* Dimension of FE space P */ 113 PetscInt dim; /* Spatial dimension */ 114 PetscInt Nc; /* Field components */ 115 PetscReal *tmpB, *tmpD, *tmpH; 116 PetscInt p, d, j, k, c; 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 121 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 122 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 123 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr); 124 /* Evaluate the prime basis functions at all points */ 125 if (B) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 126 if (D) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 127 if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 128 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);CHKERRQ(ierr); 129 /* Translate to the nodal basis */ 130 for (p = 0; p < npoints; ++p) { 131 if (B) { 132 /* Multiply by V^{-1} (pdim x pdim) */ 133 for (j = 0; j < pdim; ++j) { 134 const PetscInt i = (p*pdim + j)*Nc; 135 136 for (c = 0; c < Nc; ++c) { 137 B[i+c] = 0.0; 138 for (k = 0; k < pdim; ++k) { 139 B[i+c] += fem->invV[k*pdim+j] * tmpB[(p*pdim + k)*Nc+c]; 140 } 141 } 142 } 143 } 144 if (D) { 145 /* Multiply by V^{-1} (pdim x pdim) */ 146 for (j = 0; j < pdim; ++j) { 147 for (c = 0; c < Nc; ++c) { 148 for (d = 0; d < dim; ++d) { 149 const PetscInt i = ((p*pdim + j)*Nc + c)*dim + d; 150 151 D[i] = 0.0; 152 for (k = 0; k < pdim; ++k) { 153 D[i] += fem->invV[k*pdim+j] * tmpD[((p*pdim + k)*Nc + c)*dim + d]; 154 } 155 } 156 } 157 } 158 } 159 if (H) { 160 /* Multiply by V^{-1} (pdim x pdim) */ 161 for (j = 0; j < pdim; ++j) { 162 for (c = 0; c < Nc; ++c) { 163 for (d = 0; d < dim*dim; ++d) { 164 const PetscInt i = ((p*pdim + j)*Nc + c)*dim*dim + d; 165 166 H[i] = 0.0; 167 for (k = 0; k < pdim; ++k) { 168 H[i] += fem->invV[k*pdim+j] * tmpH[((p*pdim + k)*Nc + c)*dim*dim + d]; 169 } 170 } 171 } 172 } 173 } 174 } 175 if (B) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 176 if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 177 if (H) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 178 PetscFunctionReturn(0); 179 } 180 181 static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 182 const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 183 { 184 const PetscInt debug = 0; 185 PetscFE fe; 186 PetscPointFunc obj_func; 187 PetscQuadrature quad; 188 PetscScalar *u, *u_x, *a, *a_x; 189 const PetscScalar *constants; 190 PetscReal *x; 191 PetscReal **B, **D, **BAux = NULL, **DAux = NULL; 192 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 193 PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 194 PetscBool isAffine; 195 const PetscReal *quadPoints, *quadWeights; 196 PetscInt qNc, Nq, q; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 ierr = PetscDSGetObjective(ds, field, &obj_func);CHKERRQ(ierr); 201 if (!obj_func) PetscFunctionReturn(0); 202 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 203 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 204 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 205 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 206 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 207 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 208 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 209 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 210 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 211 ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr); 212 ierr = PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 213 ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr); 214 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 215 if (dsAux) { 216 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 217 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 218 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 219 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 220 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 221 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 222 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 223 ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 224 } 225 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 226 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 227 Np = cgeom->numPoints; 228 dE = cgeom->dimEmbed; 229 isAffine = cgeom->isAffine; 230 for (e = 0; e < Ne; ++e) { 231 PetscFEGeom fegeom; 232 233 if (isAffine) { 234 fegeom.v = x; 235 fegeom.xi = cgeom->xi; 236 fegeom.J = &cgeom->J[e*dE*dE]; 237 fegeom.invJ = &cgeom->invJ[e*dE*dE]; 238 fegeom.detJ = &cgeom->detJ[e]; 239 } 240 for (q = 0; q < Nq; ++q) { 241 PetscScalar integrand; 242 PetscReal w; 243 244 if (isAffine) { 245 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x); 246 } else { 247 fegeom.v = &cgeom->v[(e*Np+q)*dE]; 248 fegeom.J = &cgeom->J[(e*Np+q)*dE*dE]; 249 fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 250 fegeom.detJ = &cgeom->detJ[e*Np+q]; 251 } 252 w = fegeom.detJ[0]*quadWeights[q]; 253 if (debug > 1 && q < Np) { 254 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 255 #if !defined(PETSC_USE_COMPLEX) 256 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 257 #endif 258 } 259 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 260 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr); 261 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 262 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); 263 integrand *= w; 264 integral[e*Nf+field] += integrand; 265 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));CHKERRQ(ierr);} 266 } 267 cOffset += totDim; 268 cOffsetAux += totDimAux; 269 } 270 PetscFunctionReturn(0); 271 } 272 273 static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, 274 PetscBdPointFunc obj_func, 275 PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[]) 276 { 277 const PetscInt debug = 0; 278 PetscFE fe; 279 PetscQuadrature quad; 280 PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal; 281 const PetscScalar *constants; 282 PetscReal *x; 283 PetscReal **B, **D, **BAux = NULL, **DAux = NULL; 284 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 285 PetscBool isAffine, auxOnBd; 286 const PetscReal *quadPoints, *quadWeights; 287 PetscInt qNc, Nq, q, Np, dE; 288 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e; 289 PetscErrorCode ierr; 290 291 PetscFunctionBegin; 292 if (!obj_func) PetscFunctionReturn(0); 293 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 294 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 295 ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr); 296 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 297 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 298 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 299 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 300 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 301 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 302 ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr); 303 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr); 304 ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr); 305 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 306 if (dsAux) { 307 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr); 308 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 309 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 310 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 311 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 312 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 313 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 314 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 315 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 316 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 317 else {ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 318 } 319 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 320 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 321 Np = fgeom->numPoints; 322 dE = fgeom->dimEmbed; 323 isAffine = fgeom->isAffine; 324 for (e = 0; e < Ne; ++e) { 325 PetscFEGeom fegeom, cgeom; 326 const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */ 327 fegeom.n = 0; 328 fegeom.v = 0; 329 fegeom.J = 0; 330 fegeom.detJ = 0; 331 if (isAffine) { 332 fegeom.v = x; 333 fegeom.xi = fgeom->xi; 334 fegeom.J = &fgeom->J[e*dE*dE]; 335 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 336 fegeom.detJ = &fgeom->detJ[e]; 337 fegeom.n = &fgeom->n[e*dE]; 338 339 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 340 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 341 cgeom.detJ = &fgeom->suppDetJ[0][e]; 342 } 343 for (q = 0; q < Nq; ++q) { 344 PetscScalar integrand; 345 PetscReal w; 346 347 if (isAffine) { 348 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 349 } else { 350 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 351 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 352 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 353 fegeom.detJ = &fgeom->detJ[e*Np+q]; 354 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 355 356 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 357 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 358 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 359 } 360 w = fegeom.detJ[0]*quadWeights[q]; 361 if (debug > 1 && q < Np) { 362 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 363 #ifndef PETSC_USE_COMPLEX 364 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 365 #endif 366 } 367 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 368 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr); 369 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 370 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); 371 integrand *= w; 372 integral[e*Nf+field] += integrand; 373 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);} 374 } 375 cOffset += totDim; 376 cOffsetAux += totDimAux; 377 } 378 PetscFunctionReturn(0); 379 } 380 381 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 382 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 383 { 384 const PetscInt debug = 0; 385 PetscFE fe; 386 PetscPointFunc f0_func; 387 PetscPointFunc f1_func; 388 PetscQuadrature quad; 389 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 390 const PetscScalar *constants; 391 PetscReal *x; 392 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 393 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 394 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 395 PetscBool isAffine; 396 const PetscReal *quadPoints, *quadWeights; 397 PetscInt qNc, Nq, q, Np, dE; 398 PetscErrorCode ierr; 399 400 PetscFunctionBegin; 401 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 402 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 403 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 404 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 405 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 406 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 407 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 408 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 409 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 410 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr); 411 ierr = PetscDSGetResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr); 412 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 413 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr); 414 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 415 if (!f0_func && !f1_func) PetscFunctionReturn(0); 416 ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr); 417 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 418 if (dsAux) { 419 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 420 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 421 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 422 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 423 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 424 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 425 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 426 ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 427 } 428 NbI = Nb[field]; 429 NcI = Nc[field]; 430 BI = B[field]; 431 DI = D[field]; 432 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 433 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 434 Np = cgeom->numPoints; 435 dE = cgeom->dimEmbed; 436 isAffine = cgeom->isAffine; 437 for (e = 0; e < Ne; ++e) { 438 PetscFEGeom fegeom; 439 440 if (isAffine) { 441 fegeom.v = x; 442 fegeom.xi = cgeom->xi; 443 fegeom.J = &cgeom->J[e*dE*dE]; 444 fegeom.invJ = &cgeom->invJ[e*dE*dE]; 445 fegeom.detJ = &cgeom->detJ[e]; 446 } 447 ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr); 448 ierr = PetscArrayzero(f1, Nq*NcI*dim);CHKERRQ(ierr); 449 for (q = 0; q < Nq; ++q) { 450 PetscReal w; 451 PetscInt c, d; 452 453 if (isAffine) { 454 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x); 455 } else { 456 fegeom.v = &cgeom->v[(e*Np+q)*dE]; 457 fegeom.J = &cgeom->J[(e*Np+q)*dE*dE]; 458 fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 459 fegeom.detJ = &cgeom->detJ[e*Np+q]; 460 } 461 w = fegeom.detJ[0]*quadWeights[q]; 462 if (debug > 1 && q < Np) { 463 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 464 #if !defined(PETSC_USE_COMPLEX) 465 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 466 #endif 467 } 468 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 469 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr); 470 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 471 if (f0_func) { 472 f0_func(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*NcI]); 473 for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w; 474 } 475 if (f1_func) { 476 f1_func(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*NcI*dim]); 477 for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w; 478 } 479 } 480 ierr = PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, BI, DI, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr); 481 cOffset += totDim; 482 cOffsetAux += totDimAux; 483 } 484 PetscFunctionReturn(0); 485 } 486 487 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 488 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 489 { 490 const PetscInt debug = 0; 491 PetscFE fe; 492 PetscBdPointFunc f0_func; 493 PetscBdPointFunc f1_func; 494 PetscQuadrature quad; 495 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 496 const PetscScalar *constants; 497 PetscReal *x; 498 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 499 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 500 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 501 PetscBool isAffine, auxOnBd = PETSC_FALSE; 502 const PetscReal *quadPoints, *quadWeights; 503 PetscInt qNc, Nq, q, Np, dE; 504 PetscErrorCode ierr; 505 506 PetscFunctionBegin; 507 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 508 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 509 ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr); 510 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 511 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 512 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 513 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 514 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 515 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 516 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr); 517 ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr); 518 if (!f0_func && !f1_func) PetscFunctionReturn(0); 519 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 520 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr); 521 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 522 ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr); 523 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 524 if (dsAux) { 525 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr); 526 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 527 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 528 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 529 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 530 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 531 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 532 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 533 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 534 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 535 else {ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 536 } 537 NbI = Nb[field]; 538 NcI = Nc[field]; 539 BI = B[field]; 540 DI = D[field]; 541 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 542 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 543 Np = fgeom->numPoints; 544 dE = fgeom->dimEmbed; 545 isAffine = fgeom->isAffine; 546 for (e = 0; e < Ne; ++e) { 547 PetscFEGeom fegeom, cgeom; 548 const PetscInt face = fgeom->face[e][0]; 549 fegeom.n = 0; 550 fegeom.v = 0; 551 fegeom.J = 0; 552 fegeom.detJ = 0; 553 if (isAffine) { 554 fegeom.v = x; 555 fegeom.xi = fgeom->xi; 556 fegeom.J = &fgeom->J[e*dE*dE]; 557 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 558 fegeom.detJ = &fgeom->detJ[e]; 559 fegeom.n = &fgeom->n[e*dE]; 560 561 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 562 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 563 cgeom.detJ = &fgeom->suppDetJ[0][e]; 564 } 565 ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr); 566 ierr = PetscArrayzero(f1, Nq*NcI*dim);CHKERRQ(ierr); 567 for (q = 0; q < Nq; ++q) { 568 PetscReal w; 569 PetscInt c, d; 570 571 if (isAffine) { 572 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 573 } else { 574 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 575 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 576 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 577 fegeom.detJ = &fgeom->detJ[e*Np+q]; 578 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 579 580 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 581 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 582 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 583 } 584 w = fegeom.detJ[0]*quadWeights[q]; 585 if (debug > 1 && q < Np) { 586 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 587 #if !defined(PETSC_USE_COMPLEX) 588 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 589 #endif 590 } 591 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 592 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr); 593 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, auxOnBd ? q : face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 594 if (f0_func) { 595 f0_func(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]); 596 for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w; 597 } 598 if (f1_func) { 599 f1_func(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]); 600 for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w; 601 } 602 } 603 ierr = PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr); 604 cOffset += totDim; 605 cOffsetAux += totDimAux; 606 } 607 PetscFunctionReturn(0); 608 } 609 610 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom, 611 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 612 { 613 const PetscInt debug = 0; 614 PetscFE feI, feJ; 615 PetscPointJac g0_func, g1_func, g2_func, g3_func; 616 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 617 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 618 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 619 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 620 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 621 PetscQuadrature quad; 622 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 623 const PetscScalar *constants; 624 PetscReal *x; 625 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 626 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 627 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 628 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 629 PetscInt dE, Np; 630 PetscBool isAffine; 631 const PetscReal *quadPoints, *quadWeights; 632 PetscInt qNc, Nq, q; 633 PetscErrorCode ierr; 634 635 PetscFunctionBegin; 636 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr); 637 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr); 638 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr); 639 ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr); 640 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 641 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 642 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 643 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 644 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 645 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 646 switch(jtype) { 647 case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 648 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 649 case PETSCFE_JACOBIAN: ierr = PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 650 } 651 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0); 652 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 653 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr); 654 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 655 ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr); 656 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr); 657 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr); 658 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 659 if (dsAux) { 660 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 661 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 662 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 663 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 664 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 665 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 666 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 667 ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 668 } 669 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 670 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 671 BI = B[fieldI], BJ = B[fieldJ]; 672 DI = D[fieldI], DJ = D[fieldJ]; 673 /* Initialize here in case the function is not defined */ 674 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 675 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 676 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 677 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 678 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 679 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 680 Np = cgeom->numPoints; 681 dE = cgeom->dimEmbed; 682 isAffine = cgeom->isAffine; 683 for (e = 0; e < Ne; ++e) { 684 PetscFEGeom fegeom; 685 686 if (isAffine) { 687 fegeom.v = x; 688 fegeom.xi = cgeom->xi; 689 fegeom.J = &cgeom->J[e*dE*dE]; 690 fegeom.invJ = &cgeom->invJ[e*dE*dE]; 691 fegeom.detJ = &cgeom->detJ[e]; 692 } 693 for (q = 0; q < Nq; ++q) { 694 const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ]; 695 const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim]; 696 PetscReal w; 697 PetscInt c; 698 699 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 700 if (isAffine) { 701 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x); 702 } else { 703 fegeom.v = &cgeom->v[(e*Np+q)*dE]; 704 fegeom.J = &cgeom->J[(e*Np+q)*dE*dE]; 705 fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 706 fegeom.detJ = &cgeom->detJ[e*Np+q]; 707 } 708 w = fegeom.detJ[0]*quadWeights[q]; 709 if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);} 710 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 711 if (g0_func) { 712 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 713 g0_func(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); 714 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 715 } 716 if (g1_func) { 717 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 718 g1_func(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); 719 for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w; 720 } 721 if (g2_func) { 722 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 723 g2_func(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); 724 for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w; 725 } 726 if (g3_func) { 727 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 728 g3_func(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); 729 for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w; 730 } 731 732 ierr = PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr); 733 } 734 if (debug > 1) { 735 PetscInt fc, f, gc, g; 736 737 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 738 for (fc = 0; fc < NcI; ++fc) { 739 for (f = 0; f < NbI; ++f) { 740 const PetscInt i = offsetI + f*NcI+fc; 741 for (gc = 0; gc < NcJ; ++gc) { 742 for (g = 0; g < NbJ; ++g) { 743 const PetscInt j = offsetJ + g*NcJ+gc; 744 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 745 } 746 } 747 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 748 } 749 } 750 } 751 cOffset += totDim; 752 cOffsetAux += totDimAux; 753 eOffset += PetscSqr(totDim); 754 } 755 PetscFunctionReturn(0); 756 } 757 758 static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 759 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 760 { 761 const PetscInt debug = 0; 762 PetscFE feI, feJ; 763 PetscBdPointJac g0_func, g1_func, g2_func, g3_func; 764 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 765 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 766 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 767 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 768 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 769 PetscQuadrature quad; 770 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 771 const PetscScalar *constants; 772 PetscReal *x; 773 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 774 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 775 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 776 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 777 PetscBool isAffine; 778 const PetscReal *quadPoints, *quadWeights; 779 PetscInt qNc, Nq, q, Np, dE; 780 PetscErrorCode ierr; 781 782 PetscFunctionBegin; 783 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr); 784 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr); 785 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr); 786 ierr = PetscFEGetFaceQuadrature(feI, &quad);CHKERRQ(ierr); 787 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 788 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 789 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 790 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 791 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 792 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 793 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr); 794 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr); 795 ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr); 796 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0); 797 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 798 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr); 799 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 800 ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr); 801 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 802 if (dsAux) { 803 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 804 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 805 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 806 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 807 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 808 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 809 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 810 ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 811 } 812 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 813 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 814 BI = B[fieldI], BJ = B[fieldJ]; 815 DI = D[fieldI], DJ = D[fieldJ]; 816 /* Initialize here in case the function is not defined */ 817 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 818 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 819 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 820 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 821 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 822 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 823 Np = fgeom->numPoints; 824 dE = fgeom->dimEmbed; 825 isAffine = fgeom->isAffine; 826 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 827 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 828 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 829 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 830 for (e = 0; e < Ne; ++e) { 831 PetscFEGeom fegeom, cgeom; 832 const PetscInt face = fgeom->face[e][0]; 833 fegeom.n = 0; 834 fegeom.v = 0; 835 fegeom.J = 0; 836 fegeom.detJ = 0; 837 if (isAffine) { 838 fegeom.v = x; 839 fegeom.xi = fgeom->xi; 840 fegeom.J = &fgeom->J[e*dE*dE]; 841 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 842 fegeom.detJ = &fgeom->detJ[e]; 843 fegeom.n = &fgeom->n[e*dE]; 844 845 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 846 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 847 cgeom.detJ = &fgeom->suppDetJ[0][e]; 848 } 849 for (q = 0; q < Nq; ++q) { 850 const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ]; 851 const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim]; 852 PetscReal w; 853 PetscInt c; 854 855 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 856 if (isAffine) { 857 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 858 } else { 859 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 860 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 861 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 862 fegeom.detJ = &fgeom->detJ[e*Np+q]; 863 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 864 865 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 866 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 867 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 868 } 869 w = fegeom.detJ[0]*quadWeights[q]; 870 if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);} 871 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 872 if (g0_func) { 873 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 874 g0_func(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); 875 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 876 } 877 if (g1_func) { 878 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 879 g1_func(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); 880 for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w; 881 } 882 if (g2_func) { 883 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 884 g2_func(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); 885 for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w; 886 } 887 if (g3_func) { 888 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 889 g3_func(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); 890 for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w; 891 } 892 893 ierr = PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr); 894 } 895 if (debug > 1) { 896 PetscInt fc, f, gc, g; 897 898 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 899 for (fc = 0; fc < NcI; ++fc) { 900 for (f = 0; f < NbI; ++f) { 901 const PetscInt i = offsetI + f*NcI+fc; 902 for (gc = 0; gc < NcJ; ++gc) { 903 for (g = 0; g < NbJ; ++g) { 904 const PetscInt j = offsetJ + g*NcJ+gc; 905 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 906 } 907 } 908 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 909 } 910 } 911 } 912 cOffset += totDim; 913 cOffsetAux += totDimAux; 914 eOffset += PetscSqr(totDim); 915 } 916 PetscFunctionReturn(0); 917 } 918 919 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 920 { 921 PetscFunctionBegin; 922 fem->ops->setfromoptions = NULL; 923 fem->ops->setup = PetscFESetUp_Basic; 924 fem->ops->view = PetscFEView_Basic; 925 fem->ops->destroy = PetscFEDestroy_Basic; 926 fem->ops->getdimension = PetscFEGetDimension_Basic; 927 fem->ops->gettabulation = PetscFEGetTabulation_Basic; 928 fem->ops->integrate = PetscFEIntegrate_Basic; 929 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 930 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 931 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 932 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 933 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 934 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 935 PetscFunctionReturn(0); 936 } 937 938 /*MC 939 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization 940 941 Level: intermediate 942 943 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 944 M*/ 945 946 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 947 { 948 PetscFE_Basic *b; 949 PetscErrorCode ierr; 950 951 PetscFunctionBegin; 952 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 953 ierr = PetscNewLog(fem,&b);CHKERRQ(ierr); 954 fem->data = b; 955 956 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr); 957 PetscFunctionReturn(0); 958 } 959