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