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 fegeom.n = 0; 325 fegeom.v = 0; 326 fegeom.J = 0; 327 fegeom.detJ = 0; 328 if (isAffine) { 329 fegeom.v = x; 330 fegeom.xi = fgeom->xi; 331 fegeom.J = &fgeom->J[e*dE*dE]; 332 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 333 fegeom.detJ = &fgeom->detJ[e]; 334 fegeom.n = &fgeom->n[e*dE]; 335 336 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 337 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 338 cgeom.detJ = &fgeom->suppDetJ[0][e]; 339 } 340 for (q = 0; q < Nq; ++q) { 341 PetscScalar integrand; 342 PetscReal w; 343 344 if (isAffine) { 345 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 346 } else { 347 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 348 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 349 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 350 fegeom.detJ = &fgeom->detJ[e*Np+q]; 351 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 352 353 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 354 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 355 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 356 } 357 w = fegeom.detJ[0]*quadWeights[q]; 358 if (debug > 1 && q < Np) { 359 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 360 #ifndef PETSC_USE_COMPLEX 361 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 362 #endif 363 } 364 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 365 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr); 366 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);} 367 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); 368 integrand *= w; 369 integral[e*Nf+field] += integrand; 370 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);} 371 } 372 cOffset += totDim; 373 cOffsetAux += totDimAux; 374 } 375 PetscFunctionReturn(0); 376 } 377 378 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, 379 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 380 { 381 const PetscInt debug = 0; 382 PetscFE fe; 383 PetscPointFunc f0_func; 384 PetscPointFunc f1_func; 385 PetscQuadrature quad; 386 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 387 const PetscScalar *constants; 388 PetscReal *x; 389 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 390 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 391 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 392 PetscBool isAffine; 393 const PetscReal *quadPoints, *quadWeights; 394 PetscInt qNc, Nq, q, Np, dE; 395 PetscErrorCode ierr; 396 397 PetscFunctionBegin; 398 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 399 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 400 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 401 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 402 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 403 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 404 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 405 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 406 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 407 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr); 408 ierr = PetscDSGetResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr); 409 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 410 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr); 411 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 412 if (!f0_func && !f1_func) PetscFunctionReturn(0); 413 ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr); 414 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 415 if (dsAux) { 416 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 417 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 418 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 419 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 420 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 421 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 422 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 423 ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 424 } 425 NbI = Nb[field]; 426 NcI = Nc[field]; 427 BI = B[field]; 428 DI = D[field]; 429 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 430 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 431 Np = cgeom->numPoints; 432 dE = cgeom->dimEmbed; 433 isAffine = cgeom->isAffine; 434 for (e = 0; e < Ne; ++e) { 435 PetscFEGeom fegeom; 436 437 if (isAffine) { 438 fegeom.v = x; 439 fegeom.xi = cgeom->xi; 440 fegeom.J = &cgeom->J[e*dE*dE]; 441 fegeom.invJ = &cgeom->invJ[e*dE*dE]; 442 fegeom.detJ = &cgeom->detJ[e]; 443 } 444 ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr); 445 ierr = PetscArrayzero(f1, Nq*NcI*dim);CHKERRQ(ierr); 446 for (q = 0; q < Nq; ++q) { 447 PetscReal w; 448 PetscInt c, d; 449 450 if (isAffine) { 451 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x); 452 } else { 453 fegeom.v = &cgeom->v[(e*Np+q)*dE]; 454 fegeom.J = &cgeom->J[(e*Np+q)*dE*dE]; 455 fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 456 fegeom.detJ = &cgeom->detJ[e*Np+q]; 457 } 458 w = fegeom.detJ[0]*quadWeights[q]; 459 if (debug > 1 && q < Np) { 460 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 461 #if !defined(PETSC_USE_COMPLEX) 462 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 463 #endif 464 } 465 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 466 ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr); 467 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 468 if (f0_func) { 469 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]); 470 for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w; 471 } 472 if (f1_func) { 473 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]); 474 for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w; 475 } 476 } 477 ierr = PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, BI, DI, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr); 478 cOffset += totDim; 479 cOffsetAux += totDimAux; 480 } 481 PetscFunctionReturn(0); 482 } 483 484 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom, 485 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[]) 486 { 487 const PetscInt debug = 0; 488 PetscFE fe; 489 PetscBdPointFunc f0_func; 490 PetscBdPointFunc f1_func; 491 PetscQuadrature quad; 492 PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal; 493 const PetscScalar *constants; 494 PetscReal *x; 495 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI; 496 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 497 PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI; 498 PetscBool isAffine, auxOnBd = PETSC_FALSE; 499 const PetscReal *quadPoints, *quadWeights; 500 PetscInt qNc, Nq, q, Np, dE; 501 PetscErrorCode ierr; 502 503 PetscFunctionBegin; 504 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr); 505 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr); 506 ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr); 507 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 508 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 509 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 510 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 511 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 512 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 513 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr); 514 ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr); 515 if (!f0_func && !f1_func) PetscFunctionReturn(0); 516 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 517 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr); 518 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 519 ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr); 520 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 521 if (dsAux) { 522 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr); 523 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 524 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 525 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 526 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 527 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 528 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 529 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 530 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE; 531 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 532 else {ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);} 533 } 534 NbI = Nb[field]; 535 NcI = Nc[field]; 536 BI = B[field]; 537 DI = D[field]; 538 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 539 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 540 Np = fgeom->numPoints; 541 dE = fgeom->dimEmbed; 542 isAffine = fgeom->isAffine; 543 for (e = 0; e < Ne; ++e) { 544 PetscFEGeom fegeom, cgeom; 545 const PetscInt face = fgeom->face[e][0]; 546 fegeom.n = 0; 547 fegeom.v = 0; 548 fegeom.J = 0; 549 fegeom.detJ = 0; 550 if (isAffine) { 551 fegeom.v = x; 552 fegeom.xi = fgeom->xi; 553 fegeom.J = &fgeom->J[e*dE*dE]; 554 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 555 fegeom.detJ = &fgeom->detJ[e]; 556 fegeom.n = &fgeom->n[e*dE]; 557 558 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 559 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 560 cgeom.detJ = &fgeom->suppDetJ[0][e]; 561 } 562 ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr); 563 ierr = PetscArrayzero(f1, Nq*NcI*dim);CHKERRQ(ierr); 564 for (q = 0; q < Nq; ++q) { 565 PetscReal w; 566 PetscInt c, d; 567 568 if (isAffine) { 569 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 570 } else { 571 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 572 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 573 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 574 fegeom.detJ = &fgeom->detJ[e*Np+q]; 575 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 576 577 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 578 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 579 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 580 } 581 w = fegeom.detJ[0]*quadWeights[q]; 582 if (debug > 1 && q < Np) { 583 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr); 584 #if !defined(PETSC_USE_COMPLEX) 585 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr); 586 #endif 587 } 588 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 589 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); 590 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);} 591 if (f0_func) { 592 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]); 593 for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w; 594 } 595 if (f1_func) { 596 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]); 597 for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w; 598 } 599 } 600 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); 601 cOffset += totDim; 602 cOffsetAux += totDimAux; 603 } 604 PetscFunctionReturn(0); 605 } 606 607 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom, 608 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 609 { 610 const PetscInt debug = 0; 611 PetscFE feI, feJ; 612 PetscPointJac g0_func, g1_func, g2_func, g3_func; 613 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 614 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 615 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 616 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 617 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 618 PetscQuadrature quad; 619 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 620 const PetscScalar *constants; 621 PetscReal *x; 622 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 623 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 624 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 625 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 626 PetscInt dE, Np; 627 PetscBool isAffine; 628 const PetscReal *quadPoints, *quadWeights; 629 PetscInt qNc, Nq, q; 630 PetscErrorCode ierr; 631 632 PetscFunctionBegin; 633 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr); 634 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr); 635 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr); 636 ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr); 637 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 638 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 639 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 640 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 641 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 642 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 643 switch(jtype) { 644 case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 645 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 646 case PETSCFE_JACOBIAN: ierr = PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break; 647 } 648 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0); 649 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 650 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr); 651 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 652 ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr); 653 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr); 654 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr); 655 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 656 if (dsAux) { 657 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 658 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 659 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 660 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 661 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 662 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 663 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 664 ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 665 } 666 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 667 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 668 BI = B[fieldI], BJ = B[fieldJ]; 669 DI = D[fieldI], DJ = D[fieldJ]; 670 /* Initialize here in case the function is not defined */ 671 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 672 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 673 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 674 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 675 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 676 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 677 Np = cgeom->numPoints; 678 dE = cgeom->dimEmbed; 679 isAffine = cgeom->isAffine; 680 for (e = 0; e < Ne; ++e) { 681 PetscFEGeom fegeom; 682 683 if (isAffine) { 684 fegeom.v = x; 685 fegeom.xi = cgeom->xi; 686 fegeom.J = &cgeom->J[e*dE*dE]; 687 fegeom.invJ = &cgeom->invJ[e*dE*dE]; 688 fegeom.detJ = &cgeom->detJ[e]; 689 } 690 for (q = 0; q < Nq; ++q) { 691 const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ]; 692 const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim]; 693 PetscReal w; 694 PetscInt c; 695 696 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 697 if (isAffine) { 698 CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x); 699 } else { 700 fegeom.v = &cgeom->v[(e*Np+q)*dE]; 701 fegeom.J = &cgeom->J[(e*Np+q)*dE*dE]; 702 fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE]; 703 fegeom.detJ = &cgeom->detJ[e*Np+q]; 704 } 705 w = fegeom.detJ[0]*quadWeights[q]; 706 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);} 707 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);} 708 if (g0_func) { 709 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 710 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); 711 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 712 } 713 if (g1_func) { 714 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 715 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); 716 for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w; 717 } 718 if (g2_func) { 719 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 720 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); 721 for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w; 722 } 723 if (g3_func) { 724 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 725 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); 726 for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w; 727 } 728 729 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); 730 } 731 if (debug > 1) { 732 PetscInt fc, f, gc, g; 733 734 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 735 for (fc = 0; fc < NcI; ++fc) { 736 for (f = 0; f < NbI; ++f) { 737 const PetscInt i = offsetI + f*NcI+fc; 738 for (gc = 0; gc < NcJ; ++gc) { 739 for (g = 0; g < NbJ; ++g) { 740 const PetscInt j = offsetJ + g*NcJ+gc; 741 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 742 } 743 } 744 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 745 } 746 } 747 } 748 cOffset += totDim; 749 cOffsetAux += totDimAux; 750 eOffset += PetscSqr(totDim); 751 } 752 PetscFunctionReturn(0); 753 } 754 755 PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom, 756 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[]) 757 { 758 const PetscInt debug = 0; 759 PetscFE feI, feJ; 760 PetscBdPointJac g0_func, g1_func, g2_func, g3_func; 761 PetscInt cOffset = 0; /* Offset into coefficients[] for element e */ 762 PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */ 763 PetscInt eOffset = 0; /* Offset into elemMat[] for element e */ 764 PetscInt offsetI = 0; /* Offset into an element vector for fieldI */ 765 PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */ 766 PetscQuadrature quad; 767 PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal; 768 const PetscScalar *constants; 769 PetscReal *x; 770 PetscReal **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ; 771 PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL; 772 PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0; 773 PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e; 774 PetscBool isAffine; 775 const PetscReal *quadPoints, *quadWeights; 776 PetscInt qNc, Nq, q, Np, dE; 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr); 781 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr); 782 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr); 783 ierr = PetscFEGetFaceQuadrature(feI, &quad);CHKERRQ(ierr); 784 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 785 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 786 ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr); 787 ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr); 788 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr); 789 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr); 790 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr); 791 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr); 792 ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr); 793 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr); 794 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr); 795 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr); 796 ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr); 797 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr); 798 if (dsAux) { 799 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr); 800 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr); 801 ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr); 802 ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr); 803 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr); 804 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr); 805 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr); 806 ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr); 807 } 808 NbI = Nb[fieldI], NbJ = Nb[fieldJ]; 809 NcI = Nc[fieldI], NcJ = Nc[fieldJ]; 810 BI = B[fieldI], BJ = B[fieldJ]; 811 DI = D[fieldI], DJ = D[fieldJ]; 812 /* Initialize here in case the function is not defined */ 813 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 814 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 815 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 816 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 817 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 818 if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc); 819 Np = fgeom->numPoints; 820 dE = fgeom->dimEmbed; 821 isAffine = fgeom->isAffine; 822 for (e = 0; e < Ne; ++e) { 823 PetscFEGeom fegeom, cgeom; 824 const PetscInt face = fgeom->face[e][0]; 825 fegeom.n = 0; 826 fegeom.v = 0; 827 fegeom.J = 0; 828 fegeom.detJ = 0; 829 if (isAffine) { 830 fegeom.v = x; 831 fegeom.xi = fgeom->xi; 832 fegeom.J = &fgeom->J[e*dE*dE]; 833 fegeom.invJ = &fgeom->invJ[e*dE*dE]; 834 fegeom.detJ = &fgeom->detJ[e]; 835 fegeom.n = &fgeom->n[e*dE]; 836 837 cgeom.J = &fgeom->suppJ[0][e*dE*dE]; 838 cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE]; 839 cgeom.detJ = &fgeom->suppDetJ[0][e]; 840 } 841 for (q = 0; q < Nq; ++q) { 842 const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ]; 843 const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim]; 844 PetscReal w; 845 PetscInt c; 846 847 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);} 848 if (isAffine) { 849 CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x); 850 } else { 851 fegeom.v = &fgeom->v[(e*Np+q)*dE]; 852 fegeom.J = &fgeom->J[(e*Np+q)*dE*dE]; 853 fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE]; 854 fegeom.detJ = &fgeom->detJ[e*Np+q]; 855 fegeom.n = &fgeom->n[(e*Np+q)*dE]; 856 857 cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE]; 858 cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE]; 859 cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q]; 860 } 861 w = fegeom.detJ[0]*quadWeights[q]; 862 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);} 863 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);} 864 if (g0_func) { 865 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr); 866 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); 867 for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w; 868 } 869 if (g1_func) { 870 ierr = PetscArrayzero(g1, NcI*NcJ*dim);CHKERRQ(ierr); 871 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); 872 for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w; 873 } 874 if (g2_func) { 875 ierr = PetscArrayzero(g2, NcI*NcJ*dim);CHKERRQ(ierr); 876 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); 877 for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w; 878 } 879 if (g3_func) { 880 ierr = PetscArrayzero(g3, NcI*NcJ*dim*dim);CHKERRQ(ierr); 881 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); 882 for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w; 883 } 884 885 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); 886 } 887 if (debug > 1) { 888 PetscInt fc, f, gc, g; 889 890 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr); 891 for (fc = 0; fc < NcI; ++fc) { 892 for (f = 0; f < NbI; ++f) { 893 const PetscInt i = offsetI + f*NcI+fc; 894 for (gc = 0; gc < NcJ; ++gc) { 895 for (g = 0; g < NbJ; ++g) { 896 const PetscInt j = offsetJ + g*NcJ+gc; 897 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr); 898 } 899 } 900 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr); 901 } 902 } 903 } 904 cOffset += totDim; 905 cOffsetAux += totDimAux; 906 eOffset += PetscSqr(totDim); 907 } 908 PetscFunctionReturn(0); 909 } 910 911 PetscErrorCode PetscFEInitialize_Basic(PetscFE fem) 912 { 913 PetscFunctionBegin; 914 fem->ops->setfromoptions = NULL; 915 fem->ops->setup = PetscFESetUp_Basic; 916 fem->ops->view = PetscFEView_Basic; 917 fem->ops->destroy = PetscFEDestroy_Basic; 918 fem->ops->getdimension = PetscFEGetDimension_Basic; 919 fem->ops->gettabulation = PetscFEGetTabulation_Basic; 920 fem->ops->integrate = PetscFEIntegrate_Basic; 921 fem->ops->integratebd = PetscFEIntegrateBd_Basic; 922 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 923 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 924 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 925 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 926 fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic; 927 PetscFunctionReturn(0); 928 } 929 930 /*MC 931 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization 932 933 Level: intermediate 934 935 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 936 M*/ 937 938 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem) 939 { 940 PetscFE_Basic *b; 941 PetscErrorCode ierr; 942 943 PetscFunctionBegin; 944 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 945 ierr = PetscNewLog(fem,&b);CHKERRQ(ierr); 946 fem->data = b; 947 948 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr); 949 PetscFunctionReturn(0); 950 } 951