Skip to content

Commit

Permalink
Merge branch 'knepley/feature-plex-bdfunc'
Browse files Browse the repository at this point in the history 
* knepley/feature-plex-bdfunc: (56 commits)
  updated snes/ex56
  updated snes/ex56
  Plex: Goddamn it! Shadowing a variable name
  PetscDS: Must be setup in order to access field sizes and components
  fix bad merge 44ac993 wrt src/snes/examples/tutorials/ex56.c [by attempting to redo the merge with kdiff3 - and then comparing the final result with the result in the above merge]
  fixed snes/ex56 for PetscDSAddBoundary change
  dmplexsnes: Remove more unused variables
  Remove unused variables
  SNES ex77: Must also set face quadrature
  SNES ex12: No longer need boundary discretization orders
  SNES ex62: Numerical flutter
  SNES ex12: Numerical flutter
  Plex: In projection, get Nc from DS instead of calculating it again
  PetscFE: Added missing destroy
  SNEX ex77: Update boundary args and FE creation
  Plex: Update boundary integrals to new API
  SNES ex12: Tests of bc from fields
  SNES ex12: Added the ability to use finite element fields as boundary conditions
  Plex: Need invJ for field projection
  Plex: Fix bug with field projection - We must tabulate the FE on the dual basis quadrature points
  ...
  • Loading branch information
@knepley
knepley committed Oct 26, 2016
2 parents a4befa6 + 35da518 commit 66b078f
Show file tree
Hide file tree
Showing 92 changed files with 14,453 additions and 10,141 deletions.
36 changes: 28 additions & 8 deletions config/builder.py
View file

Large diffs are not rendered by default.

22 changes: 0 additions & 22 deletions include/petsc/finclude/ftn-custom/petscdt.h90
View file
Expand Up @@ -59,25 +59,3 @@
PETSCDS_HIDE prob
End Subroutine
End Interface

Interface
Subroutine PetscDSGetBdTabulation(prob,f,b,bDer,ierr)
USE_PETSCDS_HIDE
PetscInt f
PetscReal, pointer :: b(:)
PetscReal, pointer :: bDer(:)
PetscErrorCode ierr
PETSCDS_HIDE prob
End Subroutine
End Interface

