#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: disc.c,v 1.3 2000/01/10 03:54:16 knepley Exp $"; #endif /* This is the interface to discretization objects */ #include "src/grid/discretization/discimpl.h" /*I "discretization.h" I*/ /* Logging support */ int DISCRETIZATION_COOKIE; int DISCRETIZATION_EvaluateOperatorGalerkin; /*------------------------------------------------- Generic Functions ------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationDestroy" /*@ DiscretizationDestroy - Destroys a discretization object. Collective on Discretization Input Parameter: . disc - The discretization Level: beginner .keywords: discretization, destroy .seealso: DiscretizationView() @*/ int DiscretizationDestroy(Discretization disc) { int op, arg; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); if (--disc->refct > 0) PetscFunctionReturn(0); if (disc->bdDisc != PETSC_NULL) { ierr = DiscretizationDestroy(disc->bdDisc); CHKERRQ(ierr); } ierr = PetscFree(disc->quadPoints); CHKERRQ(ierr); ierr = PetscFree(disc->quadWeights); CHKERRQ(ierr); ierr = PetscFree(disc->quadShapeFuncs); CHKERRQ(ierr); ierr = PetscFree(disc->quadShapeFuncDers); CHKERRQ(ierr); ierr = PetscFree(disc->funcVal); CHKERRQ(ierr); for(op = 0; op < disc->numOps; op++) { if (disc->operators[op]->precompInt) { ierr = PetscFree(disc->operators[op]->precompInt); CHKERRQ(ierr); } ierr = PetscFree(disc->operators[op]); CHKERRQ(ierr); } ierr = PetscFree(disc->operators); CHKERRQ(ierr); for(arg = 0; arg < 2; arg++) { ierr = PetscFree(disc->fieldVal[arg]); CHKERRQ(ierr); } ierr = PetscFree(disc->fieldVal); CHKERRQ(ierr); ierr = (*disc->ops->destroy)(disc); CHKERRQ(ierr); PetscLogObjectDestroy(disc); PetscHeaderDestroy(disc); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationView" /*@ DiscretizationView - Views a discretization object. Collective on Discretization Input Parameters: + disc - The disc context to distroy - viewer - The viewer Level: beginner .keywords: discretization, view .seealso: DiscretizationDestroy() @*/ int DiscretizationView(Discretization disc, PetscViewer viewer) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); if (viewer == PETSC_NULL) { viewer = PETSC_VIEWER_STDOUT_SELF; } else { PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); } ierr = (*disc->ops->view)(disc, viewer); CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "DiscretizationSerialize_Generic" int DiscretizationSerialize_Generic(MPI_Comm comm, Discretization *disc, PetscViewer viewer, PetscTruth store) { Discretization d; int fd; int len; char *type_name = PETSC_NULL; PetscTruth match; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); PetscValidPointer(disc); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match); CHKERRQ(ierr); if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer"); ierr = PetscViewerBinaryGetDescriptor(viewer, &fd); CHKERRQ(ierr); if (store) { d = *disc; ierr = PetscStrlen(d->type_name, &len); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &len, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, d->type_name, len, PETSC_CHAR, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &d->comp, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &d->field, 1, PETSC_INT, 0); CHKERRQ(ierr); } else { ierr = DiscretizationCreate(comm, &d); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &len, 1, PETSC_INT); CHKERRQ(ierr); if (len) { ierr = PetscMalloc((len+1) * sizeof(char), &type_name); CHKERRQ(ierr); type_name[len] = 0; } ierr = PetscBinaryRead(fd, type_name, len, PETSC_CHAR); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &d->comp, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &d->field, 1, PETSC_INT); CHKERRQ(ierr); ierr = DiscretizationSetType(d, type_name); CHKERRQ(ierr); if (len) { ierr = PetscFree(type_name); CHKERRQ(ierr); } *disc = d; } PetscFunctionReturn(0); } EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetTypeFromOptions" /* DiscretizationSetTypeFromOptions - Sets the type of discretization from user options. Collective on Discretization Input Parameter: . disc - The discretization Level: intermediate .keywords: Discretization, set, options, database, type .seealso: DiscretizationSetFromOptions(), DiscretizationSetType() */ static int DiscretizationSetTypeFromOptions(Discretization disc) { PetscTruth opt; char *defaultType; char typeName[256]; int ierr; PetscFunctionBegin; if (disc->type_name != PETSC_NULL) { defaultType = disc->type_name; } else { defaultType = DISCRETIZATION_TRIANGULAR_2D_LINEAR; } if (DiscretizationRegisterAllCalled == PETSC_FALSE) { ierr = DiscretizationRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsList("-disc_type", "Discretization method"," DiscretizationSetType", DiscretizationList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = DiscretizationSetType(disc, typeName); CHKERRQ(ierr); } else { ierr = DiscretizationSetType(disc, defaultType); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetSerializeTypeFromOptions" /* DiscretizationSetSerializeTypeFromOptions - Sets the type of discretization serialization from user options. Collective on Discretization Input Parameter: . disc - The discretization Level: intermediate .keywords: Discretization, set, options, database, type .seealso: DiscretizationSetFromOptions(), DiscretizationSetSerializeType() */ static int DiscretizationSetSerializeTypeFromOptions(Discretization disc) { PetscTruth opt; char *defaultType; char typeName[256]; int ierr; PetscFunctionBegin; if (disc->serialize_name != PETSC_NULL) { defaultType = disc->serialize_name; } else { defaultType = DISCRETIZATION_SER_GENERIC; } if (DiscretizationSerializeRegisterAllCalled == PETSC_FALSE) { ierr = DiscretizationSerializeRegisterAll(PETSC_NULL); CHKERRQ(ierr); } ierr = PetscOptionsList("-disc_serialize_type", "Discretization serialization method"," DiscretizationSetSerializeType", DiscretizationSerializeList, defaultType, typeName, 256, &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = DiscretizationSetSerializeType(disc, typeName); CHKERRQ(ierr); } else { ierr = DiscretizationSetSerializeType(disc, defaultType); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetFromOptions" /*@ DiscretizationSetFromOptions - Sets various Discretization parameters from user options. Collective on Discretization Input Parameter: . disc - The discretization Notes: To see all options, run your program with the -help option, or consult the users manual. Must be called after DiscretizationCreate() but before the Discretization is used. Level: intermediate .keywords: Discretization, set, options, database .seealso: DiscretizationCreate(), DiscretizationPrintHelp(), DiscretizationSetOptionsPrefix() @*/ int DiscretizationSetFromOptions(Discretization disc) { Discretization bdDisc; PetscTruth opt; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); ierr = PetscOptionsBegin(disc->comm, disc->prefix, "Discretization options", "Discretization"); CHKERRQ(ierr); /* Handle generic disc options */ ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = DiscretizationPrintHelp(disc); CHKERRQ(ierr); } /* Handle disc type options */ ierr = DiscretizationSetTypeFromOptions(disc); CHKERRQ(ierr); /* Handle disc serialization options */ ierr = DiscretizationSetSerializeTypeFromOptions(disc); CHKERRQ(ierr); /* Handle specific disc options */ if (disc->ops->setfromoptions != PETSC_NULL) { ierr = (*disc->ops->setfromoptions)(disc); CHKERRQ(ierr); } ierr = PetscOptionsEnd(); CHKERRQ(ierr); /* Handle subobject options */ ierr = DiscretizationGetBoundaryDiscretization(disc, &bdDisc); CHKERRQ(ierr); if (bdDisc != PETSC_NULL) { ierr = DiscretizationSetFromOptions(bdDisc); CHKERRQ(ierr); } ierr = DiscretizationViewFromOptions(disc, disc->name); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationViewFromOptions" /*@ DiscretizationViewFromOptions - This function visualizes the discretization based upon user options. Collective on Discretization Input Parameter: . disc - The discretization Level: intermediate .keywords: Discretization, view, options, database .seealso: DiscretizationSetFromOptions(), DiscretizationView() @*/ int DiscretizationViewFromOptions(Discretization disc, char *title) { PetscViewer viewer; PetscDraw draw; PetscTruth opt; char *titleStr; char typeName[1024]; char fileName[PETSC_MAX_PATH_LEN]; int len; int ierr; PetscFunctionBegin; ierr = PetscOptionsHasName(disc->prefix, "-disc_view", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscOptionsGetString(disc->prefix, "-disc_view", typeName, 1024, &opt); CHKERRQ(ierr); ierr = PetscStrlen(typeName, &len); CHKERRQ(ierr); if (len > 0) { ierr = PetscViewerCreate(disc->comm, &viewer); CHKERRQ(ierr); ierr = PetscViewerSetType(viewer, typeName); CHKERRQ(ierr); ierr = PetscOptionsGetString(disc->prefix, "-disc_view_file", fileName, 1024, &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerSetFilename(viewer, fileName); CHKERRQ(ierr); } else { ierr = PetscViewerSetFilename(viewer, disc->name); CHKERRQ(ierr); } ierr = DiscretizationView(disc, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } else { ierr = DiscretizationView(disc, PETSC_NULL); CHKERRQ(ierr); } } ierr = PetscOptionsHasName(disc->prefix, "-disc_view_draw", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { ierr = PetscViewerDrawOpen(disc->comm, 0, 0, 0, 0, 300, 300, &viewer); CHKERRQ(ierr); ierr = PetscViewerDrawGetDraw(viewer, 0, &draw); CHKERRQ(ierr); if (title != PETSC_NULL) { titleStr = title; } else { ierr = PetscObjectName((PetscObject) disc); CHKERRQ(ierr) ; titleStr = disc->name; } ierr = PetscDrawSetTitle(draw, titleStr); CHKERRQ(ierr); ierr = DiscretizationView(disc, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = PetscDrawPause(draw); CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationPrintHelp" /*@ DiscretizationPrintHelp - Prints all options for the discretization. Input Parameter: . disc - The discretization Options Database Keys: $ -help, -h Level: intermediate .keywords: Discretization, help .seealso: DiscretizationSetFromOptions() @*/ int DiscretizationPrintHelp(Discretization disc) { char p[64]; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); ierr = PetscStrcpy(p, "-"); CHKERRQ(ierr); if (disc->prefix != PETSC_NULL) { ierr = PetscStrcat(p, disc->prefix); CHKERRQ(ierr); } (*PetscHelpPrintf)(disc->comm, "Discretization options ------------------------------------------------\n"); (*PetscHelpPrintf)(disc->comm," %sdisc_type : Sets the discretization type\n", p); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetOptionsPrefix" /*@C DiscretizationSetOptionsPrefix - Sets the prefix used for searching for all discretization options in the database. Input Parameters: + disc - The discretization - prefix - The prefix to prepend to all option names Notes: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. Level: intermediate .keywords: Discretization, set, options, prefix, database .seealso: DiscretizationSetFromOptions() @*/ int DiscretizationSetOptionsPrefix(Discretization disc, char *prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); ierr = PetscObjectSetOptionsPrefix((PetscObject) disc, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationAppendOptionsPrefix" /*@C DiscretizationAppendOptionsPrefix - Appends to the prefix used for searching for all discretization options in the database. Collective on Discretization Input Parameters: + disc - The discretization - prefix - The prefix to prepend to all option names Notes: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. Level: intermediate .keywords: Discretization, append, options, prefix, database .seealso: DiscretizationGetOptionsPrefix() @*/ int DiscretizationAppendOptionsPrefix(Discretization disc, char *prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); ierr = PetscObjectAppendOptionsPrefix((PetscObject) disc, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetOptionsPrefix" /*@C DiscretizationGetOptionsPrefix - Sets the prefix used for searching for all discretization options in the database. Input Parameter: . disc - The discretization Output Parameter: . prefix - A pointer to the prefix string used Level: intermediate .keywords: Discretization, get, options, prefix, database .seealso: DiscretizationAppendOptionsPrefix() @*/ int DiscretizationGetOptionsPrefix(Discretization disc, char **prefix) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); ierr = PetscObjectGetOptionsPrefix((PetscObject) disc, prefix); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetupDefaultOperators" /*@ DiscretizationSetupDefaultOperators - Adds the default set of operators to the discretization if they are not already present. Input Parameter: . disc - The discretization Level: intermediate .keywords: Discretization, default, operator .seealso: DiscretizationSetup(), DiscretizationSerialize() @*/ int DiscretizationSetupDefaultOperators(Discretization disc) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); if (disc->numOps > 0) PetscFunctionReturn(0); ierr = (*disc->ops->setupoperators)(disc); CHKERRQ(ierr); PetscFunctionReturn(0); } /*-------------------------------------------------- Query Functions -------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetNumComponents" /*@C DiscretizationSetNumComponents - This function sets the number of components in the field being discretized. Not collective Input Parameters: . disc - The discretization . numComp - The number of components in the field Level: intermediate .keywords discretization, component, field .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationSetNumComponents(Discretization disc, int numComp) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); if (numComp <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Discretization msut have at least one component"); disc->comp = numComp; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetNumComponents" /*@C DiscretizationGetNumComponents - This function returns the number of components in the field being discretized. Not collective Input Parameter: . disc - The discretization Output Parameter: . numComp - The number of components in the field Level: intermediate .keywords discretization, component, field .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetNumComponents(Discretization disc, int *numComp) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(numComp); *numComp = disc->comp; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetNumFunctions" /*@C DiscretizationGetNumFunctions - This function returns the number of functions per element in the Discretization. Not collective Input Parameter: . disc - The discretization Output Parameter: . numFunc - The number of functions per element Level: intermediate .keywords discretization, functions .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetNumFunctions(Discretization disc, int *numFunc) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(numFunc); *numFunc = disc->funcs; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetBoundaryDiscretization" /*@C DiscretizationGetBoundaryDiscretization - This function returns boundary discretization, or PETSC_NULL if none exists. Not collective Input Parameter: . disc - The discretization Output Parameter: . bdDisc - The boundary discretization, or PETSC_NULL Level: intermediate .keywords discretization, functions .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetBoundaryDiscretization(Discretization disc, Discretization *bdDisc) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidPointer(bdDisc); *bdDisc = disc->bdDisc; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationSetField" /*@C DiscretizationSetField - This function sets the corresponding field in the Grid for this Discretization. Not collective Input Parameters: . disc - The discretization . field - The corresponding field in the Grid Level: intermediate .keywords Discretization, component, field .seealso DiscretizationGetField(), DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationSetField(Discretization disc, int field) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); /* Could possibly validate with an associated Grid */ disc->field = field; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetField" /*@C DiscretizationGetField - This function gets the corresponding field in the Grid for this Discretization. Not collective Input Parameter: . disc - The discretization Output Parameter: . field - The corresponding field in the Grid Level: intermediate .keywords discretization, component, field .seealso DiscretizationSetField(), DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetField(Discretization disc, int *field) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(field); *field = disc->field; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetFieldValues" /*@C DiscretizationGetFieldValues - This function gets the values for the corresponding field in the Grid for this Discretization. Not collective Input Parameter: . disc - The discretization Output Parameter: . fieldValues - The values for the corresponding field in the Grid Level: intermediate .keywords discretization, component, field .seealso DiscretizationSetField(), DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetFieldValues(Discretization disc, PetscScalar **fieldValues) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidScalarPointer(*fieldValues); *fieldValues = disc->fieldVal[0]; PetscFunctionReturn(0); } /*--------------------------------------------- Quadrature Query Functions -------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetNumQuadraturePoints" /*@C DiscretizationGetNumQuadraturePoints - This function returns the number of points used in the element quadrature. Not collective Input Parameter: . disc - The discretization Output Parameter: . numPoints - The number of quadrature points Level: intermediate .keywords discretization, quadrature, points .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetNumQuadraturePoints(Discretization disc, int *numPoints) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(numPoints); *numPoints = disc->numQuadPoints; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetQuadraturePoints" /*@C DiscretizationGetQuadraturePoints - This function returns the coordinates of each quadrature point. Not collective Input Parameter: . disc - The discretization Output Parameter: . coords - The quadrature point coordinates Level: intermediate .keywords discretization, quadrature, coordinates, points .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetQuadraturePoints(Discretization disc, double **coords) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidPointer(coords); *coords = disc->quadPoints; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetQuadratureWeights" /*@C DiscretizationGetQuadratureWeights - This function returns the weight associated with each quadrature point. Not collective Input Parameter: . disc - The discretization Output Parameter: . weights - The quadrature weights Level: intermediate .keywords discretization, quadrature, weights .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetQuadratureWeights(Discretization disc, double **weights) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidPointer(weights); *weights = disc->quadWeights; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetQuadratureFunctions" /*@C DiscretizationGetQuadratureFunctions - This function returns each shape function evaluated at the quadrature points. Not collective Input Parameter: . disc - The discretization Output Parameter: . funcs - The shape functions evaluated at the quadrature points Level: intermediate .keywords discretization, quadrature, functions .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetQuadratureFunctions(Discretization disc, double **funcs) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidPointer(funcs); *funcs = disc->quadShapeFuncs; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationGetQuadratureDerivatives" /*@C DiscretizationGetQuadratureDerivatives - This function returns the derivatives of each shape function evaluated at the quadrature points. Not collective Input Parameter: . disc - The discretization Output Parameter: . ders - The shape function derivatives evaluated at the quadrature points Level: intermediate .keywords discretization, quadrature, derivatives .seealso DiscretizationCreate(), DiscretizationRegisterOperator() @*/ int DiscretizationGetQuadratureDerivatives(Discretization disc, double **ders) { PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidPointer(ders); *ders = disc->quadShapeFuncDers; PetscFunctionReturn(0); } /*----------------------------------------------- Registration Functions ---------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationRegisterOperator_Private" static int DiscretizationRegisterOperator_Private(Discretization disc, PetscScalar *precompInt, OperatorFunction func, ALEOperatorFunction ALEFunc, int *newOp) { Operator *tempOperators; Operator op; int ierr; PetscFunctionBegin; while (disc->numOps >= disc->maxOps) { ierr = PetscMalloc(disc->maxOps*2 * sizeof(Operator), &tempOperators); CHKERRQ(ierr); PetscLogObjectMemory(disc, disc->maxOps * sizeof(Operator)); ierr = PetscMemcpy(tempOperators, disc->operators, disc->maxOps * sizeof(Operator)); CHKERRQ(ierr); ierr = PetscFree(disc->operators); CHKERRQ(ierr); disc->operators = tempOperators; disc->maxOps *= 2; } /* Create the new operator */ ierr = PetscNew(struct _Operator, &op); CHKERRQ(ierr); PetscLogObjectMemory(disc, sizeof(struct _Operator)); op->precompInt = precompInt; op->test = PETSC_NULL; op->opFunc = func; op->ALEOpFunc = ALEFunc; disc->operators[disc->numOps] = op; *newOp = disc->numOps; disc->numOps++; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationRegisterOperator" /*@C DiscretizationRegisterOperator - This function defines a new operator Collective on Discretization Input Parameters: + disc - The discretization - func - The function defining the operator Output Parameter: . newOp - The index of the new operator Level: advanced .keywords Discretization, operator, register .seealso DiscretizationRegisterALEOperator(), GridAddMatOperator(), GridAddRhsOperator() @*/ int DiscretizationRegisterOperator(Discretization disc, OperatorFunction func, int *newOp) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(newOp); ierr = DiscretizationRegisterOperator_Private(disc, PETSC_NULL, func, PETSC_NULL, newOp); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationRegisterALEOperator" /*@C DiscretizationRegisterALEOperator - This function defines a new ALE operator Collective on Discretization Input Parameters: + disc - The discretization - func - The function deinfing the operator Output Parameter: . newOp - The index of the new operator Level: advanced .keywords Discretization, operator, register .seealso DiscretizationRegisterOperator(), GridAddMatOperator(), GridAddRhsOperator() @*/ int DiscretizationRegisterALEOperator(Discretization disc, ALEOperatorFunction func, int *newOp) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidIntPointer(newOp); ierr = DiscretizationRegisterOperator_Private(disc, PETSC_NULL, PETSC_NULL, func, newOp); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationRegisterPrecomputedOperator" /*@C DiscretizationRegisterPrecomputedOperator - This function defines a new operator for which the integral can be precomputed . This requires knowledge of the discretization itself and the storage format, and thus is only useful for developers. Collective on Discretization Input Parameters: + disc - The discretization - precompInt - The array of precomputed values Output Parameter: . newOp - The index of the new operator Level: developer .keywords Discretization, operator, register .seealso DiscretizationRegisterALEOperator(), GridAddMatOperator(), GridAddRhsOperator() @*/ int DiscretizationRegisterPrecomputedOperator(Discretization disc, PetscScalar *precompInt, int *newOp) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidScalarPointer(precompInt); PetscValidIntPointer(newOp); ierr = DiscretizationRegisterOperator_Private(disc, precompInt, PETSC_NULL, PETSC_NULL, newOp); CHKERRQ(ierr); PetscFunctionReturn(0); } /*------------------------------------------------ Evaluation Functions ----------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateFunctionGalerkin" /*@ DiscretizationEvaluateFunctionGalerkin - This function calculates the weak form of a function over a single element. Collective on Discretization Input Parameters: + disc - The discretization . mesh - The associated mesh . f - The PointFunction of which we want the weak form . alpha - A scalar multiplier . elem - The local element number - ctx - The user-supplied context Output Parameter: . v - The element vector for the given element Level: beginner .keywords: discretization, function, weak form .seealso: DiscretizationEvaluateOperatorGalerkin(), DiscretizationEvaluateNonlinearOperatorGalerkin() @*/ int DiscretizationEvaluateFunctionGalerkin(Discretization disc, Mesh mesh, PointFunction f, PetscScalar alpha, int elem, PetscScalar *v, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidScalarPointer(v); ierr = (*(disc->ops->evaluatefunctiongalerkin))(disc, mesh, f, alpha, elem, v, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateOperatorGalerkin" /*@ DiscretizationEvaluateOperatorGalerkin - This function calculates the weak form of an operator over a single element. Collective on Discretization Input Parameters: + disc - The discretization for the basis functions . mesh - The associated mesh . elemSize - The size of the element matrix . rowStart - The starting row index in the element matrix . colStart - The starting column index in the element matrix . op - The operator index (of registered operators) . alpha - The scalar multiple of the operator . elem - The local element number . field - The field values - ctx - The user-supplied context Output Parameter: . mat - The element matrix for the given element Level: beginner .keywords: discretization, operator, weak form .seealso: DiscretizationEvaluateFunctionGalerkin(), DiscretizationEvaluateALEOperatorGalerkin(), DiscretizationEvaluateOperatorGalerkinMF(), DiscretizationEvaluateNonlinearOperatorGalerkin() @*/ int DiscretizationEvaluateOperatorGalerkin(Discretization disc, Mesh mesh, int elemSize, int rowStart, int colStart, int op, PetscScalar alpha, int elem, PetscScalar *field, PetscScalar *mat, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); /* field may be PETSC_NULL */ PetscValidScalarPointer(mat); ierr = (*(disc->ops->evaluateoperatorgalerkin))(disc, mesh, elemSize, rowStart, colStart, op, alpha, elem, field, mat, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateALEOperatorGalerkin" /*@ DiscretizationEvaluateALEOperatorGalerkin - This function calculates the weak form of an operator over a single element. Collective on Discretization Input Parameters: + disc - The discretization for the basis functions . mesh - The associated mesh . elemSize - The size of the element matrix . rowStart - The starting row index in the element matrix . colStart - The starting column index in the element matrix . op - The operator index (of registered operators) . alpha - The scalar multiple of the operator . elem - The local element number . field - The field values . ALEfield - The mesh velocity field values - ctx - The user-supplied context Output Parameter: . mat - The element matrix for the given element Level: beginner .keywords: discretization, operator, ALE, weak form .seealso: DiscretizationEvaluateOperatorGalerkin() @*/ int DiscretizationEvaluateALEOperatorGalerkin(Discretization disc, Mesh mesh, int elemSize, int rowStart, int colStart, int op, PetscScalar alpha, int elem, PetscScalar *field, PetscScalar *ALEfield, PetscScalar *mat, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); /* field may be PETSC_NULL */ PetscValidScalarPointer(ALEfield); PetscValidScalarPointer(mat); ierr = (*(disc->ops->evaluatealeoperatorgalerkin))(disc, mesh, elemSize, rowStart, colStart, op, alpha, elem, field, ALEfield, mat, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateOperatorGalerkinMF" /*@ DiscretizationEvaluateOperatorGalerkinMF - This function calculates the weak form of the action of an operator over a single element. Collective on Discretization Input Parameters: + disc - The discretization . mesh - The associated mesh . elemSize - The size of the element matrix . rowStart - The starting row index in the element matrix . colStart - The starting column index in the element matrix . op - The operator index (of registered operators) . alpha - The scalar multiple of the operator . elem - The local element number . field - The element vector for the grid vector acted upon . app - The element vector to which to apply the matrix, usually identical to field . mat - A temporary element matrix for work space - ctx - The user-supplied context Output Parameter: . vec - The element vector for the given element Level: beginner .keywords: discretization, operator, MF, weak form .seealso: DiscretizationEvaluateOperatorGalerkin() @*/ int DiscretizationEvaluateOperatorGalerkinMF(Discretization disc, Mesh mesh, int elemSize, int rowStart, int colStart, int op, PetscScalar alpha, int elem, PetscScalar *field, PetscScalar *app, PetscScalar *vec, PetscScalar *mat, void *ctx) { Discretization test; /* The space of test functions */ int rowSize; /* The number of shape functions per element */ int colSize; /* The number of test functions per element */ int i, j; int ierr; PetscFunctionBegin; /* Get discretization info */ PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); /* field may be PETSC_NULL */ PetscValidScalarPointer(vec); PetscValidScalarPointer(mat); test = disc->operators[op]->test; PetscValidHeaderSpecific(test, DISCRETIZATION_COOKIE); rowSize = test->size; colSize = disc->size; ierr = DiscretizationEvaluateOperatorGalerkin(disc, mesh, elemSize, rowStart, colStart, op, alpha, elem, field, mat, ctx); CHKERRQ(ierr); /* Perform the local element matvec */ for(i = 0; i < rowSize; i++) for(j = 0; j < colSize; j++) vec[i+rowStart] += mat[(i+rowStart)*elemSize+j+colStart]*app[j]; PetscLogFlops(2*rowSize*colSize); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateALEOperatorGalerkinMF" /*@ DiscretizationEvaluateALEOperatorGalerkinMF - This function calculates the weak form of the action of an operator over a single element. Collective on Discretization Input Parameters: + disc - The discretization . mesh - The associated mesh . elemSize - The size of the element matrix . rowStart - The starting row index in the element matrix . colStart - The starting column index in the element matrix . op - The operator index (of registered operators) . alpha - The scalar multiple of the operator . elem - The local element number . field - The element vector for the grid vector acted upon . app - The element vector to which to apply the matrix, usually identical to field . ALEfield - The element vector for the ALE velocity field . mat - A temporary element matrix for work space - ctx - The user-supplied context Output Parameter: . vec - The element vector for the given element Level: beginner .keywords: discretization, operator, ALE, MF, weak form .seealso: DiscretizationEvaluateOperatorGalerkin() @*/ int DiscretizationEvaluateALEOperatorGalerkinMF(Discretization disc, Mesh mesh, int elemSize, int rowStart, int colStart, int op, PetscScalar alpha, int elem, PetscScalar *field, PetscScalar *app, PetscScalar *ALEfield, PetscScalar *vec, PetscScalar *mat, void *ctx) { Discretization test; /* The space of test functions */ int rowSize; /* The number of shape functions per element */ int colSize; /* The number of test functions per element */ int i, j; int ierr; PetscFunctionBegin; /* Get discretization info */ PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); /* field may be PETSC_NULL */ PetscValidScalarPointer(ALEfield); PetscValidScalarPointer(vec); PetscValidScalarPointer(mat); test = disc->operators[op]->test; PetscValidHeaderSpecific(test, DISCRETIZATION_COOKIE); rowSize = test->size; colSize = disc->size; ierr = DiscretizationEvaluateALEOperatorGalerkin(disc, mesh, elemSize, rowStart, colStart, op, alpha, elem, field, ALEfield, mat, ctx); CHKERRQ(ierr); /* Perform the local element matvec */ for(i = 0; i < rowSize; i++) for(j = 0; j < colSize; j++) vec[i+rowStart] += mat[(i+rowStart)*elemSize+j+colStart]*app[j]; PetscLogFlops(2*rowSize*colSize); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateNonlinearOperatorGalerkin" /*@ DiscretizationEvaluateNonlinearOperatorGalerkin - This function calculates the weak form of a nonlinear operator over a single element. Collective on Discretization Input Parameters: + disc - The discretization . mesh - The associated mesh . op - The function defining the nonlinear operator . alpha - A scalar multiplier . elem - The local element number . numArgs - THe number of input element vectors . field - The element vectors for the grid vectors acted upon - ctx - The user-supplied context Output Parameter: . vec - The global element vector for the given element Level: beginner .keywords: discretization, operator, nonlinear, weak form .seealso: DiscretizationEvaluateFunctionGalerkin(), DiscretizationEvaluateOperatorGalerkin(), DiscretizationEvaluateNonlinearALEOperatorGalerkin() @*/ int DiscretizationEvaluateNonlinearOperatorGalerkin(Discretization disc, Mesh mesh, NonlinearOperator op, PetscScalar alpha, int elem, int numArgs, PetscScalar **field, PetscScalar *vec, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidScalarPointer(field); PetscValidScalarPointer(vec); ierr = (*(disc->ops->evaluatenonlinearoperatorgalerkin))(disc, mesh, op, alpha, elem, numArgs, field, vec, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationEvaluateNonlinearALEOperatorGalerkin" /*@ DiscretizationEvaluateNonlinearALEOperatorGalerkin - This function calculates the weak form of a nonlinear operator over a single element. Input Parameters: + disc - The discretization . mesh - The associated mesh . op - The function defining the nonlinear operator . alpha - A scalar multiplier . elem - The local element number . numArgs - THe number of input element vectors . field - The element vectors for the grid vectors acted upon . ALEfield - The element vector for the mesh velocity field - ctx - The user-supplied context Output Parameter: . vec - The global element vector for the given element Level: beginner .keywords: discretization, operator, nonlinear, ALE, weak form .seealso: DiscretizationEvaluateFunctionGalerkin(), DiscretizationEvaluateOperatorGalerkin(), DiscretizationEvaluateNonlinearOperatorGalerkin() @*/ int DiscretizationEvaluateNonlinearALEOperatorGalerkin(Discretization disc, Mesh mesh, NonlinearOperator op, PetscScalar alpha, int elem, int numArgs, PetscScalar **field, PetscScalar *ALEfield, PetscScalar *vec, void *ctx) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidScalarPointer(field); PetscValidScalarPointer(ALEfield); PetscValidScalarPointer(vec); ierr = (*(disc->ops->evaluatenonlinearaleoperatorgalerkin))(disc, mesh, op, alpha, elem, numArgs, field, ALEfield, vec, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } /*---------------------------------------------- Interpolation Functions ---------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "DiscretizationInterpolateField" /*@ DiscretizationInterpolateField - This function interpolates a field at a given point. Collective on Discretization Input Parameters: + disc - The discretization . mesh - The associated mesh . elem - The element containing the point . x,y,z - The interpolation point . oldFieldVal - The element vector for 'elem' - type - The method of interpolation, e.g. INTERPOLATION_LOCAL Output Parameter: . newFieldVal - The field components at the given point Level: beginner .keywords: discretization, interpolation .seealso: DiscretizationInterpolateElementVec() @*/ int DiscretizationInterpolateField(Discretization disc, Mesh mesh, int elem, double x, double y, double z, PetscScalar *oldFieldVal, PetscScalar *newFieldVal, InterpolationType type) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(mesh, MESH_COOKIE); PetscValidScalarPointer(oldFieldVal); PetscValidScalarPointer(newFieldVal); ierr = (*disc->ops->interpolatefield)(disc, mesh, elem, x, y, z, oldFieldVal, newFieldVal, type); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "DiscretizationInterpolateElementVec" /*@ DiscretizationInterpolateElementVec - Interpolates a given element vector from one discretization to another. Input Parameters: + disc - The original discretization . vec - The original element vector - newDisc - The discretization defining the new element vector Output Parameter: . newVec - The interpolated element vector Level: beginner .keywords: discretization, interpolation .seealso: DiscretizationInterpolateField() @*/ int DiscretizationInterpolateElementVec(Discretization disc, ElementVec vec, Discretization newDisc, ElementVec newVec) { int size, newSize; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(disc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(newDisc, DISCRETIZATION_COOKIE); PetscValidHeaderSpecific(vec, ELEMENT_VEC_COOKIE); PetscValidHeaderSpecific(newVec, ELEMENT_VEC_COOKIE); if (disc->comp != newDisc->comp) SETERRQ(PETSC_ERR_ARG_INCOMP, "Fields must have the same number of components"); ierr = ElementVecGetSize(vec, &size); CHKERRQ(ierr); ierr = ElementVecGetSize(newVec, &newSize); CHKERRQ(ierr); if ((disc->size > size) || (newDisc->size > newSize)) SETERRQ(PETSC_ERR_ARG_WRONG, "Element vector is too small to contain field"); ierr = (*disc->ops->interpolateelemvec)(disc, vec, newDisc, newVec); CHKERRQ(ierr); PetscFunctionReturn(0); }