#include "private/meshimpl.h" /*I "petscmesh.h" I*/ #include #include #include "src/dm/mesh/meshvtk.h" #include "src/dm/mesh/meshpcice.h" #include "src/dm/mesh/meshpylith.h" /* Logging support */ PetscCookie PETSCDM_DLLEXPORT MESH_COOKIE = 0; PetscEvent Mesh_View = 0, Mesh_GetGlobalScatter = 0, Mesh_restrictVector = 0, Mesh_assembleVector = 0, Mesh_assembleVectorComplete = 0, Mesh_assembleMatrix = 0, Mesh_updateOperator = 0; PetscTruth MeshRegisterAllCalled = PETSC_FALSE; PetscFList MeshList; EXTERN PetscErrorCode MeshView_Mesh(Mesh, PetscViewer); EXTERN PetscErrorCode MeshRefine_Mesh(Mesh, MPI_Comm, Mesh *); EXTERN PetscErrorCode MeshCoarsenHierarchy_Mesh(Mesh, int, Mesh **); EXTERN PetscErrorCode MeshGetInterpolation_Mesh(Mesh, Mesh, Mat *, Vec *); EXTERN PetscErrorCode updateOperatorCompat(Mat, const ALE::Obj&, const ALE::Obj&, const ALECompat::Mesh::point_type&, PetscScalar[], InsertMode); EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "Mesh_DelTag" /* Private routine to delete internal tag storage when a communicator is freed. This is called by MPI, not by users. Note: this is declared extern "C" because it is passed to MPI_Keyval_create we do not use PetscFree() since it is unsafe after PetscFinalize() */ PetscMPIInt PETSC_DLLEXPORT Mesh_DelTag(MPI_Comm comm,PetscMPIInt keyval,void* attr_val,void* extra_state) { free(attr_val); return(MPI_SUCCESS); } EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "MeshFinalize" PetscErrorCode MeshFinalize() { PetscFunctionBegin; ALE::Mesh::NumberingFactory::singleton(0, 0, true); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshView_Sieve_Ascii" PetscErrorCode MeshView_Sieve_Ascii(const ALE::Obj& mesh, PetscViewer viewer) { PetscViewerFormat format; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); if (format == PETSC_VIEWER_ASCII_VTK) { ierr = VTKViewer::writeHeader(viewer);CHKERRQ(ierr); ierr = VTKViewer::writeVertices(mesh, viewer);CHKERRQ(ierr); ierr = VTKViewer::writeElements(mesh, viewer);CHKERRQ(ierr); } else if (format == PETSC_VIEWER_ASCII_PYLITH) { char *filename; char connectFilename[2048]; char coordFilename[2048]; ierr = PetscViewerFileGetName(viewer, &filename);CHKERRQ(ierr); ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);CHKERRQ(ierr); ierr = PetscStrcpy(connectFilename, filename);CHKERRQ(ierr); ierr = PetscStrcat(connectFilename, ".connect");CHKERRQ(ierr); ierr = PetscViewerFileSetName(viewer, connectFilename);CHKERRQ(ierr); ierr = ALE::PyLith::Viewer::writeElements(mesh, mesh->getIntSection("material"), viewer);CHKERRQ(ierr); ierr = PetscStrcpy(coordFilename, filename);CHKERRQ(ierr); ierr = PetscStrcat(coordFilename, ".coord");CHKERRQ(ierr); ierr = PetscViewerFileSetName(viewer, coordFilename);CHKERRQ(ierr); ierr = ALE::PyLith::Viewer::writeVertices(mesh, viewer);CHKERRQ(ierr); ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); ierr = PetscExceptionTry1(PetscViewerFileSetName(viewer, filename), PETSC_ERR_FILE_OPEN); if (PetscExceptionValue(ierr)) { /* this means that a caller above me has also tryed this exception so I don't handle it here, pass it up */ } else if (PetscExceptionCaught(ierr, PETSC_ERR_FILE_OPEN)) { ierr = 0; } CHKERRQ(ierr); } else if (format == PETSC_VIEWER_ASCII_PYLITH_LOCAL) { PetscViewer connectViewer, coordViewer; char *filename; char localFilename[2048]; int rank = mesh->commRank(); ierr = PetscViewerFileGetName(viewer, &filename);CHKERRQ(ierr); sprintf(localFilename, "%s.%d.connect", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &connectViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(connectViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(connectViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(connectViewer, localFilename);CHKERRQ(ierr); ierr = ALE::PyLith::Viewer::writeElementsLocal(mesh, mesh->getIntSection("material"), connectViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(connectViewer);CHKERRQ(ierr); sprintf(localFilename, "%s.%d.coord", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &coordViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(coordViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(coordViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(coordViewer, localFilename);CHKERRQ(ierr); ierr = ALE::PyLith::Viewer::writeVerticesLocal(mesh, coordViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(coordViewer);CHKERRQ(ierr); } else if (format == PETSC_VIEWER_ASCII_PCICE) { char *filename; char coordFilename[2048]; PetscTruth isConnect; size_t len; ierr = PetscViewerFileGetName(viewer, &filename);CHKERRQ(ierr); ierr = PetscStrlen(filename, &len);CHKERRQ(ierr); ierr = PetscStrcmp(&(filename[len-5]), ".lcon", &isConnect);CHKERRQ(ierr); if (!isConnect) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid element connectivity filename: %s", filename); } ierr = ALE::PCICE::Viewer::writeElements(mesh, viewer);CHKERRQ(ierr); ierr = PetscStrncpy(coordFilename, filename, len-5);CHKERRQ(ierr); coordFilename[len-5] = '\0'; ierr = PetscStrcat(coordFilename, ".nodes");CHKERRQ(ierr); ierr = PetscViewerFileSetName(viewer, coordFilename);CHKERRQ(ierr); ierr = ALE::PCICE::Viewer::writeVertices(mesh, viewer);CHKERRQ(ierr); } else { int dim = mesh->getDimension(); ierr = PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);CHKERRQ(ierr); for(int d = 0; d <= dim; d++) { // FIX: Need to globalize ierr = PetscViewerASCIIPrintf(viewer, " %d %d-cells\n", mesh->depthStratum(d)->size(), d);CHKERRQ(ierr); } } ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCompatView_Sieve_Ascii" PetscErrorCode MeshCompatView_Sieve_Ascii(const ALE::Obj& mesh, PetscViewer viewer) { PetscViewerFormat format; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); if (format == PETSC_VIEWER_ASCII_PYLITH) { char *filename; char connectFilename[2048]; char coordFilename[2048]; ierr = PetscViewerFileGetName(viewer, &filename);CHKERRQ(ierr); ierr = PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);CHKERRQ(ierr); ierr = PetscStrcpy(connectFilename, filename);CHKERRQ(ierr); ierr = PetscStrcat(connectFilename, ".connect");CHKERRQ(ierr); ierr = PetscViewerFileSetName(viewer, connectFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeElements(mesh, mesh->getIntSection("material"), viewer);CHKERRQ(ierr); ierr = PetscStrcpy(coordFilename, filename);CHKERRQ(ierr); ierr = PetscStrcat(coordFilename, ".coord");CHKERRQ(ierr); ierr = PetscViewerFileSetName(viewer, coordFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeVertices(mesh, viewer);CHKERRQ(ierr); ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); ierr = PetscExceptionTry1(PetscViewerFileSetName(viewer, filename), PETSC_ERR_FILE_OPEN); if (PetscExceptionValue(ierr)) { /* this means that a caller above me has also tryed this exception so I don't handle it here, pass it up */ } else if (PetscExceptionCaught(ierr, PETSC_ERR_FILE_OPEN)) { ierr = 0; } CHKERRQ(ierr); } else if (format == PETSC_VIEWER_ASCII_PYLITH_LOCAL) { PetscViewer connectViewer, coordViewer; char *filename; char localFilename[2048]; int rank = mesh->commRank(); ierr = PetscViewerFileGetName(viewer, &filename);CHKERRQ(ierr); sprintf(localFilename, "%s.%d.connect", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &connectViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(connectViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(connectViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(connectViewer, localFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeElementsLocal(mesh, mesh->getIntSection("material"), connectViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(connectViewer);CHKERRQ(ierr); sprintf(localFilename, "%s.%d.coord", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &coordViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(coordViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(coordViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(coordViewer, localFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeVerticesLocal(mesh, coordViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(coordViewer);CHKERRQ(ierr); if (mesh->hasPairSection("split")) { PetscViewer splitViewer; sprintf(localFilename, "%s.%d.split", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &splitViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(splitViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(splitViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(splitViewer, localFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeSplitLocal(mesh, mesh->getPairSection("split"), splitViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(splitViewer);CHKERRQ(ierr); } if (mesh->hasRealSection("traction")) { PetscViewer tractionViewer; sprintf(localFilename, "%s.%d.traction", filename, rank); ierr = PetscViewerCreate(PETSC_COMM_SELF, &tractionViewer);CHKERRQ(ierr); ierr = PetscViewerSetType(tractionViewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr); ierr = PetscViewerSetFormat(tractionViewer, PETSC_VIEWER_ASCII_PYLITH);CHKERRQ(ierr); ierr = PetscViewerFileSetName(tractionViewer, localFilename);CHKERRQ(ierr); ierr = ALECompat::PyLith::Viewer::writeTractionsLocal(mesh, mesh->getRealSection("traction"), tractionViewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(tractionViewer);CHKERRQ(ierr); } } else { int dim = mesh->getDimension(); ierr = PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);CHKERRQ(ierr); for(int d = 0; d <= dim; d++) { // FIX: Need to globalize ierr = PetscViewerASCIIPrintf(viewer, " %d %d-cells\n", mesh->getTopology()->depthStratum(0, d)->size(), d);CHKERRQ(ierr); } } ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshView_Sieve" PetscErrorCode MeshView_Sieve(const ALE::Obj& mesh, PetscViewer viewer) { PetscTruth iascii, isbinary, isdraw; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &isbinary);CHKERRQ(ierr); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_DRAW, &isdraw);CHKERRQ(ierr); if (iascii){ ierr = MeshView_Sieve_Ascii(mesh, viewer);CHKERRQ(ierr); } else if (isbinary) { SETERRQ(PETSC_ERR_SUP, "Binary viewer not implemented for Mesh"); } else if (isdraw){ SETERRQ(PETSC_ERR_SUP, "Draw viewer not implemented for Mesh"); } else { SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshView_Mesh" PetscErrorCode PETSCDM_DLLEXPORT MeshView_Mesh(Mesh mesh, PetscViewer viewer) { PetscErrorCode ierr; PetscFunctionBegin; if (!mesh->mcompat.isNull()) { ierr = MeshCompatView_Sieve_Ascii(mesh->mcompat, viewer);CHKERRQ(ierr); } else { ierr = MeshView_Sieve(mesh->m, viewer);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshView" /*@C MeshView - Views a Mesh object. Collective on Mesh Input Parameters: + mesh - the mesh - viewer - an optional visualization context Notes: The available visualization contexts include + PETSC_VIEWER_STDOUT_SELF - standard output (default) - PETSC_VIEWER_STDOUT_WORLD - synchronized standard output where only the first processor opens the file. All other processors send their data to the first processor to print. You can change the format the mesh is printed using the option PetscViewerSetFormat(). The user can open alternative visualization contexts with + PetscViewerASCIIOpen() - Outputs mesh to a specified file . PetscViewerBinaryOpen() - Outputs mesh in binary to a specified file; corresponding input uses MeshLoad() . PetscViewerDrawOpen() - Outputs mesh to an X window display The user can call PetscViewerSetFormat() to specify the output format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include + PETSC_VIEWER_ASCII_DEFAULT - default, prints mesh information - PETSC_VIEWER_ASCII_VTK - outputs a VTK file describing the mesh Level: beginner Concepts: mesh^printing Concepts: mesh^saving to disk .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscViewerBinaryOpen(), MeshLoad(), PetscViewerCreate() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshView(Mesh mesh, PetscViewer viewer) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); PetscValidType(mesh, 1); if (!viewer) { ierr = PetscViewerASCIIGetStdout(mesh->comm,&viewer);CHKERRQ(ierr); } PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE, 2); PetscCheckSameComm(mesh, 1, viewer, 2); ierr = PetscLogEventBegin(Mesh_View,0,0,0,0);CHKERRQ(ierr); ierr = (*mesh->ops->view)(mesh, viewer);CHKERRQ(ierr); ierr = PetscLogEventEnd(Mesh_View,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshLoad" /*@C MeshLoad - Create a mesh topology from the saved data in a viewer. Collective on Viewer Input Parameter: . viewer - The viewer containing the data Output Parameters: . mesh - the mesh object Level: advanced .seealso MeshView() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshLoad(PetscViewer viewer, Mesh *mesh) { SETERRQ(PETSC_ERR_SUP, ""); } #undef __FUNCT__ #define __FUNCT__ "MeshGetMesh" /*@C MeshGetMesh - Gets the internal mesh object Not collective Input Parameter: . mesh - the mesh object Output Parameter: . m - the internal mesh object Level: advanced .seealso MeshCreate(), MeshSetMesh() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshGetMesh(Mesh mesh, ALE::Obj& m) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); m = mesh->m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetMesh" /*@C MeshSetMesh - Sets the internal mesh object Not collective Input Parameters: + mesh - the mesh object - m - the internal mesh object Level: advanced .seealso MeshCreate(), MeshGetMesh() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshSetMesh(Mesh mesh, const ALE::Obj& m) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); mesh->m = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreateMatrix" template PetscErrorCode PETSCDM_DLLEXPORT MeshCreateMatrix(const Obj& mesh, const Obj
& section, MatType mtype, Mat *J) { const ALE::Obj& order = mesh->getFactory()->getGlobalOrder(mesh, "default", section); int localSize = order->getLocalSize(); int globalSize = order->getGlobalSize(); PetscTruth isShell, isBlock, isSeqBlock, isMPIBlock; PetscErrorCode ierr; PetscFunctionBegin; ierr = MatCreate(mesh->comm(), J);CHKERRQ(ierr); ierr = MatSetSizes(*J, localSize, localSize, globalSize, globalSize);CHKERRQ(ierr); ierr = MatSetType(*J, mtype);CHKERRQ(ierr); ierr = MatSetFromOptions(*J);CHKERRQ(ierr); ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr); ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr); ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr); ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr); if (!isShell) { int bs = 1; if (isBlock || isSeqBlock || isMPIBlock) { bs = section->getFiberDimension(*section->getChart().begin()); } ierr = preallocateOperator(mesh, bs, section->getAtlas(), order, *J);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreateMatrix" PetscErrorCode PETSCDM_DLLEXPORT MeshCreateMatrix(Mesh mesh, SectionReal section, MatType mtype, Mat *J) { ALE::Obj m; ALE::Obj s; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionRealGetSection(section, s);CHKERRQ(ierr); ierr = MeshCreateMatrix(m, s, mtype, J);CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *J, "mesh", (PetscObject) mesh);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetVertexMatrix" PetscErrorCode PETSCDM_DLLEXPORT MeshGetVertexMatrix(Mesh mesh, MatType mtype, Mat *J) { SectionReal section; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetVertexSectionReal(mesh, 1, §ion);CHKERRQ(ierr); ierr = MeshCreateMatrix(mesh, section, mtype, J);CHKERRQ(ierr); ierr = SectionRealDestroy(section);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetMatrix" /*@C MeshGetMatrix - Creates a matrix with the correct parallel layout required for computing the Jacobian on a function defined using the informatin in Mesh. Collective on Mesh Input Parameter: + mesh - the mesh object - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.). Output Parameters: . J - matrix with the correct nonzero preallocation (obviously without the correct Jacobian values) Level: advanced Notes: This properly preallocates the number of nonzeros in the sparse matrix so you do not need to do it yourself. .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshGetMatrix(Mesh mesh, MatType mtype, Mat *J) { ALE::Obj m; PetscTruth flag; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshHasSectionReal(mesh, "default", &flag);CHKERRQ(ierr); if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section"); ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = MeshCreateMatrix(m, m->getRealSection("default"), mtype, J);CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *J, "mesh", (PetscObject) mesh);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreate" /*@C MeshCreate - Creates a DM object, used to manage data for an unstructured problem described by a Sieve. Collective on MPI_Comm Input Parameter: . comm - the processors that will share the global vector Output Parameters: . mesh - the mesh object Level: advanced .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshGetGlobalIndices() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshCreate(MPI_Comm comm,Mesh *mesh) { PetscErrorCode ierr; Mesh p; PetscFunctionBegin; PetscValidPointer(mesh,2); *mesh = PETSC_NULL; #ifndef PETSC_USE_DYNAMIC_LIBRARIES ierr = DMInitializePackage(PETSC_NULL);CHKERRQ(ierr); #endif ierr = PetscHeaderCreate(p,_p_Mesh,struct _MeshOps,MESH_COOKIE,0,"Mesh",comm,MeshDestroy,MeshView);CHKERRQ(ierr); p->ops->view = MeshView_Mesh; p->ops->destroy = PETSC_NULL; p->ops->createglobalvector = MeshCreateGlobalVector; p->ops->getcoloring = PETSC_NULL; p->ops->getmatrix = MeshGetMatrix; p->ops->getinterpolation = MeshGetInterpolation_Mesh; p->ops->getinjection = PETSC_NULL; p->ops->refine = MeshRefine_Mesh; p->ops->coarsen = PETSC_NULL; p->ops->refinehierarchy = PETSC_NULL; p->ops->coarsenhierarchy = MeshCoarsenHierarchy_Mesh; ierr = PetscObjectChangeTypeName((PetscObject) p, "sieve");CHKERRQ(ierr); p->m = PETSC_NULL; p->globalScatter = PETSC_NULL; p->lf = PETSC_NULL; p->lj = PETSC_NULL; p->mcompat = PETSC_NULL; p->data = PETSC_NULL; *mesh = p; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshDestroy" /*@ MeshDestroy - Destroys a mesh. Collective on Mesh Input Parameter: . mesh - the mesh object Level: advanced .seealso MeshCreate(), MeshCreateGlobalVector(), MeshGetGlobalIndices() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshDestroy(Mesh mesh) { PetscErrorCode ierr; PetscFunctionBegin; if (--mesh->refct > 0) PetscFunctionReturn(0); if (mesh->globalScatter) {ierr = VecScatterDestroy(mesh->globalScatter);CHKERRQ(ierr);} mesh->m = PETSC_NULL; ierr = PetscHeaderDestroy(mesh);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetType" /*@C MeshSetType - Sets the Mesh type Collective on Mesh Input Parameters: + mesh - the Mesh context - type - the type Options Database Key: . -mesh_type - Sets the type; use -help for a list of available types (for instance, cartesian or sieve) Notes: See "petsc/include/petscmesh.h" for available types (for instance, MESHCARTESIAN or MESHSIEVE). Level: intermediate .keywords: Mesh, set, typr .seealso: MeshGetType(), MeshType @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshSetType(Mesh mesh, MeshType type) { PetscErrorCode ierr,(*r)(Mesh); PetscTruth match; PetscFunctionBegin; PetscValidHeaderSpecific(mesh,MESH_COOKIE,1); PetscValidCharPointer(type,2); ierr = PetscTypeCompare((PetscObject)mesh,type,&match);CHKERRQ(ierr); if (match) PetscFunctionReturn(0); ierr = PetscFListFind(MeshList,mesh->comm,type,(void (**)(void)) &r);CHKERRQ(ierr); if (!r) SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested Mesh type %s",type); /* Destroy the previous private Mesh context */ if (mesh->ops->destroy) { ierr = (*mesh->ops->destroy)(mesh);CHKERRQ(ierr); } /* Reinitialize function pointers in MeshOps structure */ ierr = PetscMemzero(mesh->ops, sizeof(struct _MeshOps));CHKERRQ(ierr); /* Call the MeshCreate_XXX routine for this particular mesh */ ierr = (*r)(mesh);CHKERRQ(ierr); ierr = PetscObjectChangeTypeName((PetscObject) mesh, type);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetType" /*@C MeshGetType - Gets the Mesh type as a string from the Mesh object. Not Collective Input Parameter: . mesh - Mesh context Output Parameter: . name - name of Mesh type Level: intermediate .keywords: Mesh, get, type .seealso: MeshSetType() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshGetType(Mesh mesh,MeshType *type) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh,MESH_COOKIE,1); PetscValidPointer(type,2); *type = mesh->type_name; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRegister" /*@C MeshRegister - See MeshRegisterDynamic() Level: advanced @*/ PetscErrorCode PETSCKSP_DLLEXPORT MeshRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mesh)) { PetscErrorCode ierr; char fullname[PETSC_MAX_PATH_LEN]; PetscFunctionBegin; ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); ierr = PetscFListAdd(&MeshList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_BEGIN EXTERN PetscErrorCode PETSCKSP_DLLEXPORT MeshCreate_Cartesian(Mesh); EXTERN_C_END #undef __FUNCT__ #define __FUNCT__ "MeshRegisterAll" /*@C MeshRegisterAll - Registers all of the Mesh types in the Mesh package. Not Collective Level: advanced .keywords: Mesh, register, all .seealso: MeshRegisterDestroy() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshRegisterAll(const char path[]) { PetscErrorCode ierr; PetscFunctionBegin; MeshRegisterAllCalled = PETSC_TRUE; ierr = MeshRegisterDynamic(MESHCARTESIAN, path, "MeshCreate_Cartesian", MeshCreate_Cartesian);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRegisterDestroy" /*@ MeshRegisterDestroy - Frees the list of Mesh types that were registered by MeshRegister(). Not Collective Level: advanced .keywords: Mesh, register, destroy .seealso: MeshRegister(), MeshRegisterAll() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshRegisterDestroy(void) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFListDestroy(&MeshList);CHKERRQ(ierr); MeshRegisterAllCalled = PETSC_FALSE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreateGlobalVector" /*@C MeshCreateGlobalVector - Creates a vector of the correct size to be gathered into by the mesh. Collective on Mesh Input Parameter: . mesh - the mesh object Output Parameters: . gvec - the global vector Level: advanced Notes: Once this has been created you cannot add additional arrays or vectors to be packed. .seealso MeshDestroy(), MeshCreate(), MeshGetGlobalIndices() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshCreateGlobalVector(Mesh mesh, Vec *gvec) { ALE::Obj m; PetscTruth flag; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshHasSectionReal(mesh, "default", &flag);CHKERRQ(ierr); if (!flag) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Must set default section"); ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); const ALE::Obj& order = m->getFactory()->getGlobalOrder(m, "default", m->getRealSection("default")); ierr = VecCreate(m->comm(), gvec);CHKERRQ(ierr); ierr = VecSetSizes(*gvec, order->getLocalSize(), order->getGlobalSize());CHKERRQ(ierr); ierr = VecSetFromOptions(*gvec);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetGlobalIndices" /*@C MeshGetGlobalIndices - Gets the global indices for all the local entries Collective on Mesh Input Parameter: . mesh - the mesh object Output Parameters: . idx - the individual indices for each packed vector/array Level: advanced Notes: The idx parameters should be freed by the calling routine with PetscFree() .seealso MeshDestroy(), MeshCreateGlobalVector(), MeshCreate() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshGetGlobalIndices(Mesh mesh,PetscInt *idx[]) { SETERRQ(PETSC_ERR_SUP, ""); } #undef __FUNCT__ #define __FUNCT__ "MeshCreateGlobalScatter" PetscErrorCode PETSCDM_DLLEXPORT MeshCreateGlobalScatter(Mesh mesh, SectionReal section, VecScatter *scatter) { ALE::Obj m; ALE::Obj s; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionRealGetSection(section, s);CHKERRQ(ierr); ierr = MeshCreateGlobalScatter(m, s, scatter);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreateGlobalScatter" template PetscErrorCode PETSCDM_DLLEXPORT MeshCreateGlobalScatter(const ALE::Obj& m, const ALE::Obj
& s, VecScatter *scatter) { typedef ALE::Mesh::real_section_type::index_type index_type; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr); const ALE::Mesh::real_section_type::chart_type& chart = s->getChart(); const ALE::Obj& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s); int *localIndices, *globalIndices; int localSize = s->size(); int localIndx = 0, globalIndx = 0; Vec globalVec, localVec; IS localIS, globalIS; ierr = VecCreate(m->comm(), &globalVec);CHKERRQ(ierr); ierr = VecSetSizes(globalVec, globalOrder->getLocalSize(), PETSC_DETERMINE);CHKERRQ(ierr); ierr = VecSetFromOptions(globalVec);CHKERRQ(ierr); // Loop over all local points ierr = PetscMalloc(localSize*sizeof(int), &localIndices); CHKERRQ(ierr); ierr = PetscMalloc(localSize*sizeof(int), &globalIndices); CHKERRQ(ierr); for(ALE::Mesh::real_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) { // Map local indices to global indices s->getIndices(*p_iter, localIndices, &localIndx, 0, true); s->getIndices(*p_iter, globalOrder, globalIndices, &globalIndx, 0, true); } if (localIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should be %d", localIndx, localSize); if (globalIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should be %d", globalIndx, localSize); if (m->debug()) { globalOrder->view("Global Order"); for(int i = 0; i < localSize; ++i) { printf("[%d] localIndex[%d]: %d globalIndex[%d]: %d\n", m->commRank(), i, localIndices[i], i, globalIndices[i]); } } ierr = ISCreateGeneral(PETSC_COMM_SELF, localSize, localIndices, &localIS);CHKERRQ(ierr); ierr = ISCreateGeneral(PETSC_COMM_SELF, localSize, globalIndices, &globalIS);CHKERRQ(ierr); ierr = PetscFree(localIndices);CHKERRQ(ierr); ierr = PetscFree(globalIndices);CHKERRQ(ierr); ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, s->sizeWithBC(), s->restrict(), &localVec);CHKERRQ(ierr); ierr = VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);CHKERRQ(ierr); ierr = ISDestroy(globalIS);CHKERRQ(ierr); ierr = ISDestroy(localIS);CHKERRQ(ierr); ierr = VecDestroy(localVec);CHKERRQ(ierr); ierr = VecDestroy(globalVec);CHKERRQ(ierr); ierr = PetscLogEventEnd(Mesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetGlobalScatter" PetscErrorCode PETSCDM_DLLEXPORT MeshGetGlobalScatter(Mesh mesh, VecScatter *scatter) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); PetscValidPointer(scatter, 2); if (!mesh->globalScatter) { SectionReal section; ierr = MeshGetSectionReal(mesh, "default", §ion);CHKERRQ(ierr); ierr = MeshCreateGlobalScatter(mesh, section, &mesh->globalScatter);CHKERRQ(ierr); ierr = SectionRealDestroy(section);CHKERRQ(ierr); } *scatter = mesh->globalScatter; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetLocalFunction" PetscErrorCode PETSCDM_DLLEXPORT MeshSetLocalFunction(Mesh mesh, PetscErrorCode (*lf)(Mesh, SectionReal, SectionReal, void *)) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); mesh->lf = lf; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetLocalJacobian" PetscErrorCode PETSCDM_DLLEXPORT MeshSetLocalJacobian(Mesh mesh, PetscErrorCode (*lj)(Mesh, SectionReal, Mat, void *)) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); mesh->lj = lj; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshFormFunction" PetscErrorCode PETSCDM_DLLEXPORT MeshFormFunction(Mesh mesh, SectionReal X, SectionReal F, void *ctx) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); PetscValidHeaderSpecific(X, SECTIONREAL_COOKIE, 2); PetscValidHeaderSpecific(F, SECTIONREAL_COOKIE, 3); if (mesh->lf) { ierr = (*mesh->lf)(mesh, X, F, ctx);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshFormJacobian" PetscErrorCode PETSCDM_DLLEXPORT MeshFormJacobian(Mesh mesh, SectionReal X, Mat J, void *ctx) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); PetscValidHeaderSpecific(X, SECTIONREAL_COOKIE, 2); PetscValidHeaderSpecific(J, MAT_COOKIE, 3); if (mesh->lj) { ierr = (*mesh->lj)(mesh, X, J, ctx);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshInterpolatePoints" // Here we assume: // - Assumes 3D and tetrahedron // - The section takes values on vertices and is P1 // - Points have the same dimension as the mesh // - All values have the same dimension PetscErrorCode PETSCDM_DLLEXPORT MeshInterpolatePoints(Mesh mesh, SectionReal section, int numPoints, double *points, double **values) { Obj m; Obj s; double *v0, *J, *invJ, detJ; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionRealGetSection(section, s);CHKERRQ(ierr); const Obj& coordinates = m->getRealSection("coordinates"); int embedDim = coordinates->getFiberDimension(*m->depthStratum(0)->begin()); int dim = s->getFiberDimension(*m->depthStratum(0)->begin()); ierr = PetscMalloc3(embedDim,double,&v0,embedDim*embedDim,double,&J,embedDim*embedDim,double,&invJ);CHKERRQ(ierr); ierr = PetscMalloc(numPoints*dim * sizeof(double), &values);CHKERRQ(ierr); for(int p = 0; p < numPoints; p++) { double *point = &points[p*embedDim]; ALE::Mesh::point_type e = m->locatePoint(point); const ALE::Mesh::real_section_type::value_type *coeff = s->restrictPoint(e); m->computeElementGeometry(coordinates, e, v0, J, invJ, detJ); double xi = (invJ[0*embedDim+0]*(point[0] - v0[0]) + invJ[0*embedDim+1]*(point[1] - v0[1]) + invJ[0*embedDim+2]*(point[2] - v0[2]))*0.5; double eta = (invJ[1*embedDim+0]*(point[0] - v0[0]) + invJ[1*embedDim+1]*(point[1] - v0[1]) + invJ[1*embedDim+2]*(point[2] - v0[2]))*0.5; double zeta = (invJ[2*embedDim+0]*(point[0] - v0[0]) + invJ[2*embedDim+1]*(point[1] - v0[1]) + invJ[2*embedDim+2]*(point[2] - v0[2]))*0.5; for(int d = 0; d < dim; d++) { (*values)[p*dim+d] = coeff[0*dim+d]*(1 - xi - eta - zeta) + coeff[1*dim+d]*xi + coeff[2*dim+d]*eta + coeff[3*dim+d]*zeta; } } ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetMaximumDegree" /*@ MeshGetMaximumDegree - Return the maximum degree of any mesh vertex Collective on mesh Input Parameter: . mesh - The Mesh Output Parameter: . maxDegree - The maximum number of edges at any vertex Level: beginner .seealso: MeshCreate() @*/ PetscErrorCode MeshGetMaximumDegree(Mesh mesh, PetscInt *maxDegree) { Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); const ALE::Obj& vertices = m->depthStratum(0); const ALE::Obj& sieve = m->getSieve(); PetscInt maxDeg = -1; for(ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) { maxDeg = PetscMax(maxDeg, (PetscInt) sieve->support(*v_iter)->size()); } *maxDegree = maxDeg; PetscFunctionReturn(0); } EXTERN PetscErrorCode assembleFullField(VecScatter, Vec, Vec, InsertMode); #undef __FUNCT__ #define __FUNCT__ "restrictVector" /*@C restrictVector - Insert values from a global vector into a local ghosted vector Collective on g Input Parameters: + g - The global vector . l - The local vector - mode - either ADD_VALUES or INSERT_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values Level: beginner .seealso: MatSetOption() @*/ PetscErrorCode restrictVector(Vec g, Vec l, InsertMode mode) { VecScatter injection; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_restrictVector,0,0,0,0);CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);CHKERRQ(ierr); if (injection) { ierr = VecScatterBegin(injection, g, l, mode, SCATTER_REVERSE); ierr = VecScatterEnd(injection, g, l, mode, SCATTER_REVERSE); } else { if (mode == INSERT_VALUES) { ierr = VecCopy(g, l);CHKERRQ(ierr); } else { ierr = VecAXPY(l, 1.0, g);CHKERRQ(ierr); } } ierr = PetscLogEventEnd(Mesh_restrictVector,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "assembleVectorComplete" /*@C assembleVectorComplete - Insert values from a local ghosted vector into a global vector Collective on g Input Parameters: + g - The global vector . l - The local vector - mode - either ADD_VALUES or INSERT_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values Level: beginner .seealso: MatSetOption() @*/ PetscErrorCode assembleVectorComplete(Vec g, Vec l, InsertMode mode) { VecScatter injection; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_assembleVectorComplete,0,0,0,0);CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) g, "injection", (PetscObject *) &injection);CHKERRQ(ierr); if (injection) { ierr = VecScatterBegin(injection, l, g, mode, SCATTER_FORWARD);CHKERRQ(ierr); ierr = VecScatterEnd(injection, l, g, mode, SCATTER_FORWARD);CHKERRQ(ierr); } else { if (mode == INSERT_VALUES) { ierr = VecCopy(l, g);CHKERRQ(ierr); } else { ierr = VecAXPY(g, 1.0, l);CHKERRQ(ierr); } } ierr = PetscLogEventEnd(Mesh_assembleVectorComplete,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "assembleVector" /*@C assembleVector - Insert values into a vector Collective on A Input Parameters: + b - the vector . e - The element number . v - The values - mode - either ADD_VALUES or INSERT_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values Level: beginner .seealso: VecSetOption() @*/ PetscErrorCode assembleVector(Vec b, PetscInt e, PetscScalar v[], InsertMode mode) { Mesh mesh; ALE::Obj m; PetscInt firstElement; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_assembleVector,0,0,0,0);CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) b, "mesh", (PetscObject *) &mesh);CHKERRQ(ierr); ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); //firstElement = elementBundle->getLocalSizes()[bundle->getCommRank()]; firstElement = 0; // Must relate b to field if (mode == INSERT_VALUES) { m->update(m->getRealSection(std::string("x")), ALE::Mesh::point_type(e + firstElement), v); } else { m->updateAdd(m->getRealSection(std::string("x")), ALE::Mesh::point_type(e + firstElement), v); } ierr = PetscLogEventEnd(Mesh_assembleVector,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "updateOperator" PetscErrorCode updateOperator(Mat A, const ALE::Obj& m, const ALE::Obj& section, const ALE::Obj& globalOrder, const ALE::Mesh::point_type& e, PetscScalar array[], InsertMode mode) { PetscErrorCode ierr; PetscFunctionBegin; const ALE::Mesh::indices_type indicesBlock = m->getIndices(section, e, globalOrder); const PetscInt *indices = indicesBlock.first; const int& numIndices = indicesBlock.second; ierr = PetscLogEventBegin(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); if (section->debug()) { printf("[%d]mat for element %d\n", section->commRank(), e); for(int i = 0; i < numIndices; i++) { printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]); } for(int i = 0; i < numIndices; i++) { printf("[%d]", section->commRank()); for(int j = 0; j < numIndices; j++) { printf(" %g", array[i*numIndices+j]); } printf("\n"); } } ierr = MatSetValues(A, numIndices, indices, numIndices, indices, array, mode); if (ierr) { printf("[%d]ERROR in updateOperator: point %d\n", section->commRank(), e); for(int i = 0; i < numIndices; i++) { printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]); } CHKERRQ(ierr); } ierr = PetscLogEventEnd(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "updateOperatorGeneral" PetscErrorCode updateOperatorGeneral(Mat A, const ALE::Obj& rowM, const ALE::Obj& rowSection, const ALE::Obj& rowGlobalOrder, const ALE::Mesh::point_type& rowE, const ALE::Obj& colM, const ALE::Obj& colSection, const ALE::Obj& colGlobalOrder, const ALE::Mesh::point_type& colE, PetscScalar array[], InsertMode mode) { PetscErrorCode ierr; PetscFunctionBegin; const ALE::Mesh::indices_type rowIndicesBlock = rowM->getIndices(rowSection, rowE, rowGlobalOrder); const PetscInt *rowIndices = rowIndicesBlock.first; const int& numRowIndices = rowIndicesBlock.second; const ALE::Mesh::indices_type colIndicesBlock = colM->getIndices(colSection, colE, colGlobalOrder); const PetscInt *colIndices = colIndicesBlock.first; const int& numColIndices = colIndicesBlock.second; ierr = PetscLogEventBegin(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); if (rowSection->debug()) { printf("[%d]mat for elements %d %d\n", rowSection->commRank(), rowE, colE); for(int i = 0; i < numRowIndices; i++) { printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]); } } if (colSection->debug()) { for(int i = 0; i < numColIndices; i++) { printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]); } for(int i = 0; i < numRowIndices; i++) { printf("[%d]", rowSection->commRank()); for(int j = 0; j < numColIndices; j++) { printf(" %g", array[i*numColIndices+j]); } printf("\n"); } } ierr = MatSetValues(A, numRowIndices, rowIndices, numColIndices, colIndices, array, mode); if (ierr) { printf("[%d]ERROR in updateOperator: points %d %d\n", colSection->commRank(), rowE, colE); for(int i = 0; i < numRowIndices; i++) { printf("[%d]mat row indices[%d] = %d\n", rowSection->commRank(), i, rowIndices[i]); } for(int i = 0; i < numColIndices; i++) { printf("[%d]mat col indices[%d] = %d\n", colSection->commRank(), i, colIndices[i]); } CHKERRQ(ierr); } ierr = PetscLogEventEnd(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "assembleMatrix" /*@C assembleMatrix - Insert values into a matrix Collective on A Input Parameters: + A - the matrix . e - The element number . v - The values - mode - either ADD_VALUES or INSERT_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values Level: beginner .seealso: MatSetOption() @*/ PetscErrorCode assembleMatrix(Mat A, PetscInt e, PetscScalar v[], InsertMode mode) { Mesh mesh; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_assembleMatrix,0,0,0,0);CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) A, "mesh", (PetscObject *) &mesh);CHKERRQ(ierr); try { if (!mesh->mcompat.isNull()) { Obj m; ierr = MeshCompatGetMesh(mesh, m);CHKERRQ(ierr); const ALE::Obj& topology = m->getTopology(); const ALE::Obj& cNumbering = m->getFactory()->getLocalNumbering(topology, 0, topology->depth()); const ALE::Obj& s = m->getRealSection("default"); const ALE::Obj& globalOrder = m->getFactory()->getGlobalOrder(topology, 0, "default", s->getAtlas()); if (m->debug()) { std::cout << "Assembling matrix for element number " << e << " --> point " << cNumbering->getPoint(e) << std::endl; } ierr = updateOperatorCompat(A, s, globalOrder, cNumbering->getPoint(e), v, mode);CHKERRQ(ierr); } else { Obj m; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); const ALE::Obj& cNumbering = m->getFactory()->getLocalNumbering(m, m->depth()); const ALE::Obj& s = m->getRealSection("default"); const ALE::Obj& globalOrder = m->getFactory()->getGlobalOrder(m, "default", s); if (m->debug()) { std::cout << "Assembling matrix for element number " << e << " --> point " << cNumbering->getPoint(e) << std::endl; } ierr = updateOperator(A, m, s, globalOrder, cNumbering->getPoint(e), v, mode);CHKERRQ(ierr); } } catch (ALE::Exception e) { std::cout << e.msg() << std::endl; } ierr = PetscLogEventEnd(Mesh_assembleMatrix,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "preallocateOperator" template PetscErrorCode preallocateOperator(const ALE::Obj& mesh, const int bs, const ALE::Obj& atlas, const ALE::Obj& globalOrder, Mat A) { typedef ALE::SieveAlg sieve_alg_type; MPI_Comm comm = mesh->comm(); const ALE::Obj adjBundle = new ALE::Mesh(comm, mesh->debug()); const ALE::Obj adjGraph = new ALE::Mesh::sieve_type(comm, mesh->debug()); PetscInt numLocalRows, firstRow; PetscInt *dnz, *onz; PetscErrorCode ierr; PetscFunctionBegin; adjBundle->setSieve(adjGraph); numLocalRows = globalOrder->getLocalSize(); firstRow = globalOrder->getGlobalOffsets()[mesh->commRank()]; ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr); /* Create local adjacency graph */ /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */ const typename Atlas::chart_type& chart = atlas->getChart(); for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const Obj& star = sieve_alg_type::star(mesh, *c_iter); for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) { const Obj& closure = sieve_alg_type::closure(mesh, *s_iter); for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) { adjGraph->addCone(*cl_iter, *c_iter); } } } /* Distribute adjacency graph */ adjBundle->constructOverlap(); typedef typename ALE::Mesh::sieve_type::point_type point_type; typedef typename ALE::Mesh::send_overlap_type send_overlap_type; typedef typename ALE::Mesh::recv_overlap_type recv_overlap_type; typedef typename ALE::Field > send_section_type; typedef typename ALE::Field > recv_section_type; const Obj& vertexSendOverlap = mesh->getSendOverlap(); const Obj& vertexRecvOverlap = mesh->getRecvOverlap(); const Obj nbrSendOverlap = new send_overlap_type(comm, mesh->debug()); const Obj nbrRecvOverlap = new recv_overlap_type(comm, mesh->debug()); const Obj sendSection = new send_section_type(comm, mesh->debug()); const Obj recvSection = new recv_section_type(comm, sendSection->getTag(), mesh->debug()); ALE::Distribution::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjBundle, sendSection, recvSection); /* Distribute indices for new points */ ALE::Distribution::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap); mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, true); /* Read out adjacency graph */ const ALE::Obj graph = adjBundle->getSieve(); ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr); ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr); for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const typename Atlas::point_type& point = *c_iter; if (globalOrder->isLocal(point)) { const ALE::Obj& adj = graph->cone(point); const ALE::Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0]; const int row = rIdx.prefix; const int rSize = rIdx.index/bs; if (rSize == 0) continue; for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) { const ALE::Mesh::point_type& neighbor = *v_iter; const ALE::Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0]; const int& cSize = cIdx.index/bs; if (cSize > 0) { if (globalOrder->isLocal(neighbor)) { for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;} } else { for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;} } } } } } if (mesh->debug()) { int rank = mesh->commRank(); for(int r = 0; r < numLocalRows/bs; r++) { std::cout << "["<& mesh, const int bs, const ALE::Obj& atlas, const ALE::Obj& globalOrder, Mat A) { return preallocateOperator(mesh, bs, atlas, globalOrder, A); } /******************************** C Wrappers **********************************/ #undef __FUNCT__ #define __FUNCT__ "WriteVTKHeader" PetscErrorCode WriteVTKHeader(PetscViewer viewer) { return VTKViewer::writeHeader(viewer); } #undef __FUNCT__ #define __FUNCT__ "WriteVTKVertices" PetscErrorCode WriteVTKVertices(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return VTKViewer::writeVertices(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WriteVTKElements" PetscErrorCode WriteVTKElements(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return VTKViewer::writeElements(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePCICEVertices" PetscErrorCode WritePCICEVertices(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PCICE::Viewer::writeVertices(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePCICEElements" PetscErrorCode WritePCICEElements(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PCICE::Viewer::writeElements(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePCICERestart" PetscErrorCode WritePCICERestart(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PCICE::Viewer::writeRestart(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "MeshCreatePCICE" /*@C MeshCreatePCICE - Create a Mesh from PCICE files. Not Collective Input Parameters: + dim - The topological mesh dimension . coordFilename - The file containing vertex coordinates . adjFilename - The file containing the vertices for each element . interpolate - The flag for construction of intermediate elements . bcFilename - The file containing the boundary topology and conditions . numBdFaces - The number of boundary faces (or edges) - numBdVertices - The number of boundary vertices Output Parameter: . mesh - The Mesh object Level: beginner .keywords: mesh, PCICE .seealso: MeshCreate() @*/ PetscErrorCode MeshCreatePCICE(MPI_Comm comm, const int dim, const char coordFilename[], const char adjFilename[], PetscTruth interpolate, const char bcFilename[], Mesh *mesh) { ALE::Obj m; PetscInt debug = 0; PetscTruth flag; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshCreate(comm, mesh);CHKERRQ(ierr); ierr = PetscOptionsGetInt(PETSC_NULL, "-debug", &debug, &flag);CHKERRQ(ierr); try { m = ALE::PCICE::Builder::readMesh(comm, dim, std::string(coordFilename), std::string(adjFilename), false, interpolate, debug); if (debug) {m->view("Mesh");} } catch(ALE::Exception e) { SETERRQ(PETSC_ERR_FILE_OPEN, e.message()); } if (bcFilename) { ALE::PCICE::Builder::readBoundary(m, std::string(bcFilename)); } ierr = MeshSetMesh(*mesh, m);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCreatePyLith" /*@C MeshCreatePyLith - Create a Mesh from PyLith files. Not Collective Input Parameters: + dim - The topological mesh dimension . baseFilename - The basename for mesh files . zeroBase - Use 0 to start numbering - interpolate - The flag for mesh interpolation Output Parameter: . mesh - The Mesh object Level: beginner .keywords: mesh, PCICE .seealso: MeshCreate() @*/ PetscErrorCode MeshCreatePyLith(MPI_Comm comm, const int dim, const char baseFilename[], PetscTruth zeroBase, PetscTruth interpolate, Mesh *mesh) { ALE::Obj m; PetscInt debug = 0; PetscTruth flag; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshCreate(comm, mesh);CHKERRQ(ierr); ierr = PetscOptionsGetInt(PETSC_NULL, "-debug", &debug, &flag);CHKERRQ(ierr); try { m = ALE::PyLith::Builder::readMesh(comm, dim, std::string(baseFilename), zeroBase, interpolate, debug); } catch(ALE::Exception e) { SETERRQ(PETSC_ERR_FILE_OPEN, e.message()); } ierr = MeshSetMesh(*mesh, m);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetCoordinates" /*@C MeshGetCoordinates - Creates an array holding the coordinates. Not Collective Input Parameter: + mesh - The Mesh object - columnMajor - Flag for column major order Output Parameter: + numVertices - The number of vertices . dim - The embedding dimension - coords - The array holding local coordinates Level: intermediate .keywords: mesh, coordinates .seealso: MeshCreate() @*/ PetscErrorCode MeshGetCoordinates(Mesh mesh, PetscTruth columnMajor, PetscInt *numVertices, PetscInt *dim, PetscReal *coords[]) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ALE::PCICE::Builder::outputVerticesLocal(m, numVertices, dim, coords, columnMajor); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetElements" /*@C MeshGetElements - Creates an array holding the vertices on each element. Not Collective Input Parameters: + mesh - The Mesh object - columnMajor - Flag for column major order Output Parameters: + numElements - The number of elements . numCorners - The number of vertices per element - vertices - The array holding vertices on each local element Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshGetElements(Mesh mesh, PetscTruth columnMajor, PetscInt *numElements, PetscInt *numCorners, PetscInt *vertices[]) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ALE::PCICE::Builder::outputElementsLocal(m, numElements, numCorners, vertices, columnMajor); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshDistribute" /*@C MeshDistribute - Distributes the mesh and any associated sections. Not Collective Input Parameter: + serialMesh - The original Mesh object - partitioner - The partitioning package, or NULL for the default Output Parameter: . parallelMesh - The distributed Mesh object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate(), MeshDistributeByFace() @*/ PetscErrorCode MeshDistribute(Mesh serialMesh, const char partitioner[], Mesh *parallelMesh) { ALE::Obj oldMesh; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(serialMesh, oldMesh);CHKERRQ(ierr); ierr = MeshCreate(oldMesh->comm(), parallelMesh);CHKERRQ(ierr); if (partitioner == NULL) { ALE::Obj newMesh = ALE::Distribution::distributeMesh(oldMesh); ierr = MeshSetMesh(*parallelMesh, newMesh);CHKERRQ(ierr); } else { ALE::Obj newMesh = ALE::Distribution::distributeMesh(oldMesh, 0, partitioner); ierr = MeshSetMesh(*parallelMesh, newMesh);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshDistributeByFace" /*@C MeshDistribute - Distributes the mesh and any associated sections. Not Collective Input Parameter: + serialMesh - The original Mesh object - partitioner - The partitioning package, or NULL for the default Output Parameter: . parallelMesh - The distributed Mesh object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate(), MeshDistribute() @*/ PetscErrorCode MeshDistributeByFace(Mesh serialMesh, const char partitioner[], Mesh *parallelMesh) { ALE::Obj oldMesh; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(serialMesh, oldMesh);CHKERRQ(ierr); ierr = MeshCreate(oldMesh->comm(), parallelMesh);CHKERRQ(ierr); if (partitioner == NULL) { ALE::Obj newMesh = ALE::Distribution::distributeMesh(oldMesh, 1); ierr = MeshSetMesh(*parallelMesh, newMesh);CHKERRQ(ierr); } else { ALE::Obj newMesh = ALE::Distribution::distributeMesh(oldMesh, 1, partitioner); ierr = MeshSetMesh(*parallelMesh, newMesh);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGenerate" /*@C MeshGenerate - Generates a mesh. Not Collective Input Parameters: + boundary - The Mesh boundary object - interpolate - Flag to create intermediate mesh elements Output Parameter: . mesh - The Mesh object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate(), MeshRefine() @*/ PetscErrorCode MeshGenerate(Mesh boundary, PetscTruth interpolate, Mesh *mesh) { ALE::Obj mB; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(boundary, mB);CHKERRQ(ierr); ierr = MeshCreate(mB->comm(), mesh);CHKERRQ(ierr); ALE::Obj m = ALE::Generator::generateMesh(mB, interpolate); ierr = MeshSetMesh(*mesh, m);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRefine" /*@C MeshRefine - Refines the mesh. Not Collective Input Parameters: + mesh - The original Mesh object . refinementLimit - The maximum size of any cell - interpolate - Flag to create intermediate mesh elements Output Parameter: . refinedMesh - The refined Mesh object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate(), MeshGenerate() @*/ PetscErrorCode MeshRefine(Mesh mesh, double refinementLimit, PetscTruth interpolate, Mesh *refinedMesh) { ALE::Obj oldMesh; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, oldMesh);CHKERRQ(ierr); ierr = MeshCreate(oldMesh->comm(), refinedMesh);CHKERRQ(ierr); ALE::Obj newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, interpolate); ierr = MeshSetMesh(*refinedMesh, newMesh);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRefine_Mesh" PetscErrorCode MeshRefine_Mesh(Mesh mesh, MPI_Comm comm, Mesh *refinedMesh) { ALE::Obj oldMesh; double refinementLimit; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, oldMesh);CHKERRQ(ierr); ierr = MeshCreate(comm, refinedMesh);CHKERRQ(ierr); refinementLimit = oldMesh->getMaxVolume()/2.0; ALE::Obj newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, true); ierr = MeshSetMesh(*refinedMesh, newMesh);CHKERRQ(ierr); const ALE::Obj& s = newMesh->getRealSection("default"); newMesh->setDiscretization(oldMesh->getDiscretization()); newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition()); newMesh->setupField(s); PetscFunctionReturn(0); } #include "Hierarchy.hh" #undef __FUNCT__ #define __FUNCT__ "MeshCoarsenHierarchy" /*@C MeshCoarsenHierarchy - Coarsens the mesh into a hierarchy. Not Collective Input Parameters: + mesh - The original Mesh object . numLevels - The number of . coarseningFactor - The expansion factor for coarse meshes - interpolate - Flag to create intermediate mesh elements Output Parameter: . coarseHierarchy - The coarse Mesh objects Level: intermediate .keywords: mesh, elements .seealso: MeshCreate(), MeshGenerate() @*/ PetscErrorCode MeshCoarsenHierarchy(Mesh mesh, int numLevels, double coarseningFactor, PetscTruth interpolate, Mesh **coarseHierarchy) { ALE::Obj oldMesh; PetscErrorCode ierr; PetscFunctionBegin; if (numLevels < 1) { *coarseHierarchy = PETSC_NULL; PetscFunctionReturn(0); } ierr = MeshGetMesh(mesh, oldMesh);CHKERRQ(ierr); ierr = PetscMalloc((numLevels+1) * sizeof(Mesh), coarseHierarchy);CHKERRQ(ierr); for (int i = 0; i < numLevels+1; i++) { ierr = MeshCreate(oldMesh->comm(), &(*coarseHierarchy)[i]);CHKERRQ(ierr); } ierr = MeshSpacingFunction(mesh);CHKERRQ(ierr); ierr = MeshCreateHierarchyLabel(mesh, coarseningFactor, numLevels+1, *coarseHierarchy); #if 0 if (oldMesh->getDimension() != 2) SETERRQ(PETSC_ERR_SUP, "Coarsening only works in two dimensions right now"); ierr = ALE::Coarsener::IdentifyBoundary(oldMesh, 2);CHKERRQ(ierr); ierr = ALE::Coarsener::make_coarsest_boundary(oldMesh, 2, numLevels+1);CHKERRQ(ierr); ierr = ALE::Coarsener::CreateSpacingFunction(oldMesh, 2);CHKERRQ(ierr); ierr = ALE::Coarsener::CreateCoarsenedHierarchyNew(oldMesh, 2, numLevels, coarseningFactor);CHKERRQ(ierr); ierr = PetscMalloc(numLevels * sizeof(Mesh),coarseHierarchy);CHKERRQ(ierr); for(int l = 0; l < numLevels; l++) { ALE::Obj newMesh = new ALE::Mesh(oldMesh->comm(), oldMesh->debug()); const ALE::Obj& s = newMesh->getRealSection("default"); ierr = MeshCreate(oldMesh->comm(), &(*coarseHierarchy)[l]);CHKERRQ(ierr); newMesh->getTopology()->setPatch(0, oldMesh->getTopology()->getPatch(l+1)); newMesh->setDiscretization(oldMesh->getDiscretization()); newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition()); newMesh->setupField(s); ierr = MeshSetMesh((*coarseHierarchy)[l], newMesh);CHKERRQ(ierr); } #endif PetscFunctionReturn(0); } PetscErrorCode MeshCoarsenHierarchy_Mesh(Mesh mesh, int numLevels, Mesh **coarseHierarchy) { PetscErrorCode ierr; double cfactor = 1.5; PetscFunctionBegin; ierr = PetscOptionsReal("-dmmg_coarsen_factor", "The coarsening factor", PETSC_NULL, cfactor, &cfactor, PETSC_NULL);CHKERRQ(ierr); ierr = MeshCoarsenHierarchy(mesh, numLevels, cfactor, PETSC_FALSE, coarseHierarchy);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetInterpolation_Mesh" /* This method only handle P_1 discretizations at present. */ PetscErrorCode MeshGetInterpolation_Mesh(Mesh dmCoarse, Mesh dmFine, Mat *interpolation, Vec *scaling) { ALE::Obj coarse; ALE::Obj fine; Mesh coarseMesh = (Mesh) dmCoarse; Mesh fineMesh = (Mesh) dmFine; Mat P; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(fineMesh, fine);CHKERRQ(ierr); ierr = MeshGetMesh(coarseMesh, coarse);CHKERRQ(ierr); const ALE::Obj& coarseCoordinates = coarse->getRealSection("coordinates"); const ALE::Obj& fineCoordinates = fine->getRealSection("coordinates"); const ALE::Obj& vertices = fine->depthStratum(0); const ALE::Obj& sCoarse = coarse->getRealSection("default"); const ALE::Obj& sFine = fine->getRealSection("default"); const ALE::Obj& coarseOrder = coarse->getFactory()->getGlobalOrder(coarse, "default", sCoarse); const ALE::Obj& fineOrder = fine->getFactory()->getGlobalOrder(fine, "default", sFine); const int dim = coarse->getDimension(); double *v0, *J, *invJ, detJ, *refCoords, *values; ierr = MatCreate(fine->comm(), &P);CHKERRQ(ierr); ierr = MatSetSizes(P, sFine->size(), sCoarse->size(), PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); ierr = MatSetFromOptions(P);CHKERRQ(ierr); ierr = MatSetUpPreallocation(P);CHKERRQ(ierr); ierr = PetscMalloc5(dim,double,&v0,dim*dim,double,&J,dim*dim,double,&invJ,dim,double,&refCoords,dim+1,double,&values);CHKERRQ(ierr); for(ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) { const ALE::Mesh::real_section_type::value_type *coords = fineCoordinates->restrictPoint(*v_iter); ALE::Mesh::point_type coarseCell; ALE::Mesh::point_type cellguess = -1; if (fine->hasLabel("prolongation")) { cellguess = fine->getValue(fine->getLabel("prolongation"), *v_iter); coarseCell = coarse->locatePoint(coords, cellguess); } else { coarseCell = coarse->locatePoint(coords); } // coarseCell = coarse->locatePoint(coords); if (coarseCell == -1) { // do NO CORRECTION! } else { coarse->computeElementGeometry(coarseCoordinates, coarseCell, v0, J, invJ, detJ); for(int d = 0; d < dim; ++d) { refCoords[d] = 0.0; for(int e = 0; e < dim; ++e) { refCoords[d] += invJ[d*dim+e]*(coords[e] - v0[e]); } refCoords[d] -= 1.0; } values[0] = -(refCoords[0] + refCoords[1])/2.0; values[1] = 0.5*(refCoords[0] + 1.0); values[2] = 0.5*(refCoords[1] + 1.0); ierr = updateOperatorGeneral(P, fine, sFine, fineOrder, *v_iter, coarse, sCoarse, coarseOrder, coarseCell, values, INSERT_VALUES);CHKERRQ(ierr); } } ierr = PetscFree5(v0,J,invJ,refCoords,values);CHKERRQ(ierr); ierr = MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); *interpolation = P; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshHasSectionReal" /*@C MeshHasSectionReal - Determines whether this mesh has a SectionReal with the given name. Not Collective Input Parameters: + mesh - The Mesh object - name - The section name Output Parameter: . flag - True if the SectionReal is present in the Mesh Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshHasSectionReal(Mesh mesh, const char name[], PetscTruth *flag) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); *flag = (PetscTruth) m->hasRealSection(std::string(name)); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetSectionReal" /*@C MeshGetSectionReal - Returns a SectionReal of the given name from the Mesh. Collective on Mesh Input Parameters: + mesh - The Mesh object - name - The section name Output Parameter: . section - The SectionReal Note: The section is a new object, and must be destroyed by the user Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshGetSectionReal(Mesh mesh, const char name[], SectionReal *section) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionRealCreate(m->comm(), section);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) *section, name);CHKERRQ(ierr); ierr = SectionRealSetSection(*section, m->getRealSection(std::string(name)));CHKERRQ(ierr); ierr = SectionRealSetBundle(*section, m);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetSectionReal" /*@C MeshSetSectionReal - Puts a SectionReal of the given name into the Mesh. Collective on Mesh Input Parameters: + mesh - The Mesh object - section - The SectionReal Note: This takes the section name from the PETSc object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshSetSectionReal(Mesh mesh, SectionReal section) { ALE::Obj m; ALE::Obj s; const char *name; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = PetscObjectGetName((PetscObject) section, &name);CHKERRQ(ierr); ierr = SectionRealGetSection(section, s);CHKERRQ(ierr); m->setRealSection(std::string(name), s); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshHasSectionInt" /*@C MeshHasSectionInt - Determines whether this mesh has a SectionInt with the given name. Not Collective Input Parameters: + mesh - The Mesh object - name - The section name Output Parameter: . flag - True if the SectionInt is present in the Mesh Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshHasSectionInt(Mesh mesh, const char name[], PetscTruth *flag) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); *flag = (PetscTruth) m->hasIntSection(std::string(name)); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetSectionInt" /*@C MeshGetSectionInt - Returns a SectionInt of the given name from the Mesh. Collective on Mesh Input Parameters: + mesh - The Mesh object - name - The section name Output Parameter: . section - The SectionInt Note: The section is a new object, and must be destroyed by the user Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshGetSectionInt(Mesh mesh, const char name[], SectionInt *section) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionIntCreate(m->comm(), section);CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) *section, name);CHKERRQ(ierr); ierr = SectionIntSetSection(*section, m->getIntSection(std::string(name)));CHKERRQ(ierr); ierr = SectionIntSetBundle(*section, m);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetSectionInt" /*@C MeshSetSectionInt - Puts a SectionInt of the given name into the Mesh. Collective on Mesh Input Parameters: + mesh - The Mesh object - section - The SectionInt Note: This takes the section name from the PETSc object Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode MeshSetSectionInt(Mesh mesh, SectionInt section) { ALE::Obj m; ALE::Obj s; const char *name; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = PetscObjectGetName((PetscObject) section, &name);CHKERRQ(ierr); ierr = SectionIntGetSection(section, s);CHKERRQ(ierr); m->setIntSection(std::string(name), s); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "SectionGetArray" /*@C SectionGetArray - Returns the array underlying the Section. Not Collective Input Parameters: + mesh - The Mesh object - name - The section name Output Parameters: + numElements - The number of mesh element with values . fiberDim - The number of values per element - array - The array Level: intermediate .keywords: mesh, elements .seealso: MeshCreate() @*/ PetscErrorCode SectionGetArray(Mesh mesh, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[]) { ALE::Obj m; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); const Obj& section = m->getRealSection(std::string(name)); if (section->size() == 0) { *numElements = 0; *fiberDim = 0; *array = NULL; PetscFunctionReturn(0); } const ALE::Mesh::real_section_type::chart_type& chart = section->getChart(); /* const int depth = m->depth(*chart.begin()); */ /* *numElements = m->depthStratum(depth)->size(); */ /* *fiberDim = section->getFiberDimension(*chart.begin()); */ /* *array = (PetscScalar *) m->restrict(section); */ int fiberDimMin = section->getFiberDimension(*chart.begin()); int numElem = 0; for(ALE::Mesh::real_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const int fiberDim = section->getFiberDimension(*c_iter); if (fiberDim < fiberDimMin) fiberDimMin = fiberDim; } for(ALE::Mesh::real_section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const int fiberDim = section->getFiberDimension(*c_iter); numElem += fiberDim/fiberDimMin; } *numElements = numElem; *fiberDim = fiberDimMin; *array = (PetscScalar *) section->restrict(); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "WritePyLithVertices" PetscErrorCode WritePyLithVertices(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PyLith::Viewer::writeVertices(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePyLithElements" PetscErrorCode WritePyLithElements(Mesh mesh, SectionInt material, PetscViewer viewer) { ALE::Obj m; ALE::Obj s; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionIntGetSection(material, s);CHKERRQ(ierr); return ALE::PyLith::Viewer::writeElements(m, s, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePyLithVerticesLocal" PetscErrorCode WritePyLithVerticesLocal(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PyLith::Viewer::writeVerticesLocal(m, viewer); } #undef __FUNCT__ #define __FUNCT__ "WritePyLithElementsLocal" PetscErrorCode WritePyLithElementsLocal(Mesh mesh, SectionInt material, PetscViewer viewer) { ALE::Obj m; ALE::Obj s; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); ierr = SectionIntGetSection(material, s);CHKERRQ(ierr); return ALE::PyLith::Viewer::writeElementsLocal(m, s, viewer); } #if 0 #undef __FUNCT__ #define __FUNCT__ "WritePyLithTractionsLocal" PetscErrorCode WritePyLithTractionsLocal(Mesh mesh, PetscViewer viewer) { ALE::Obj m; PetscErrorCode ierr; ierr = MeshGetMesh(mesh, m);CHKERRQ(ierr); return ALE::PyLith::Viewer::writeTractionsLocal(m, m->getRealSection("tractions"), viewer); } #endif #undef __FUNCT__ #define __FUNCT__ "MeshCompatGetMesh" /*@C MeshCompatGetMesh - Gets the internal mesh object Not collective Input Parameter: . mesh - the mesh object Output Parameter: . m - the internal mesh object Notes: This is part of the PyLith 0.8 compatibility layer. DO NOT USE unless you are developing for that tool. Level: developer .seealso MeshCreate(), MeshSetMesh() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshCompatGetMesh(Mesh mesh, ALE::Obj& m) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); m = mesh->mcompat; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCompatSetMesh" /*@C MeshCompatSetMesh - Sets the internal mesh object Not collective Input Parameters: + mesh - the mesh object - m - the internal mesh object Notes: This is part of the PyLith 0.8 compatibility layer. DO NOT USE unless you are developing for that tool. Level: developer .seealso MeshCreate(), MeshGetMesh() @*/ PetscErrorCode PETSCDM_DLLEXPORT MeshCompatSetMesh(Mesh mesh, const ALE::Obj& m) { PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); mesh->mcompat = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "ExpandInterval" inline void ExpandInterval(const ALE::Point& interval, int indices[], int& indx) { const int end = interval.prefix + interval.index; for(int i = interval.index; i < end; i++) { indices[indx++] = i; } } #undef __FUNCT__ #define __FUNCT__ "ExpandInterval_New" inline void ExpandInterval_New(ALE::Point interval, PetscInt indices[], PetscInt *indx) { for(int i = 0; i < interval.prefix; i++) { indices[(*indx)++] = interval.index + i; } for(int i = 0; i < -interval.prefix; i++) { indices[(*indx)++] = -1; } } #undef __FUNCT__ #define __FUNCT__ "ExpandIntervals" PetscErrorCode ExpandIntervals(ALE::Obj intervals, PetscInt *indices) { int k = 0; PetscFunctionBegin; for(ALECompat::Mesh::real_section_type::IndexArray::iterator i_itor = intervals->begin(); i_itor != intervals->end(); i_itor++) { ExpandInterval_New(*i_itor, indices, &k); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCompatCreateGlobalScatter" template PetscErrorCode PETSCDM_DLLEXPORT MeshCompatCreateGlobalScatter(const ALE::Obj& m, const ALE::Obj
& s, VecScatter *scatter) { typedef ALECompat::Mesh::real_section_type::index_type index_type; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscLogEventBegin(Mesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr); const ALE::Obj& topology = m->getTopology(); const ALE::Obj& atlas = s->getAtlas(); const ALECompat::Mesh::real_section_type::patch_type patch = 0; const ALECompat::Mesh::real_section_type::atlas_type::chart_type& chart = atlas->getPatch(patch); const ALE::Obj& globalOrder = m->getFactory()->getGlobalOrder(topology, patch, s->getName(), atlas); int *localIndices, *globalIndices; int localSize = s->size(patch); int localIndx = 0, globalIndx = 0; Vec globalVec, localVec; IS localIS, globalIS; ierr = VecCreate(m->comm(), &globalVec);CHKERRQ(ierr); ierr = VecSetSizes(globalVec, globalOrder->getLocalSize(), PETSC_DETERMINE);CHKERRQ(ierr); ierr = VecSetFromOptions(globalVec);CHKERRQ(ierr); // Loop over all local points ierr = PetscMalloc(localSize*sizeof(int), &localIndices); CHKERRQ(ierr); ierr = PetscMalloc(localSize*sizeof(int), &globalIndices); CHKERRQ(ierr); for(ALECompat::Mesh::real_section_type::atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) { const ALECompat::Mesh::real_section_type::index_type& idx = atlas->restrictPoint(patch, *p_iter)[0]; // Map local indices to global indices ExpandInterval(idx, localIndices, localIndx); ExpandInterval(index_type(idx.prefix, globalOrder->getIndex(*p_iter)), globalIndices, globalIndx); } if (localIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should be %d", localIndx, localSize); if (globalIndx != localSize) SETERRQ2(PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should be %d", globalIndx, localSize); if (m->debug()) { globalOrder->view("Global Order"); for(int i = 0; i < localSize; ++i) { printf("[%d] localIndex[%d]: %d globalIndex[%d]: %d\n", m->commRank(), i, localIndices[i], i, globalIndices[i]); } } ierr = ISCreateGeneral(PETSC_COMM_SELF, localSize, localIndices, &localIS);CHKERRQ(ierr); ierr = ISCreateGeneral(PETSC_COMM_SELF, localSize, globalIndices, &globalIS);CHKERRQ(ierr); ierr = PetscFree(localIndices);CHKERRQ(ierr); ierr = PetscFree(globalIndices);CHKERRQ(ierr); ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, localSize, s->restrict(patch), &localVec);CHKERRQ(ierr); ierr = VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);CHKERRQ(ierr); ierr = ISDestroy(globalIS);CHKERRQ(ierr); ierr = ISDestroy(localIS);CHKERRQ(ierr); ierr = VecDestroy(localVec);CHKERRQ(ierr); ierr = VecDestroy(globalVec);CHKERRQ(ierr); ierr = PetscLogEventEnd(Mesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCompatGetGlobalScatter" PetscErrorCode PETSCDM_DLLEXPORT MeshCompatGetGlobalScatter(Mesh mesh, VecScatter *scatter) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE, 1); PetscValidPointer(scatter, 2); if (!mesh->globalScatter) { ALE::Obj m; ierr = MeshCompatGetMesh(mesh, m);CHKERRQ(ierr); ierr = MeshCompatCreateGlobalScatter(m, m->getRealSection("default"), &mesh->globalScatter);CHKERRQ(ierr); } *scatter = mesh->globalScatter; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "preallocateOperatorCompat" template PetscErrorCode preallocateOperatorCompat(const ALE::Obj& topology, const ALE::Obj& atlas, const ALE::Obj& globalOrder, Mat A) { typedef ALECompat::New::NumberingFactory NumberingFactory; const ALE::Obj adjGraph = new ALECompat::Mesh::sieve_type(topology->comm(), topology->debug()); const ALE::Obj adjTopology = new ALECompat::Mesh::topology_type(topology->comm(), topology->debug()); const ALECompat::Mesh::real_section_type::patch_type patch = 0; const ALE::Obj& sieve = topology->getPatch(patch); PetscInt numLocalRows, firstRow; PetscInt *dnz, *onz; PetscErrorCode ierr; PetscFunctionBegin; adjTopology->setPatch(patch, adjGraph); numLocalRows = globalOrder->getLocalSize(); firstRow = globalOrder->getGlobalOffsets()[topology->commRank()]; ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr); /* Create local adjacency graph */ /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */ const ALECompat::Mesh::real_section_type::atlas_type::chart_type& chart = atlas->getPatch(patch); for(ALECompat::Mesh::real_section_type::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const ALECompat::Mesh::real_section_type::atlas_type::point_type& point = *c_iter; adjGraph->addCone(sieve->cone(sieve->support(point)), point); } /* Distribute adjacency graph */ topology->constructOverlap(patch); const Obj& vertexSendOverlap = topology->getSendOverlap(); const Obj& vertexRecvOverlap = topology->getRecvOverlap(); const Obj nbrSendOverlap = new ALECompat::Mesh::send_overlap_type(topology->comm(), topology->debug()); const Obj nbrRecvOverlap = new ALECompat::Mesh::recv_overlap_type(topology->comm(), topology->debug()); const Obj sendSection = new ALECompat::Mesh::send_section_type(topology->comm(), topology->debug()); const Obj recvSection = new ALECompat::Mesh::recv_section_type(topology->comm(), sendSection->getTag(), topology->debug()); ALECompat::New::Distribution::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjTopology, sendSection, recvSection); /* Distribute indices for new points */ ALECompat::New::Distribution::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap); NumberingFactory::singleton(topology->debug())->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, patch, true); /* Read out adjacency graph */ const ALE::Obj graph = adjTopology->getPatch(patch); ierr = PetscMemzero(dnz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr); ierr = PetscMemzero(onz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr); for(ALECompat::Mesh::real_section_type::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) { const ALECompat::Mesh::real_section_type::atlas_type::point_type& point = *c_iter; if (globalOrder->isLocal(point)) { const ALE::Obj& adj = graph->cone(point); const ALECompat::Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(patch, point)[0]; const int row = rIdx.prefix; const int rSize = rIdx.index; for(ALECompat::Mesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) { const ALECompat::Mesh::real_section_type::atlas_type::point_type& neighbor = *v_iter; const ALECompat::Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(patch, neighbor)[0]; const int& cSize = cIdx.index; if (cSize > 0) { if (globalOrder->isLocal(neighbor)) { for(int r = 0; r < rSize; ++r) {dnz[row - firstRow + r] += cSize;} } else { for(int r = 0; r < rSize; ++r) {onz[row - firstRow + r] += cSize;} } } } } } if (topology->debug()) { int rank = topology->commRank(); for(int r = 0; r < numLocalRows; r++) { std::cout << "["<& topology, const ALE::Obj& atlas, const ALE::Obj& globalOrder, Mat A) { return preallocateOperatorCompat(topology, atlas, globalOrder, A); } #undef __FUNCT__ #define __FUNCT__ "updateOperatorCompat" PetscErrorCode updateOperatorCompat(Mat A, const ALE::Obj& section, const ALE::Obj& globalOrder, const ALECompat::Mesh::point_type& e, PetscScalar array[], InsertMode mode) { ALECompat::Mesh::real_section_type::patch_type patch = 0; static PetscInt indicesSize = 0; static PetscInt *indices = NULL; PetscInt numIndices = 0; PetscErrorCode ierr; PetscFunctionBegin; const ALE::Obj intervals = section->getIndices(patch, e, globalOrder); ierr = PetscLogEventBegin(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); if (section->debug()) {printf("[%d]mat for element %d\n", section->commRank(), e);} for(ALECompat::Mesh::real_section_type::IndexArray::iterator i_iter = intervals->begin(); i_iter != intervals->end(); ++i_iter) { numIndices += std::abs(i_iter->prefix); if (section->debug()) { printf("[%d]mat interval (%d, %d)\n", section->commRank(), i_iter->prefix, i_iter->index); } } if (indicesSize && (indicesSize != numIndices)) { ierr = PetscFree(indices); CHKERRQ(ierr); indices = NULL; } if (!indices) { indicesSize = numIndices; ierr = PetscMalloc(indicesSize * sizeof(PetscInt), &indices); CHKERRQ(ierr); } ierr = ExpandIntervals(intervals, indices); CHKERRQ(ierr); if (section->debug()) { for(int i = 0; i < numIndices; i++) { printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]); } for(int i = 0; i < numIndices; i++) { printf("[%d]", section->commRank()); for(int j = 0; j < numIndices; j++) { printf(" %g", array[i*numIndices+j]); } printf("\n"); } } ierr = MatSetValues(A, numIndices, indices, numIndices, indices, array, mode); if (ierr) { printf("[%d]ERROR in updateOperator: point %d\n", section->commRank(), e); for(int i = 0; i < numIndices; i++) { printf("[%d]mat indices[%d] = %d\n", section->commRank(), i, indices[i]); } CHKERRQ(ierr); } ierr = PetscLogEventEnd(Mesh_updateOperator,0,0,0,0);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCompatCreatePyLith" /*@C MeshCompatCreatePyLith - Create a Mesh from PyLith files. Not Collective Input Parameters: + dim - The topological mesh dimension . baseFilename - The basename for mesh files . zeroBase - Use 0 to start numbering - interpolate - The flag for mesh interpolation Output Parameter: . mesh - The Mesh object Notes: This is part of the PyLith 0.8 compatibility layer. DO NOT USE unless you are developing for that tool. Level: developer .keywords: mesh, PCICE .seealso: MeshCreate() @*/ PetscErrorCode MeshCompatCreatePyLith(MPI_Comm comm, const int dim, const char baseFilename[], PetscTruth zeroBase, PetscTruth interpolate, Mesh *mesh) { ALE::Obj m; PetscInt debug = 0; PetscTruth flag; PetscErrorCode ierr; PetscFunctionBegin; ierr = MeshCreate(comm, mesh);CHKERRQ(ierr); ierr = PetscOptionsGetInt(PETSC_NULL, "-debug", &debug, &flag);CHKERRQ(ierr); try { m = ALECompat::PyLith::Builder::readMesh(comm, dim, std::string(baseFilename), zeroBase, interpolate, debug); } catch(ALE::Exception e) { SETERRQ(PETSC_ERR_FILE_OPEN, e.message()); } ierr = MeshCompatSetMesh(*mesh, m);CHKERRQ(ierr); PetscFunctionReturn(0); }