#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: gvec2dView.c,v 1.22 2000/10/08 00:27:05 knepley Exp $"; #endif #include "src/gvec/gvecimpl.h" /*I "gvec.h" I*/ #include "gvec2dView.h" extern int ViewerSiloCheckMesh(PetscViewer, Mesh); #undef __FUNCT__ #define __FUNCT__ "GVecView_Triangular_2D_Draw" int GVecView_Triangular_2D_Draw(GVec gvec, PetscViewer v) { int ierr; PetscFunctionBegin; ierr = VecView(gvec, v); PetscFunctionReturn(ierr); } #ifdef PETSC_HAVE_SILO #include "src/sys/src/viewer/impls/silo/vsilo.h" #undef __FUNCT__ #define __FUNCT__ "GVecView_Triangular_2D_Silo" int GVecView_Triangular_2D_Silo(GVec gvec, PetscViewer viewer) { Viewer_Silo *vsilo = (Viewer_Silo *) viewer->data; Grid grid; Mesh mesh; Mesh_Triangular *tri; VarOrdering order; FieldClassMap cm; Vec ghostVec; VecScatter ghostScatter; ElementVec elemVec; PetscTruth reduceSystem; PetscTruth reduceElement; int numElements, numCorners; int *elements; char **subNames; float **newArrays; int **localStart; PetscScalar *array; DBfile *fp; DBoptlist *optlist; int *classes; int *offsets; char *fieldName; char name[256]; char buf[1024]; int numFields, numNodes; int *elemStart; int size, elemSize, comp, len; int f, field, elem, func, node, c, index; int ierr; PetscFunctionBegin; ierr = GVecGetGrid(gvec, &grid); CHKERRQ(ierr); ierr = GVecGetOrder(gvec, &order); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(order, &cm); CHKERRQ(ierr); ierr = VecGetSize(gvec, &size); CHKERRQ(ierr); offsets = order->offsets; localStart = order->localStart; numFields = cm->numFields; classes = cm->classes; ierr = ViewerSiloGetFilePointer(viewer, &fp); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = ViewerSiloCheckMesh(viewer, mesh); CHKERRQ(ierr); ierr = MeshGetInfo(mesh, PETSC_NULL, &numNodes, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr); /* Construct names */ for(f = 0; f < numFields; f++) { field = cm->fields[f]; ierr = GridGetFieldComponents(grid, field, &comp); CHKERRQ(ierr); ierr = GridGetFieldName(grid, field, &fieldName); CHKERRQ(ierr); /* ierr = GridGetFieldUnits(grid, field, &fieldUnits); CHKERRQ(ierr); */ if (vsilo->objName != PETSC_NULL) { ierr = PetscStrncpy(name, vsilo->objName, 240); CHKERRQ(ierr); } else { ierr = PetscStrcpy(name, ""); CHKERRQ(ierr); } if (fieldName == PETSC_NULL) { sprintf(buf, "field%d", field); ierr = PetscStrcat(name, buf); CHKERRQ(ierr); } else { ierr = PetscStrcat(name, fieldName); CHKERRQ(ierr); } ierr = PetscMalloc(comp * sizeof(char *), &subNames); CHKERRQ(ierr); ierr = PetscMalloc(comp * sizeof(float *), &newArrays); CHKERRQ(ierr); ierr = PetscStrlen(name, &len); CHKERRQ(ierr); len += (int) log10((double) comp) + 6; for(c = 0; c < comp; c++) { ierr = PetscMalloc(len * sizeof(char), &subNames[c]); CHKERRQ(ierr); ierr = PetscMalloc(numNodes * sizeof(float), &newArrays[c]); CHKERRQ(ierr); sprintf(subNames[c], "%scomp%d", name, c); } #if 1 /* Setup reduction */ ierr = (*grid->ops->gridsetupghostscatter)(grid, order, &ghostVec, &ghostScatter); CHKERRQ(ierr); /* Fill the local vector */ ierr = GridGlobalToLocalGeneral(grid, gvec, ghostVec, INSERT_VALUES, ghostScatter); CHKERRQ(ierr); /* Setup element vector and matrix */ elemSize = grid->locOrder->elemSize; elemStart = grid->locOrder->elemStart; if (cm->isConstrained == PETSC_TRUE) { for(f = 0; f < numFields; f++) { if (grid->isFieldConstrained[cm->fields[f]]) elemSize += grid->disc[cm->fields[f]]->funcs*grid->constraintCompDiff[cm->fields[f]]; } } ierr = ElementVecCreate(gvec->comm, elemSize, &elemVec); CHKERRQ(ierr); array = elemVec->array; reduceSystem = grid->reduceSystem; reduceElement = grid->reduceElement; tri = (Mesh_Triangular *) mesh->data; numElements = tri->numFaces; numCorners = tri->numCorners; elements = tri->faces; for(elem = 0; elem < numElements; elem++) { /* Initialize element vector */ ierr = ElementVecZero(elemVec); CHKERRQ(ierr); elemVec->reduceSize = grid->locOrder->elemSize; /* Setup local row indices */ ierr = GridCalcGeneralElementVecIndices(grid, elem, order, PETSC_TRUE, elemVec); CHKERRQ(ierr); /* Setup local vectors */ ierr = GridLocalToElementGeneral(grid, ghostVec, reduceSystem, reduceElement, elemVec); CHKERRQ(ierr); /* Must transform to unconstrained variables */ ierr = GridProjectElementVec(grid, mesh, elem, order, PETSC_FALSE, elemVec); CHKERRQ(ierr); /* Put values into new arrays */ index = elemStart[field]; for(func = 0; func < grid->disc[field]->funcs; func++) { node = elements[elem*numCorners+func]; for(c = 0; c < comp; c++, index++) { #ifdef PETSC_USE_BOPT_g if ((index < 0) || (index >= elemSize)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Bad index %d should be in [0,%d)", index, elemSize); } #endif newArrays[c][node] = array[index]; } } } /* Cleanup */ ierr = ElementVecDestroy(elemVec); CHKERRQ(ierr); ierr = VecDestroy(ghostVec); CHKERRQ(ierr); ierr = VecScatterDestroy(ghostScatter); CHKERRQ(ierr); #else ierr = VecGetArray(gvec, &array); CHKERRQ(ierr); for(node = 0; node < numNodes; node++) { nclass = classes[node]; if (localStart[field] == PETSC_NULL) { continue; } else if (localStart[field][nclass] == -1) { for(c = 0; c < comp; c++) { newArrays[c][node] = 0.0; } } else { for(c = 0; c < comp; c++) { index = offsets[node]+localStart[field][nclass]+c; #ifdef PETSC_USE_BOPT_g if ((index < 0) || (index >= size)) { SETERRQ4(PETSC_ERR_ARG_OUTOFRANGE, "Bad index %d+%d= %d should be in [0,%d)", offsets[node], localStart[field][nclass]+c, index, size); } #endif newArrays[c][node] = array[offsets[node]+localStart[field][nclass]+c]; } } } ierr = VecRestoreArray(gvec, &array); CHKERRQ(ierr); #endif /* Put the vector into the archive */ for(c = 0; c < comp; c++) { optlist = DBMakeOptlist(3); ierr = DBAddOption(optlist, DBOPT_LABEL, name); CHKERRQ(ierr); /* ierr = DBAddOption(optlist, DBOPT_UNITS, fieldUnits); CHKERRQ(ierr);*/ ierr = DBAddOption(optlist, DBOPT_UNITS, "units"); CHKERRQ(ierr); if (vsilo->meshName == PETSC_NULL) { ierr = DBPutUcdvar1(fp, subNames[c], "PetscMesh", newArrays[c], numNodes, PETSC_NULL, 0, DB_FLOAT, DB_NODECENT, optlist); CHKERRQ(ierr); } else { ierr = DBPutUcdvar1(fp, subNames[c], vsilo->meshName, newArrays[c], numNodes, PETSC_NULL, 0, DB_FLOAT, DB_NODECENT, optlist); CHKERRQ(ierr); } ierr = DBFreeOptlist(optlist); CHKERRQ(ierr); } /* Define a vecor variable for multicomponent fields */ if(comp > 1) { ierr = PetscMemzero(buf, 1024 * sizeof(char)); CHKERRQ(ierr); len = DBGetVarLength(fp, "_meshtv_defvars"); if (len > 0) { if (DBGetVarType(fp, "_meshtv_defvars") != DB_CHAR) { SETERRQ(PETSC_ERR_FILE_READ, "Invalid type for variable _meshtv_defvars"); } if (len > 1024) SETERRQ(PETSC_ERR_SUP, "Need to do dyanmic allocation here"); ierr = DBReadVar(fp, "_meshtv_defvars", buf); CHKERRQ(ierr); PetscStrcat(buf, ";vec"); CHKERRQ(ierr); } else { PetscStrcpy(buf, "vec"); CHKERRQ(ierr); } PetscStrcat(buf, name); CHKERRQ(ierr); PetscStrcat(buf, " vector {"); CHKERRQ(ierr); for(c = 0; c < comp-1; c++) { PetscStrcat(buf, subNames[c]); CHKERRQ(ierr); PetscStrcat(buf, ","); CHKERRQ(ierr); } PetscStrcat(buf, subNames[c]); CHKERRQ(ierr); PetscStrcat(buf, "}"); CHKERRQ(ierr); ierr = PetscStrlen(buf, &len); CHKERRQ(ierr); if (len > 1024) SETERRQ(PETSC_ERR_SUP, "Need to do dyanmic allocation here"); len = 1024; ierr = DBWrite(fp, "_meshtv_defvars", buf, &len, 1, DB_CHAR); CHKERRQ(ierr); } /* Cleanup */ for(c = 0; c < comp; c++) { ierr = PetscFree(newArrays[c]); CHKERRQ(ierr); } ierr = PetscFree(newArrays); CHKERRQ(ierr); } PetscFunctionReturn(0); } #endif #undef __FUNCT__ #define __FUNCT__ "GVecCreateElementVec_Private" static int GVecCreateElementVec_Private(GVec gvec, LocalVarOrdering locOrder, ElementVec *elemVec) { MPI_Comm comm; Grid grid; VarOrdering order; FieldClassMap cm; Discretization disc; PetscTruth isConstrained, isFieldConstrained; int numFields, numFuncs, elemSize; int f, field; int ierr; PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject) gvec, &comm); CHKERRQ(ierr); ierr = GVecGetGrid(gvec, &grid); CHKERRQ(ierr); ierr = GVecGetOrder(gvec, &order); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(order, &cm); CHKERRQ(ierr); ierr = LocalVarOrderingGetSize(locOrder, &elemSize); CHKERRQ(ierr); ierr = FieldClassMapIsConstrained(cm, &isConstrained); CHKERRQ(ierr); ierr = FieldClassMapGetNumFields(cm, &numFields); CHKERRQ(ierr); if (isConstrained == PETSC_TRUE) { for(f = 0; f < numFields; f++) { ierr = FieldClassMapGetField(cm, f, &field); CHKERRQ(ierr); ierr = GridIsFieldConstrained(grid, field, &isFieldConstrained); CHKERRQ(ierr); if (isFieldConstrained == PETSC_TRUE) { ierr = GridGetDiscretization(grid, field, &disc); CHKERRQ(ierr); ierr = DiscretizationGetNumFunctions(disc, &numFuncs); CHKERRQ(ierr); elemSize += numFuncs*grid->fields[field].constraintCompDiff; } } } ierr = ElementVecCreate(comm, elemSize, elemVec); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNumLocalVertices_Private" static int MeshGetNumLocalVertices_Private(Mesh mesh, int *numLocalVertices) { Partition part; PetscTruth *seen; int numOverlapElements, numCorners, numLocalNodes; int elem, corner, node; int ierr; PetscFunctionBegin; *numLocalVertices = 0; ierr = MeshGetPartition(mesh, &part); CHKERRQ(ierr); ierr = MeshGetNumCorners(mesh, &numCorners); CHKERRQ(ierr); ierr = PartitionGetNumOverlapElements(part, &numOverlapElements); CHKERRQ(ierr); ierr = PartitionGetNumNodes(part, &numLocalNodes); CHKERRQ(ierr); if (numCorners == 3) { *numLocalVertices = numLocalNodes; PetscFunctionReturn(0); } ierr = PetscMalloc(numLocalNodes * sizeof(PetscTruth), &seen); CHKERRQ(ierr); ierr = PetscMemzero(seen, numLocalNodes * sizeof(PetscTruth)); CHKERRQ(ierr); for(elem = 0; elem < numOverlapElements; elem++) { for(corner = 0; corner < 3; corner++) { ierr = MeshGetNodeFromElement(mesh, elem, corner, &node); CHKERRQ(ierr); if ((node < numLocalNodes) && (seen[node] == PETSC_FALSE)) { seen[node] = PETSC_TRUE; *numLocalVertices += 1; } } } ierr = PetscFree(seen); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNumVertices_Private" static int MeshGetNumVertices_Private(Mesh mesh, int *numVertices) { int ierr; PetscFunctionBegin; ierr = MeshGetInfo(mesh, numVertices, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecCanonicalizeField_Private" static int GVecCanonicalizeField_Private(GVec gvec, int canonicalField, LocalVarOrdering locOrder, Vec *canonicalVec) { MPI_Comm comm; Grid grid; Mesh mesh; Partition part; VarOrdering order, reductionOrder; AO meshOrdering, partOrdering; Discretization disc; Vec ghostVec, reductionVec; VecScatter ghostScatter; ElementVec elemVec; PetscTruth reduceSystem; PetscScalar reduceAlpha; PetscScalar minusOne = -1.0; int *nodes, *indices; PetscScalar *array; int numLocNodes, numNodes, fieldStart, numComp, numLocElements, numCorners, numFuncs, elemSize; int elem, func, node, c; int ierr; PetscFunctionBegin; PetscValidPointer(canonicalVec); ierr = PetscObjectGetComm((PetscObject) gvec, &comm); CHKERRQ(ierr); ierr = GVecGetGrid(gvec, &grid); CHKERRQ(ierr); ierr = GVecGetOrder(gvec, &order); CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) gvec, "reductionOrder", (PetscObject *) &reductionOrder); CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) gvec, "reductionVec", (PetscObject *) &reductionVec); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetPartition(mesh, &part); CHKERRQ(ierr); if (reductionVec == PETSC_NULL) reductionVec = grid->bdReduceVecCur; /* Create canonical vector */ ierr = MeshGetNumCorners(mesh, &numCorners); CHKERRQ(ierr); ierr = GridGetDiscretization(grid, canonicalField, &disc); CHKERRQ(ierr); ierr = DiscretizationGetNumFunctions(disc, &numFuncs); CHKERRQ(ierr); ierr = GridGetFieldComponents(grid, canonicalField, &numComp); CHKERRQ(ierr); if (numFuncs < numCorners) { ierr = MeshGetNumLocalVertices_Private(mesh, &numLocNodes); CHKERRQ(ierr); ierr = MeshGetNumVertices_Private(mesh, &numNodes); CHKERRQ(ierr); } else { ierr = PartitionGetNumNodes(part, &numLocNodes); CHKERRQ(ierr); ierr = PartitionGetTotalNodes(part, &numNodes); CHKERRQ(ierr); } ierr = VecCreate(comm, canonicalVec); CHKERRQ(ierr); ierr = VecSetSizes(*canonicalVec, numLocNodes*numComp, numNodes*numComp); CHKERRQ(ierr); ierr = VecSetType(*canonicalVec, VECMPI); CHKERRQ(ierr); /* Setup reduction */ ierr = (*grid->ops->gridsetupghostscatter)(grid, order, &ghostVec, &ghostScatter); CHKERRQ(ierr); /* Fill the local vector */ ierr = GridGlobalToLocalGeneral(grid, gvec, ghostVec, INSERT_VALUES, ghostScatter); CHKERRQ(ierr); /* Setup element vector */ ierr = GVecCreateElementVec_Private(gvec, locOrder, &elemVec); CHKERRQ(ierr); ierr = LocalVarOrderingGetFieldStart(locOrder, canonicalField, &fieldStart); CHKERRQ(ierr); array = &elemVec->array[fieldStart]; ierr = GridGetReduceSystem(grid, &reduceSystem); CHKERRQ(ierr); ierr = GridGetBCMultiplier(grid, &reduceAlpha); CHKERRQ(ierr); ierr = GridSetBCMultiplier(grid, minusOne); CHKERRQ(ierr); ierr = MeshGetNodeOrdering(mesh, &meshOrdering); CHKERRQ(ierr); ierr = PartitionGetNodeOrdering(part, &partOrdering); CHKERRQ(ierr); ierr = LocalVarOrderingGetSize(locOrder, &elemSize); CHKERRQ(ierr); ierr = PetscMalloc(numFuncs * sizeof(int), &nodes); CHKERRQ(ierr); ierr = PetscMalloc(numFuncs*numComp * sizeof(int), &indices); CHKERRQ(ierr); ierr = PartitionGetNumElements(part, &numLocElements); CHKERRQ(ierr); for(elem = 0; elem < numLocElements; elem++) { /* Initialize element vector */ ierr = ElementVecZero(elemVec); CHKERRQ(ierr); elemVec->reduceSize = elemSize; /* Setup local row indices */ ierr = GridCalcGeneralElementVecIndices(grid, elem, order, reductionOrder, PETSC_TRUE, elemVec); CHKERRQ(ierr); /* Setup local vectors */ ierr = GridLocalToElementGeneral(grid, ghostVec, reductionVec, reduceSystem, PETSC_TRUE, elemVec); CHKERRQ(ierr); /* Must transform to unconstrained variables */ ierr = GridProjectElementVec(grid, mesh, elem, order, PETSC_FALSE, elemVec); CHKERRQ(ierr); /* Put values into canonical vector */ for(func = 0; func < numFuncs; func++) { ierr = MeshGetNodeFromElement(mesh, elem, func, &node); CHKERRQ(ierr); ierr = PartitionLocalToGlobalNodeIndex(part, node, &nodes[func]); CHKERRQ(ierr); } /* We must globally renumber so that midnodes come after vertices */ if (partOrdering != PETSC_NULL) { ierr = AOPetscToApplication(partOrdering, numFuncs, nodes); CHKERRQ(ierr); } if (meshOrdering != PETSC_NULL) { ierr = AOPetscToApplication(meshOrdering, numFuncs, nodes); CHKERRQ(ierr); } for(func = 0; func < numFuncs; func++) { for(c = 0; c < numComp; c++) { indices[func*numComp+c] = nodes[func]*numComp+c; } } ierr = VecSetValues(*canonicalVec, numFuncs*numComp, indices, array, INSERT_VALUES); CHKERRQ(ierr); /* Usually reset by ElementVecSetValues() */ elemVec->reduceSize = elemVec->size; } ierr = VecAssemblyBegin(*canonicalVec); CHKERRQ(ierr); ierr = VecAssemblyEnd(*canonicalVec); CHKERRQ(ierr); ierr = GridSetBCMultiplier(grid, reduceAlpha); CHKERRQ(ierr); /* Cleanup */ ierr = PetscFree(nodes); CHKERRQ(ierr); ierr = PetscFree(indices); CHKERRQ(ierr); ierr = ElementVecDestroy(elemVec); CHKERRQ(ierr); ierr = VecDestroy(ghostVec); CHKERRQ(ierr); ierr = VecScatterDestroy(ghostScatter); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecView_Triangular_2D_VU" int GVecView_Triangular_2D_VU(GVec gvec, PetscViewer viewer) { MPI_Comm comm; Grid grid; Mesh mesh; Partition part; VarOrdering order; LocalVarOrdering locOrder; FieldClassMap cm; Discretization disc; Vec locCanonicalVec, canonicalVec; PetscTruth vecSeen; FILE *fp; PetscScalar *values; char *fieldName; int *fields; int rank, numFields, numComp, numFuncs, numNodes; int f, field, node, c; int ierr; PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject) gvec, &comm); CHKERRQ(ierr); ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); ierr = GVecGetGrid(gvec, &grid); CHKERRQ(ierr); ierr = GVecGetOrder(gvec, &order); CHKERRQ(ierr); ierr = GridGetMesh(grid, &mesh); CHKERRQ(ierr); ierr = MeshGetPartition(mesh, &part); CHKERRQ(ierr); ierr = VarOrderingGetClassMap(order, &cm); CHKERRQ(ierr); ierr = FieldClassMapGetNumFields(cm, &numFields); CHKERRQ(ierr); ierr = PetscViewerVUGetPointer(viewer, &fp); CHKERRQ(ierr); /* Create local ordering */ ierr = PetscMalloc(numFields * sizeof(int), &fields); CHKERRQ(ierr); for(f = 0; f < numFields; f++) { ierr = FieldClassMapGetField(cm, f, &fields[f]); CHKERRQ(ierr); } ierr = LocalVarOrderingCreate(grid, numFields, fields, &locOrder); CHKERRQ(ierr); ierr = PetscFree(fields); CHKERRQ(ierr); /* Write each field */ PetscFPrintf(comm, fp, "// Field values\n"); for(f = 0; f < numFields; f++) { ierr = FieldClassMapGetField(cm, f, &field); CHKERRQ(ierr); ierr = GridGetFieldComponents(grid, field, &numComp); CHKERRQ(ierr); ierr = GridGetFieldName(grid, field, &fieldName); CHKERRQ(ierr); if (fieldName == PETSC_NULL) SETERRQ(PETSC_ERR_ARG_WRONG, "You must name the fields for VU"); ierr = GridGetDiscretization(grid, field, &disc); CHKERRQ(ierr); ierr = DiscretizationGetNumFunctions(disc, &numFuncs); CHKERRQ(ierr); ierr = GVecCanonicalizeField_Private(gvec, field, locOrder, &canonicalVec); CHKERRQ(ierr); ierr = VecConvertMPIToMPIZero(canonicalVec, &locCanonicalVec); CHKERRQ(ierr); ierr = VecGetSize(locCanonicalVec, &numNodes); CHKERRQ(ierr); ierr = VecGetArray(locCanonicalVec, &values); CHKERRQ(ierr); if (rank == 0) { if (numNodes%numComp) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid size for canonical field vector"); numNodes /= numComp; for(c = 0; c < numComp; c++) { PetscFPrintf(comm, fp, "FIELD %s%d( ) = {\n", fieldName, c); for(node = 0; node < numNodes; node++) { PetscFPrintf(comm, fp, "%g\n", PetscRealPart(values[node*numComp+c])); } PetscFPrintf(comm, fp, "};\n\n"); } } ierr = VecRestoreArray(locCanonicalVec, &values); CHKERRQ(ierr); ierr = VecDestroy(locCanonicalVec); CHKERRQ(ierr); ierr = VecDestroy(canonicalVec); CHKERRQ(ierr); ierr = PetscFree(fieldName); CHKERRQ(ierr); } ierr = LocalVarOrderingDestroy(locOrder); CHKERRQ(ierr); /* Write solution as field union */ ierr = PetscViewerVUGetVecSeen(viewer, &vecSeen); CHKERRQ(ierr); if (vecSeen == PETSC_FALSE) { ierr = PetscViewerVUPrintDeferred(viewer, "// The full solution\nSOLUTION VecSolution( ) =\n{\n"); CHKERRQ(ierr); ierr = PetscViewerVUSetVecSeen(viewer, PETSC_TRUE); CHKERRQ(ierr); } for(f = 0; f < numFields; f++) { ierr = FieldClassMapGetField(cm, f, &field); CHKERRQ(ierr); ierr = GridGetFieldName(grid, field, &fieldName); CHKERRQ(ierr); ierr = GridGetFieldComponents(grid, field, &numComp); CHKERRQ(ierr); ierr = GridGetDiscretization(grid, field, &disc); CHKERRQ(ierr); ierr = DiscretizationGetNumFunctions(disc, &numFuncs); CHKERRQ(ierr); for(c = 0; c < numComp; c++) { switch(numFuncs) { case 3: ierr = PetscViewerVUPrintDeferred(viewer, " VARIABLE %s%d(LagrTrian03, %s%d, Connectivity, Zone1);\n", fieldName, c, fieldName, c);CHKERRQ(ierr); break; case 6: ierr = PetscViewerVUPrintDeferred(viewer, " VARIABLE %s%d(LagrTrian06, %s%d, FullConnectivity, Zone1);\n", fieldName, c, fieldName, c);CHKERRQ(ierr); break; default: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of function in discretization %d", numFuncs); } } ierr = PetscFree(fieldName); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GVecView_Triangular_2D" int GVecView_Triangular_2D(GVec gvec, PetscViewer viewer) { Grid grid; PetscTruth isascii, isdraw, isvu, ismathematica, issilo; int ierr; PetscFunctionBegin; ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii); CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_DRAW, &isdraw); CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_VU, &isvu); CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_MATHEMATICA, &ismathematica); CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_SILO, &issilo); CHKERRQ(ierr); if (isascii == PETSC_TRUE) { ierr = GVecGetGrid(gvec, &grid); CHKERRQ(ierr); ierr = GridView(grid, viewer); CHKERRQ(ierr); ierr = PetscViewerFlush(viewer); CHKERRQ(ierr); ierr = VecView(gvec, viewer); CHKERRQ(ierr); } else if (isdraw == PETSC_TRUE) { ierr = GVecView_Triangular_2D_Draw(gvec, viewer); CHKERRQ(ierr); } else if (isvu == PETSC_TRUE) { ierr = GVecView_Triangular_2D_VU(gvec, viewer); CHKERRQ(ierr); #ifdef PETSC_HAVE_MATHEMATICA } else if (ismathematica == PETSC_TRUE) { ierr = ViewerMathematica_GVec_Triangular_2D(viewer, gvec); CHKERRQ(ierr); #endif #ifdef PETSC_HAVE_MATHEMATICA } else if (issilo == PETSC_TRUE) { ierr = GVecView_Triangular_2D_Silo(gvec, viewer); CHKERRQ(ierr); #endif } PetscFunctionReturn(0); }