Interface
Subroutine PetscDSRestoreBdTabulation(prob,f,b,bDer,ierr)
USE_PETSCDS_HIDE
PetscInt f
PetscReal, pointer :: b(:)
PetscReal, pointer :: bDer(:)
PetscErrorCode ierr
PETSCDS_HIDE prob
End Subroutine
End Interface
3 changes: 2 additions & 1 deletion include/petsc/private/dmimpl.h
View file
Expand Up @@ -63,7 +63,8 @@ struct _DMOps {

PetscErrorCode (*projectfunctionlocal)(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
PetscErrorCode (*projectfunctionlabellocal)(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
PetscErrorCode (*projectfieldlocal)(DM,Vec,void(**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec);
PetscErrorCode (*projectfieldlocal)(DM,PetscReal,Vec,void(**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec);
PetscErrorCode (*projectfieldlabellocal)(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],Vec,void(**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscScalar[]),InsertMode,Vec);
PetscErrorCode (*computel2diff)(DM,PetscReal,PetscErrorCode(**)(PetscInt, PetscReal,const PetscReal [], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
PetscErrorCode (*computel2gradientdiff)(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal [],const PetscReal[],PetscInt, PetscScalar *,void *),void **,Vec,const PetscReal[],PetscReal *);
PetscErrorCode (*computel2fielddiff)(DM,PetscReal,PetscErrorCode(**)(PetscInt, PetscReal,const PetscReal [], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
Expand Down
67 changes: 35 additions & 32 deletions include/petsc/private/petscdsimpl.h
View file
Expand Up @@ -12,7 +12,7 @@ typedef struct _n_DSBoundary *DSBoundary;
struct _n_DSBoundary {
const char *name;
const char *labelname;
PetscBool essential;
DMBoundaryConditionType type;
PetscInt field;
PetscInt numcomps;
PetscInt *comps;
Expand All @@ -33,39 +33,42 @@ struct _PetscDSOps {

struct _p_PetscDS {
PETSCHEADER(struct _PetscDSOps);
void *data; /* Implementation object */
PetscBool setup; /* Flag for setup */
PetscInt Nf; /* The number of solution fields */
PetscBool *implicit; /* Flag for implicit or explicit solve */
PetscBool *adjacency; /* Flag for variable influence */
PetscObject *disc; /* The discretization for each solution field (PetscFE, PetscFV, etc.) */
PetscObject *discBd; /* The boundary discretization for each solution field (PetscFE, PetscFV, etc.) */
PetscPointFunc *obj; /* Scalar integral (like an objective function) */
PetscPointFunc *f; /* Weak form integrands for F, f_0, f_1 */
PetscPointJac *g; /* Weak form integrands for J = dF/du, g_0, g_1, g_2, g_3 */
PetscPointJac *gp; /* Weak form integrands for preconditioner for J, g_0, g_1, g_2, g_3 */
PetscPointJac *gt; /* Weak form integrands for dF/du_t, g_0, g_1, g_2, g_3 */
PetscBdPointFunc *fBd; /* Weak form boundary integrands F_bd, f_0, f_1 */
PetscBdPointJac *gBd; /* Weak form boundary integrands J_bd = dF_bd/du, g_0, g_1, g_2, g_3 */
PetscRiemannFunc *r; /* Riemann solvers */
void **ctx; /* User contexts for each field */
PetscInt dim; /* The spatial dimension */
void *data; /* Implementation object */
PetscBool setup; /* Flag for setup */
PetscInt Nf; /* The number of solution fields */
PetscBool *implicit; /* Flag for implicit or explicit solve */
PetscBool *adjacency; /* Flag for variable influence */
PetscObject *disc; /* The discretization for each solution field (PetscFE, PetscFV, etc.) */
PetscPointFunc *obj; /* Scalar integral (like an objective function) */
PetscPointFunc *f; /* Weak form integrands for F, f_0, f_1 */
PetscPointJac *g; /* Weak form integrands for J = dF/du, g_0, g_1, g_2, g_3 */
PetscPointJac *gp; /* Weak form integrands for preconditioner for J, g_0, g_1, g_2, g_3 */
PetscPointJac *gt; /* Weak form integrands for dF/du_t, g_0, g_1, g_2, g_3 */
PetscBdPointFunc *fBd; /* Weak form boundary integrands F_bd, f_0, f_1 */
PetscBdPointJac *gBd; /* Weak form boundary integrands J_bd = dF_bd/du, g_0, g_1, g_2, g_3 */
PetscRiemannFunc *r; /* Riemann solvers */
void **ctx; /* User contexts for each field */
PetscInt dim; /* The spatial dimension */
/* Computed sizes */
PetscInt totDim, totDimBd; /* Total system dimension */
PetscInt totComp; /* Total field components */
PetscInt totDim; /* Total system dimension */
PetscInt totComp; /* Total field components */
PetscInt *Nc; /* Number of components for each field */
PetscInt *Nb; /* Number of basis functions for each field */
PetscInt *off; /* Offsets for each field */
PetscInt *offDer; /* Derivative offsets for each field */
PetscReal **basis; /* Default basis tabulation for each field */
PetscReal **basisDer; /* Default basis derivative tabulation for each field */
PetscReal **basisFace; /* Basis tabulation for each local face and field */
PetscReal **basisDerFace; /* Basis derivative tabulation for each local face and field */
/* Work space */
PetscInt *off, *offBd; /* Offsets for each field */
PetscInt *offDer, *offDerBd; /* Derivative offsets for each field */
PetscReal **basis, **basisBd; /* Default basis tabulation for each field */
PetscReal **basisDer, **basisDerBd; /* Default basis derivative tabulation for each field */
PetscScalar *u; /* Field evaluation */
PetscScalar *u_t; /* Field time derivative evaluation */
PetscScalar *u_x; /* Field gradient evaluation */
PetscScalar *refSpaceDer; /* Workspace for computing derivative in the reference coordinates */
PetscReal *x; /* Workspace for computing real coordinates */
PetscScalar *f0, *f1; /* Point evaluations of weak form residual integrands */
PetscScalar *g0, *g1, *g2, *g3; /* Point evaluations of weak form Jacobian integrands */
DSBoundary boundary; /* Linked list of boundary conditions */
PetscScalar *u; /* Field evaluation */
PetscScalar *u_t; /* Field time derivative evaluation */
PetscScalar *u_x; /* Field gradient evaluation */
PetscScalar *refSpaceDer; /* Workspace for computing derivative in the reference coordinates */
PetscReal *x; /* Workspace for computing real coordinates */
PetscScalar *f0, *f1; /* Point evaluations of weak form residual integrands */
PetscScalar *g0, *g1, *g2, *g3; /* Point evaluations of weak form Jacobian integrands */
DSBoundary boundary; /* Linked list of boundary conditions */
};

typedef struct {
Expand Down
123 changes: 48 additions & 75 deletions include/petsc/private/petscfeimpl.h
View file
Expand Up @@ -94,23 +94,25 @@ struct _PetscFEOps {
/* Element integration */
PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscReal[]);
PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEFaceGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFEFaceGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
};

struct _p_PetscFE {
PETSCHEADER(struct _PetscFEOps);
void *data; /* Implementation object */
PetscSpace basisSpace; /* The basis space P */
PetscDualSpace dualSpace; /* The dual space P' */
PetscInt numComponents; /* The number of field components */
PetscQuadrature quadrature; /* Suitable quadrature on K */
PetscInt *numDof; /* The number of dof on mesh points of each depth */
PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
PetscReal *F; /* Tabulation of basis at face centroids */
void *data; /* Implementation object */
PetscSpace basisSpace; /* The basis space P */
PetscDualSpace dualSpace; /* The dual space P' */
PetscInt numComponents; /* The number of field components */
PetscQuadrature quadrature; /* Suitable quadrature on K */
PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
PetscInt *numDof; /* The number of dof on mesh points of each depth */
PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
PetscReal *Bf, *Df, *Hf; /* Tabulation of basis and derivatives at quadrature points on each face */
PetscReal *F; /* Tabulation of basis at face centroids */
PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
};
Expand Down Expand Up @@ -181,83 +183,51 @@ PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef,

#undef __FUNCT__
#define __FUNCT__ "EvaluateFieldJets"
PETSC_STATIC_INLINE PetscErrorCode EvaluateFieldJets(PetscDS prob, PetscBool bd, PetscInt q, const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
PETSC_STATIC_INLINE void EvaluateFieldJets(PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscScalar refSpaceDer[], const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
{
PetscScalar *refSpaceDer;
PetscReal **basisField, **basisFieldDer;
PetscInt dOffset = 0, fOffset = 0;
PetscInt Nf, Nc, dimRef, dimReal, d, f;
PetscErrorCode ierr;
PetscInt dOffset = 0, fOffset = 0, f;

if (!prob) return 0;
ierr = PetscDSGetSpatialDimension(prob, &dimReal);CHKERRQ(ierr);
dimRef = dimReal;
if (bd) dimRef -= 1;
ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
ierr = PetscDSGetTotalComponents(prob, &Nc);CHKERRQ(ierr);
ierr = PetscDSGetRefCoordArrays(prob, NULL, &refSpaceDer);CHKERRQ(ierr);
if (bd) {ierr = PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);CHKERRQ(ierr);}
else {ierr = PetscDSGetTabulation(prob, &basisField, &basisFieldDer);CHKERRQ(ierr);}
for (d = 0; d < Nc; ++d) {u[d] = 0.0;}
for (d = 0; d < dimReal*Nc; ++d) {u_x[d] = 0.0;}
if (u_t) for (d = 0; d < Nc; ++d) {u_t[d] = 0.0;}
for (f = 0; f < Nf; ++f) {
const PetscReal *basis = basisField[f];
const PetscReal *basisDer = basisFieldDer[f];
PetscObject obj;
PetscClassId id;
PetscInt Nb, Ncf, b, c, e;

if (bd) {ierr = PetscDSGetBdDiscretization(prob, f, &obj);CHKERRQ(ierr);}
else {ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);}
ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
if (id == PETSCFE_CLASSID) {
PetscFE fe = (PetscFE) obj;

ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
ierr = PetscFEGetNumComponents(fe, &Ncf);CHKERRQ(ierr);
} else if (id == PETSCFV_CLASSID) {
PetscFV fv = (PetscFV) obj;

/* TODO Should also support reconstruction here */
Nb = 1;
ierr = PetscFVGetNumComponents(fv, &Ncf);CHKERRQ(ierr);
} else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
for (d = 0; d < dimRef*Ncf; ++d) refSpaceDer[d] = 0.0;
for (b = 0; b < Nb; ++b) {
const PetscInt Nbf = Nb[f], Ncf = Nc[f];
const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
PetscInt b, c, d, e;

for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
for (b = 0; b < Nbf; ++b) {
for (c = 0; c < Ncf; ++c) {
const PetscInt cidx = b*Ncf+c;
const PetscInt cidx = b*Ncf+c;

u[fOffset+c] += coefficients[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
for (d = 0; d < dimRef; ++d) refSpaceDer[c*dimRef+d] += coefficients[dOffset+cidx]*basisDer[(q*Nb*Ncf+cidx)*dimRef+d];
u[fOffset+c] += Bq[cidx]*coefficients[dOffset+cidx];
for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+cidx];
}
}
for (c = 0; c < Ncf; ++c) for (d = 0; d < dimReal; ++d) for (e = 0; e < dimRef; ++e) u_x[(fOffset+c)*dimReal+d] += invJ[e*dimReal+d]*refSpaceDer[c*dimRef+e];
for (c = 0; c < Ncf; ++c) for (d = 0; d < dim; ++d) for (e = 0, u_x[(fOffset+c)*dim+d] = 0.0; e < dim; ++e) u_x[(fOffset+c)*dim+d] += invJ[e*dim+d]*refSpaceDer[c*dim+e];
if (u_t) {
for (b = 0; b < Nb; ++b) {
for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
for (b = 0; b < Nbf; ++b) {
for (c = 0; c < Ncf; ++c) {
const PetscInt cidx = b*Ncf+c;

u_t[fOffset+c] += coefficients_t[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+cidx];
}
}
}
#if 0
for (c = 0; c < Ncf; ++c) {
ierr = PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));CHKERRQ(ierr);
if (u_t) {ierr = PetscPrintf(PETSC_COMM_SELF, " u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));CHKERRQ(ierr);}
for (d = 0; d < dimRef; ++d) {
ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dimReal+d]));CHKERRQ(ierr);
for (d = 0; d < dim; ++d) {
ierr = PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dim+d]));CHKERRQ(ierr);
}
}
#endif
fOffset += Ncf;
dOffset += Nb*Ncf;
dOffset += Nbf*Ncf;
}
return 0;
}


#undef __FUNCT__
#define __FUNCT__ "EvaluateFaceFields"
PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
Expand All @@ -271,7 +241,7 @@ PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt fie
ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);CHKERRQ(ierr);
ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
ierr = PetscFEGetFaceTabulation(fe, &faceBasis);CHKERRQ(ierr);
ierr = PetscFEGetFaceCentroidTabulation(fe, &faceBasis);CHKERRQ(ierr);
for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
for (b = 0; b < Nb; ++b) {
for (c = 0; c < Nc; ++c) {
Expand All @@ -285,24 +255,27 @@ PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt fie

#undef __FUNCT__
#define __FUNCT__ "TransformF"
PETSC_STATIC_INLINE void TransformF(PetscInt dimReal, PetscInt dimRef, PetscInt Nc, PetscInt q, const PetscReal invJ[], PetscReal detJ, const PetscReal quadWeights[], PetscScalar refSpaceDer[], PetscScalar f0[], PetscScalar f1[])
PETSC_STATIC_INLINE void TransformF(PetscInt dim, PetscInt Nc, PetscInt q, const PetscReal invJ[], PetscReal detJ, const PetscReal quadWeights[], PetscScalar refSpaceDer[], PetscScalar f0[], PetscScalar f1[])
{
PetscInt c, d, e;

for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= detJ*quadWeights[q];
for (c = 0; c < Nc; ++c) {
for (d = 0; d < dimRef; ++d) {
f1[(q*Nc + c)*dimRef+d] = 0.0;
for (e = 0; e < dimReal; ++e) f1[(q*Nc + c)*dimRef+d] += invJ[d*dimReal+e]*refSpaceDer[c*dimReal+e];
f1[(q*Nc + c)*dimRef+d] *= detJ*quadWeights[q];
const PetscReal w = detJ*quadWeights[q];
PetscInt c, d, e;

if (f0) for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= w;
if (f1) {
for (c = 0; c < Nc; ++c) {
for (d = 0; d < dim; ++d) {
f1[(q*Nc + c)*dim+d] = 0.0;
for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
f1[(q*Nc + c)*dim+d] *= w;
}
}
}
#if 0
if (debug > 1) {
for (c = 0; c < Nc; ++c) {
ierr = PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));CHKERRQ(ierr);
for (d = 0; d < dimRef; ++d) {
ierr = PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dimRef+d]));CHKERRQ(ierr);
for (d = 0; d < dim; ++d) {
ierr = PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));CHKERRQ(ierr);
}
}
}
Expand Down

0 comments on commit 66b078f

Please sign in to comment.