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