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