#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: tri2d.c,v 1.38 2000/10/17 13:48:55 knepley Exp $"; #endif /* Implements 2d triangular grids */ #include "src/mesh/impls/triangular/2d/2dimpl.h" /*I "mesh.h" I*/ #include "tri2dView.h" #include "tri2dCSR.h" #ifdef PETSC_HAVE_TRIANGLE #include "tri2d_Triangle.h" #endif extern int PartitionGhostNodeIndex_Private(Partition, int, int *); EXTERN_C_BEGIN int MeshSerialize_Triangular_2D(MPI_Comm, Mesh *, PetscViewer, PetscTruth); EXTERN_C_END /*--------------------------------------------- Public Functions ----------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "MeshTriangular2DCalcAreas" /*@C MeshTriangular2DCalcAreas - Calculate and store the area of each element Collective on Mesh Input Parameters: + mesh - The mesh - create - The flag which creates the storage if it does not already exist Note: If areas have not been calculated before, and 'create' is PETSC_FALSE, then no action taken. Level: beginner .keywords: mesh, element, area .seealso: MeshTriangular2DCalcAspectRatios() @*/ int MeshTriangular2DCalcAreas(Mesh mesh, PetscTruth create) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int numFaces = mesh->numFaces; int numCorners = mesh->numCorners; int *faces = tri->faces; double *nodes = tri->nodes; double x21, x31, y21, y31; int elem; int ierr; PetscFunctionBegin; if (tri->areas || create == PETSC_TRUE) { if (tri->areas) { ierr = PetscFree(tri->areas); CHKERRQ(ierr); } ierr = PetscMalloc(numFaces * sizeof(double), &tri->areas); CHKERRQ(ierr); for(elem = 0; elem < numFaces; elem++) { if (mesh->isPeriodic == PETSC_TRUE) { x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]); y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]); x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]); y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]); } else { x21 = nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]; y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]; x31 = nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]; y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]; } tri->areas[elem] = PetscAbsReal(x21*y31 - x31*y21)/2.0; } PetscLogFlops(numFaces*8); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshTriangular2DCalcAspectRatios" /*@C MeshTriangular2DCalcAspectRatios - Calculate and store the aspect ratio of each element Collective on Mesh Input Parameters: + mesh - The mesh - create - The flag which creates the storage if it does not already exist Note: If aspect ratios have not been calculated before, and 'create' is PETSC_FALSE, then no action taken. Level: beginner .keywords: mesh, element, aspect ratio .seealso: MeshTriangular2DCalcAreas() @*/ int MeshTriangular2DCalcAspectRatios(Mesh mesh, PetscTruth create) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int numFaces = mesh->numFaces; int numCorners = mesh->numCorners; int *faces = tri->faces; double *nodes = tri->nodes; double x21, x31, x32, y21, y31, y32; double area; int elem; int ierr; PetscFunctionBegin; if (tri->aspectRatios || create == PETSC_TRUE) { if (tri->aspectRatios) { ierr = PetscFree(tri->aspectRatios); CHKERRQ(ierr); } ierr = PetscMalloc(numFaces * sizeof(double), &tri->aspectRatios); CHKERRQ(ierr); for(elem = 0; elem < numFaces; elem++) { if (mesh->isPeriodic == PETSC_TRUE) { x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]); y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]); x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]); y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]); x32 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners+1]*2]); y32 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]); } else { x21 = nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]; y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]; x31 = nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]; y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]; x32 = nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners+1]*2]; y32 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]; } area = PetscAbsReal(x21*y31 - x31*y21); tri->aspectRatios[elem] = PetscMax(x32*x32 + y32*y32, PetscMax(x21*x21 + y21*y21, x31*x31 + y31*y31)) / area; } PetscLogFlops(numFaces*19); } PetscFunctionReturn(0); } /*--------------------------------------------- Private Functions ---------------------------------------------------*/ #undef __FUNCT__ #define __FUNCT__ "MeshDestroy_Triangular_2D" static int MeshDestroy_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int ierr; PetscFunctionBegin; if (tri->nodes != PETSC_NULL) { ierr = PetscFree(tri->nodes); CHKERRQ(ierr); ierr = PetscFree(tri->markers); CHKERRQ(ierr); ierr = PetscFree(tri->degrees); CHKERRQ(ierr); } if (tri->nodesOld != PETSC_NULL) { ierr = PetscFree(tri->nodesOld); CHKERRQ(ierr); } if (tri->edges != PETSC_NULL) { ierr = PetscFree(tri->edges); CHKERRQ(ierr); ierr = PetscFree(tri->edgemarkers); CHKERRQ(ierr); } if (tri->faces != PETSC_NULL) { ierr = PetscFree(tri->faces); CHKERRQ(ierr); } if (tri->neighbors != PETSC_NULL) { ierr = PetscFree(tri->neighbors); CHKERRQ(ierr); } if (tri->bdNodes != PETSC_NULL) { ierr = PetscFree(tri->bdNodes); CHKERRQ(ierr); ierr = PetscFree(tri->bdEdges); CHKERRQ(ierr); ierr = PetscFree(tri->bdMarkers); CHKERRQ(ierr); ierr = PetscFree(tri->bdBegin); CHKERRQ(ierr); ierr = PetscFree(tri->bdEdgeBegin); CHKERRQ(ierr); } #if 0 if (tri->bdCtx != PETSC_NULL) { ierr = MeshBoundaryDestroy(tri->bdCtx); CHKERRQ(ierr); } if (tri->bdCtxNew != PETSC_NULL) { ierr = MeshBoundaryDestroy(tri->bdCtxNew); CHKERRQ(ierr); } #endif if (tri->areas != PETSC_NULL) { ierr = PetscFree(tri->areas); CHKERRQ(ierr); } if (tri->aspectRatios != PETSC_NULL) { ierr = PetscFree(tri->aspectRatios); CHKERRQ(ierr); } ierr = PetscFree(tri); CHKERRQ(ierr); if (mesh->support != PETSC_NULL) { ierr = PetscFree(mesh->support); CHKERRQ(ierr); } ierr = PartitionDestroy(mesh->part); CHKERRQ(ierr); if (mesh->nodeOrdering != PETSC_NULL) { ierr = AODestroy(mesh->nodeOrdering); CHKERRQ(ierr); } if (mesh->coarseMap != PETSC_NULL) { ierr = AODestroy(mesh->coarseMap); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshDebug_Triangular_2D" int MeshDebug_Triangular_2D(Mesh mesh, PetscTruth distributed) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; Partition p = mesh->part; Partition_Triangular_2D *q = PETSC_NULL; int numBd = mesh->numBd; int numElements = mesh->numFaces; int numCorners = mesh->numCorners; int numNodes = mesh->numNodes; int *elements = tri->faces; int *markers = tri->markers; int *bdMarkers = tri->bdMarkers; int degree; int *support; int *degrees; int *supports; int *sizes; int numProcs, rank; PetscTruth isMidnode; int proc, bd, elem, nElem, sElem, sElem2, corner, node, newNode, size; int ierr; PetscFunctionBegin; ierr = MPI_Comm_size(mesh->comm, &numProcs); CHKERRQ(ierr); ierr = MPI_Comm_rank(mesh->comm, &rank); CHKERRQ(ierr); if (p != PETSC_NULL) { q = (Partition_Triangular_2D *) p->data; } /* Check basic sizes */ ierr = PetscMalloc(numProcs * sizeof(int), &sizes); CHKERRQ(ierr); ierr = MPI_Allgather(&mesh->numBd, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of boundaries is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); #if 0 /* I made vertices distributed in linear meshes */ ierr = MPI_Allgather(&tri->numVertices, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of vertices is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); #endif ierr = MPI_Allgather(&mesh->numCorners, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of face corners is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); ierr = MPI_Allgather(&mesh->numBdEdges, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of boundary edges is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); if (distributed == PETSC_FALSE) { ierr = MPI_Allgather(&mesh->numNodes, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of nodes is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); ierr = MPI_Allgather(&mesh->numBdNodes, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of boundary nodes is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); ierr = MPI_Allgather(&mesh->numEdges, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of edges is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); ierr = MPI_Allgather(&mesh->numFaces, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) { if (sizes[proc] != sizes[proc+1]) { PetscPrintf(mesh->comm, "The number of faces is different on proc %d(%d) and proc %d(%d)\n", proc, sizes[proc], proc+1, sizes[proc+1]); ierr = 1; } } CHKERRQ(ierr); } else { ierr = MPI_Allreduce(&mesh->numNodes, &size, 1, MPI_INT, MPI_SUM, mesh->comm); CHKERRQ(ierr); if (size != q->numNodes) { SETERRQ2(PETSC_ERR_PLIB, "The total number of nodes %d should be %d\n", size, q->numNodes); } #if 0 ierr = MPI_Allreduce(&tri->numBdNodes, &size, 1, MPI_INT, MPI_SUM, mesh->comm); CHKERRQ(ierr); #else size = mesh->numBdNodes; #endif if (size != q->numBdNodes) { SETERRQ2(PETSC_ERR_PLIB, "The total number of boundary nodes %d should be %d\n", size, q->numBdNodes); } ierr = MPI_Allreduce(&mesh->numEdges, &size, 1, MPI_INT, MPI_SUM, mesh->comm); CHKERRQ(ierr); if (size != q->numEdges) { SETERRQ2(PETSC_ERR_PLIB, "The total number of edges %d should be %d\n", size, q->numEdges); } ierr = MPI_Allreduce(&mesh->numFaces, &size, 1, MPI_INT, MPI_SUM, mesh->comm); CHKERRQ(ierr); if (size != mesh->part->numElements) { SETERRQ2(PETSC_ERR_PLIB, "The total number of faces %d should be %d\n", size, mesh->part->numElements); } } ierr = PetscFree(sizes); CHKERRQ(ierr); if (distributed == PETSC_FALSE) { /* Check markers */ for(node = 0; node < numNodes; node++) { if (!markers[node]) continue; for(bd = 0; bd < numBd; bd++) if(bdMarkers[bd] == markers[node]) break; if (bd == numBd) { SETERRQ2(PETSC_ERR_PLIB, "The marker %d for node %d is invalid\n", markers[node], node); } } /* Check mesh connectivity */ for(node = 0; node < numNodes; node++) { ierr = MeshGetNodeSupport(mesh, node, 0, °ree, &support); CHKERRQ(ierr); /* Check node degree */ ierr = PetscMalloc(numProcs * sizeof(int), °rees); CHKERRQ(ierr); ierr = MPI_Allgather(°ree, 1, MPI_INT, degrees, 1, MPI_INT, mesh->comm); CHKERRQ(ierr); for(proc = 0, ierr = 0; proc < numProcs-1; proc++) if (degrees[proc] != degrees[proc+1]) { PetscPrintf(mesh->comm, "Degree of node %d is different on proc %d(%d) and proc %d(%d)\n", node, proc, degrees[proc], proc+1, degrees[proc+1]); PetscPrintf(mesh->comm, "Node Information:\n"); ierr = PetscSequentialPhaseBegin(mesh->comm, 1); CHKERRQ(ierr); PetscPrintf(PETSC_COMM_SELF, " Processor %d\n", rank); if ((rank == proc) || (rank == proc+1)) for(sElem = 0; sElem < degree; sElem++) PetscPrintf(PETSC_COMM_SELF, " Support : %d\n", support[sElem]); PetscPrintf(PETSC_COMM_SELF, " Marker : %d\n", tri->markers[node]); PetscPrintf(PETSC_COMM_SELF, " Location : (%g,%g)\n", tri->nodes[node*2], tri->nodes[node*2+1]); PetscPrintf(PETSC_COMM_SELF, " In Elements:"); for(elem = 0, isMidnode = PETSC_FALSE; elem < numElements; elem++) for(corner = 0; corner < numCorners; corner++) if (elements[elem*numCorners+corner] == node) { PetscPrintf(PETSC_COMM_SELF, " %d", elem); if (corner >= 3) isMidnode = PETSC_TRUE; } PetscPrintf(PETSC_COMM_SELF, "\n"); if (isMidnode == PETSC_TRUE) PetscPrintf(PETSC_COMM_SELF, " Midnode\n"); else PetscPrintf(PETSC_COMM_SELF, " Vertex\n"); ierr = PetscSequentialPhaseEnd(mesh->comm, 1); CHKERRQ(ierr); ierr = 1; } CHKERRQ(ierr); ierr = PetscFree(degrees); CHKERRQ(ierr); /* Check support elements */ ierr = PetscMalloc(degree*numProcs * sizeof(int), &supports); CHKERRQ(ierr); ierr = MPI_Allgather(support, degree, MPI_INT, supports, degree, MPI_INT, mesh->comm); CHKERRQ(ierr); for(sElem = 0, ierr = 0; sElem < degree; sElem++) { nElem = supports[sElem]; for(proc = 1; proc < numProcs; proc++) { for(sElem2 = 0; sElem2 < degree; sElem2++) if (supports[proc*degree+sElem2] == nElem) break; if (sElem2 == degree) { PetscPrintf(mesh->comm, "Support element %d of node %d on proc %d(%d) is not present on proc %d\n", sElem, node, 0, supports[sElem], proc); PetscPrintf(mesh->comm, " Support of node %d on proc %d:\n ", node, proc); for(sElem2 = 0; sElem2 < degree; sElem2++) PetscPrintf(mesh->comm, " %d", supports[proc*degree+sElem2]); PetscPrintf(mesh->comm, "\n"); ierr = 1; } } } CHKERRQ(ierr); ierr = PetscFree(supports); CHKERRQ(ierr); /* Check that node only appears inside elements in the support */ for(elem = 0, ierr = 0; elem < numElements; elem++) for(corner = 0; corner < numCorners; corner++) { newNode = elements[elem*numCorners+corner]; if (node == newNode) { for(sElem = 0; sElem < degree; sElem++) if (support[sElem] == elem) break; if (sElem == degree) { PetscPrintf(mesh->comm, "Node %d found in element %d which is not present in the support\n", node, elem); ierr = 1; } } } CHKERRQ(ierr); ierr = MeshRestoreNodeSupport(mesh, node, 0, °ree, &support); CHKERRQ(ierr); } } else { /* Check markers */ for(node = 0; node < q->numOverlapNodes; node++) { if (!markers[node]) continue; for(bd = 0; bd < numBd; bd++) if(bdMarkers[bd] == markers[node]) break; if (bd == numBd) { SETERRQ2(PETSC_ERR_PLIB, "The marker %d for node %d is invalid\n", markers[node], node); } } /* Check mesh connectivity */ for(node = 0; node < q->numLocNodes; node++) { ierr = MeshGetNodeSupport(mesh, node, 0, °ree, &support); CHKERRQ(ierr); /* Check that node only appears inside elements in the support */ for(elem = 0, ierr = 0; elem < p->numOverlapElements; elem++) { for(corner = 0; corner < numCorners; corner++) { newNode = elements[elem*numCorners+corner]; if (node == newNode) { for(sElem = 0; sElem < degree; sElem++) { if (support[sElem] == elem) break; } if (sElem == degree) { PetscPrintf(PETSC_COMM_SELF, "[%d]Node %d found in element %d which is not present in the support\n", p->rank, node, elem); ierr = 1; } } } } CHKERRQ(ierr); ierr = MeshRestoreNodeSupport(mesh, node, 0, °ree, &support); CHKERRQ(ierr); } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetupSupport_Triangular_2D" int MeshSetupSupport_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int edge, node; int ierr; PetscFunctionBegin; /* Calculate maximum degree of vertices */ if (mesh->numNodes > 0) { ierr = PetscMalloc(mesh->numNodes * sizeof(int), &tri->degrees); CHKERRQ(ierr); } ierr = PetscMemzero(tri->degrees, mesh->numNodes * sizeof(int)); CHKERRQ(ierr); for(edge = 0; edge < mesh->numEdges*2; edge++) { tri->degrees[tri->edges[edge]]++; } for(node = 0, mesh->maxDegree = 0; node < mesh->numNodes; node++) { mesh->maxDegree = PetscMax(mesh->maxDegree, tri->degrees[node]); } if (mesh->maxDegree > 0) { ierr = PetscMalloc(mesh->maxDegree * sizeof(int), &mesh->support); CHKERRQ(ierr); } mesh->supportTaken = PETSC_FALSE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCheckBoundary_Triangular_2D" int MeshCheckBoundary_Triangular_2D(Mesh mesh) { MeshBoundary2D *bdCtx = mesh->bdCtx; int *markers = bdCtx->markers; int *segments = bdCtx->segments; int *segMarkers = bdCtx->segMarkers; int numBd = 0; int numBdVertices = 0; int numBdSegments = 0; int *bdMarkers; int *bdBegin; int *bdSegmentBegin; int bd, vertex, segment, marker, rank; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mesh, MESH_COOKIE); ierr = MPI_Comm_rank(mesh->comm, &rank); CHKERRQ(ierr); if (rank != 0) PetscFunctionReturn(0); ierr = PetscMalloc(bdCtx->numBd * sizeof(int), &bdMarkers); CHKERRQ(ierr); ierr = PetscMalloc((bdCtx->numBd+1) * sizeof(int), &bdBegin); CHKERRQ(ierr); ierr = PetscMalloc((bdCtx->numBd+1) * sizeof(int), &bdSegmentBegin); CHKERRQ(ierr); ierr = PetscMemzero(bdBegin, (bdCtx->numBd+1) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(bdSegmentBegin, (bdCtx->numBd+1) * sizeof(int)); CHKERRQ(ierr); for(vertex = 0; vertex < bdCtx->numVertices; vertex++) { if (markers[vertex] == 0) continue; numBdVertices++; /* Check for new marker */ for(bd = 0; bd < numBd; bd++) { if (markers[vertex] == bdMarkers[bd]) break; } if (bd == numBd) { if (numBd >= bdCtx->numBd) SETERRQ1(PETSC_ERR_ARG_CORRUPT, "More markers present than declared: %d", bdCtx->numBd); /* Insert new marker */ for(bd = 0; bd < numBd; bd++) { if (markers[vertex] < bdMarkers[bd]) break; } if (bd < numBd) { ierr = PetscMemmove(&bdMarkers[bd+1], &bdMarkers[bd], (numBd-bd) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemmove(&bdBegin[bd+2], &bdBegin[bd+1], (numBd-bd) * sizeof(int)); CHKERRQ(ierr); bdBegin[bd+1] = 0; bdSegmentBegin[bd+1] = 0; } bdMarkers[bd] = markers[vertex]; numBd++; bdBegin[bd+1]++; } else { bdBegin[bd+1]++; } } for(segment = 0; segment < bdCtx->numSegments; segment++) { if (segMarkers[segment] == 0) continue; marker = segMarkers[segment]; for(bd = 0; bd < numBd; bd++) { if (marker == bdMarkers[bd]) break; } if (bd == numBd) SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid segment marker %d on segment %d", marker, segment); numBdSegments++; bdSegmentBegin[bd+1]++; if (markers[segments[segment*2]] != marker) { SETERRQ4(PETSC_ERR_PLIB, "Marker %d on left vertex %d does not match marker %d on its segment %d", markers[segments[segment*2]], segments[segment*2], marker, segment); } if (markers[segments[segment*2+1]] != marker) { SETERRQ4(PETSC_ERR_PLIB, "Marker %d on right vertex %d does not match marker %d on its segment %d", markers[segments[segment*2+1]], segments[segment*2+1], marker, segment); } } /* Do prefix sums to get position offsets */ for(bd = 2; bd <= numBd; bd++) { bdBegin[bd] = bdBegin[bd-1] + bdBegin[bd]; bdSegmentBegin[bd] = bdSegmentBegin[bd-1] + bdSegmentBegin[bd]; } if (numBd != bdCtx->numBd) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundaries %d should be %d", numBd, bdCtx->numBd); } if (bdBegin[numBd] != numBdVertices) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary vertices %d should be %d", bdBegin[numBd], numBdVertices); } if (bdSegmentBegin[numBd] != numBdSegments) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary segments %d should be %d", bdSegmentBegin[numBd], numBdSegments); } ierr = PetscFree(bdMarkers); CHKERRQ(ierr); ierr = PetscFree(bdBegin); CHKERRQ(ierr); ierr = PetscFree(bdSegmentBegin); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetupBoundary_Triangular_2D" int MeshSetupBoundary_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int numBd = mesh->numBd; int numNodes = mesh->numNodes; int *markers = tri->markers; int numEdges = mesh->numEdges; int *edges = tri->edges; int *edgemarkers = tri->edgemarkers; int numFaces = mesh->numFaces; int numCorners = mesh->numCorners; int *faces = tri->faces; int newNumBd; /* Current number of different boundary markers */ int numBdEdges; /* Current offset into bdEdges[] */ int *bdNodes; /* Node numbers for boundary nodes ordered by boundary */ int *bdEdges; /* Node numbers for boundary edges ordered by boundary */ int *bdBeginOff; /* Current offset into the bdNodes or bdEdges array */ int *bdSeen; /* Flags for boundaries already processed */ PetscTruth closed; /* Indicates whether a boundary has been closed */ int degree; int *support; int rank, bd, marker, node, nextNode, midnode, bdNode, nextBdNode, midBdNode, edge, bdEdge; int elem, cElem, supElem, corner, supCorner, tmp; int ierr; PetscFunctionBegin; ierr = MPI_Comm_rank(mesh->comm, &rank); CHKERRQ(ierr); ierr = PetscLogEventBegin(MESH_SetupBoundary, mesh, 0, 0, 0); CHKERRQ(ierr); tri->areas = PETSC_NULL; tri->aspectRatios = PETSC_NULL; mesh->numBdNodes = 0; mesh->numBdEdges = 0; if (numBd == 0) { tri->bdMarkers = PETSC_NULL; tri->bdBegin = PETSC_NULL; tri->bdEdgeBegin = PETSC_NULL; tri->bdNodes = PETSC_NULL; tri->bdEdges = PETSC_NULL; PetscFunctionReturn(0); } ierr = PetscMalloc(numBd * sizeof(int), &tri->bdMarkers); CHKERRQ(ierr); ierr = PetscMalloc((numBd+1) * sizeof(int), &tri->bdBegin); CHKERRQ(ierr); ierr = PetscMalloc((numBd+1) * sizeof(int), &tri->bdEdgeBegin); CHKERRQ(ierr); ierr = PetscMalloc((numBd+1) * sizeof(int), &bdBeginOff); CHKERRQ(ierr); PetscLogObjectMemory(mesh, (numBd + (numBd+1)*2) * sizeof(int)); ierr = PetscMemzero(tri->bdMarkers, numBd * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(tri->bdBegin, (numBd+1) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(tri->bdEdgeBegin, (numBd+1) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemzero(bdBeginOff, (numBd+1) * sizeof(int)); CHKERRQ(ierr); for(node = 0, newNumBd = 0; node < numNodes; node++) { /* Get number of boundary nodes and markers */ if (markers[node]) { mesh->numBdNodes++; /* Check for new marker */ for(bd = 0; bd < newNumBd; bd++) { if (markers[node] == tri->bdMarkers[bd]) break; } if (bd == newNumBd) { /* Insert new marker */ for(bd = 0; bd < newNumBd; bd++) { if (markers[node] < tri->bdMarkers[bd]) break; } if (bd < newNumBd) { ierr = PetscMemmove(&tri->bdMarkers[bd+1], &tri->bdMarkers[bd], (newNumBd-bd) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemmove(&tri->bdBegin[bd+2], &tri->bdBegin[bd+1], (newNumBd-bd) * sizeof(int)); CHKERRQ(ierr); tri->bdBegin[bd+1] = 0; } tri->bdMarkers[bd] = markers[node]; newNumBd++; tri->bdBegin[bd+1]++; } else { tri->bdBegin[bd+1]++; } } } /* Do prefix sums to get position offsets */ for(bd = 2; bd <= numBd; bd++) { tri->bdBegin[bd] = tri->bdBegin[bd-1] + tri->bdBegin[bd]; } /* This is shameful -- I will write the code to look over edges soon */ if (numCorners == 3) { ierr = PetscMemcpy(tri->bdEdgeBegin, tri->bdBegin, (numBd+1) * sizeof(int)); CHKERRQ(ierr); } else if (numCorners == 6) { /* We know that there are an even number of nodes in every boundary */ for(bd = 0; bd <= numBd; bd++) { if (tri->bdBegin[bd]%2) { SETERRQ2(PETSC_ERR_PLIB, "There are %d nodes on boundary %d (not divisible by two)", tri->bdBegin[bd]-tri->bdBegin[bd-1], tri->bdMarkers[bd]); } else { tri->bdEdgeBegin[bd] = tri->bdBegin[bd]/2; } } } else { SETERRQ1(PETSC_ERR_SUP, "Number of local nodes %d not supported", numCorners); } /* Get number of boundary edges */ for(edge = 0; edge < numEdges; edge++) { marker = edgemarkers[edge]; if (marker) { mesh->numBdEdges++; ierr = MeshGetMidnodeFromEdge(mesh, edge, &midnode); CHKERRQ(ierr); if (markers[edges[edge*2]] != marker) { if ((markers[edges[edge*2+1]] == marker) && ((markers[midnode] == markers[edges[edge*2]]) || (markers[midnode] == markers[edges[edge*2+1]]))) { /* I assume here that the generator mistakenly included an edge between two boundaries */ PetscPrintf(PETSC_COMM_SELF, "[%d]Removing edge %d between boundaries %d and %d from boundary\n", rank, edge, markers[edges[edge*2]], markers[edges[edge*2+1]]); edgemarkers[edge] = 0; markers[midnode] = 0; continue; } else { SETERRQ5(PETSC_ERR_PLIB, "Marker %d on left node %d does not match marker %d on its edge %d, right marker is %d", markers[edges[edge*2]], edges[edge*2], marker, edge, markers[edges[edge*2+1]]); } } if (markers[edges[edge*2+1]] != marker) { if ((markers[edges[edge*2]] == marker) && ((markers[midnode] == markers[edges[edge*2]]) || (markers[midnode] == markers[edges[edge*2+1]]))) { /* I assume here that the generator mistakenly included an edge between two boundaries */ PetscPrintf(PETSC_COMM_SELF, "[%d]Removing edge %d between boundaries %d and %d from boundary\n", rank, edge, markers[edges[edge*2]], markers[edges[edge*2+1]]); edgemarkers[edge] = 0; markers[midnode] = 0; continue; } else { SETERRQ5(PETSC_ERR_PLIB, "Marker %d on right node %d does not match marker %d on its edge %d, left marker is %d", markers[edges[edge*2+1]], edges[edge*2+1], marker, edge, markers[edges[edge*2]]); } } if (markers[midnode] != marker) { SETERRQ5(PETSC_ERR_PLIB, "Marker %d on midnode %d does not match marker %d on its edge %d, left marker is %d", markers[midnode], midnode, marker, edge, markers[edges[edge*2]]); } } } /* Check boundary information consistency */ if (newNumBd != numBd) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundaries %d should be %d", newNumBd, numBd); } if (tri->bdBegin[numBd] != mesh->numBdNodes) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary nodes %d should be %d", tri->bdBegin[numBd], mesh->numBdNodes); } if (tri->bdEdgeBegin[numBd] != mesh->numBdEdges) { SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary edges %d should be %d", tri->bdEdgeBegin[numBd], mesh->numBdEdges); } ierr = PetscMalloc(mesh->numBdNodes * sizeof(int), &bdNodes); CHKERRQ(ierr); ierr = PetscMalloc(mesh->numBdEdges * sizeof(int), &bdEdges); CHKERRQ(ierr); PetscLogObjectMemory(mesh, (mesh->numBdNodes + mesh->numBdEdges) * sizeof(int)); /* Split nodes by marker */ ierr = PetscMemcpy(bdBeginOff, tri->bdBegin, (numBd+1) * sizeof(int)); CHKERRQ(ierr); for(node = 0; node < numNodes; node++) { for(bd = 0; bd < numBd; bd++) { if (markers[node] == tri->bdMarkers[bd]) { #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(PETSC_COMM_WORLD, "bd: %d bdNode[%d] = %d\n", bd, bdBeginOff[bd], node); #endif bdNodes[bdBeginOff[bd]++] = node; } } } for(bd = 0; bd < numBd; bd++) { if (tri->bdBegin[bd+1] != bdBeginOff[bd]) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node marker information"); } /* Split edges by marker */ ierr = PetscMemcpy(bdBeginOff, tri->bdEdgeBegin, (numBd+1) * sizeof(int)); CHKERRQ(ierr); for(edge = 0; edge < numEdges; edge++) { for(bd = 0; bd < numBd; bd++) { if (edgemarkers[edge] == tri->bdMarkers[bd]) { #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(PETSC_COMM_WORLD, "bd: %d bdEdge[%d] = %d\n", bd, bdBeginOff[bd], edge); #endif bdEdges[bdBeginOff[bd]++] = edge; } } } for(bd = 0; bd < numBd; bd++) { if (tri->bdEdgeBegin[bd+1] != bdBeginOff[bd]) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge marker information"); } ierr = PetscFree(bdBeginOff); CHKERRQ(ierr); /* Order boundary counter-clockwise */ ierr = PetscMalloc(numBd * sizeof(int), &bdSeen); CHKERRQ(ierr); ierr = PetscMemzero(bdSeen, numBd * sizeof(int)); /* Loop over elements */ for(elem = 0; elem < numFaces; elem++) { /* Find an element with a node on the boundary */ for(corner = 0; corner < 3; corner++) { if (markers[faces[elem*numCorners+corner]]) { /* Find boundary index */ cElem = elem; node = faces[elem*numCorners+corner]; for(bd = 0; bd < numBd; bd++) if (markers[node] == tri->bdMarkers[bd]) break; if (bd == numBd) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary marker"); /* Check for processed boundaries */ if (bdSeen[bd]) continue; numBdEdges = tri->bdEdgeBegin[bd]; /* Start building this boundary */ #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, "Starting boundary %d beginning at %d ending before %d\n", bd, tri->bdBegin[bd], tri->bdBegin[bd+1]); #endif /* Find initial node */ for(bdNode = tri->bdBegin[bd]; bdNode < tri->bdBegin[bd+1]; bdNode++) { if (bdNodes[bdNode] == node) { /* Move this node to the head of the list */ #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " moving node %d from %d to %d\n", bdNodes[bdNode], bdNode, tri->bdBegin[bd]); #endif tmp = bdNodes[tri->bdBegin[bd]]; bdNodes[tri->bdBegin[bd]] = bdNodes[bdNode]; bdNodes[bdNode] = tmp; break; } } /* Order edges counterclockwise around a boundary */ /* I do not currently check the orientation of the constructed boundary */ for(bdNode = tri->bdBegin[bd], closed = PETSC_FALSE; bdNode < tri->bdBegin[bd+1]; bdNode++) { #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, "\n At boundary node %d x:%lf y: %lf\n", bdNode, tri->nodes[bdNodes[bdNode]*2], tri->nodes[bdNodes[bdNode]*2+1]); #ifdef MESH_TRACE_BOUNDARY_CREATE_DETAIL for(node = tri->bdBegin[bd]; node < tri->bdBegin[bd+1]; node++) PetscPrintf(mesh->comm, " bdNode[%d]: %d\n", node, bdNodes[node]); #endif #endif /* Find a neighbor of the point -- Could maybe do better than linear search */ for(bdEdge = numBdEdges; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) { edge = bdEdges[bdEdge]; if ((edgemarkers[edge] == tri->bdMarkers[bd]) && (((edges[edge*2] == bdNodes[bdNode]) && (markers[edges[edge*2+1]] == tri->bdMarkers[bd])) || ((edges[edge*2+1] == bdNodes[bdNode]) && (markers[edges[edge*2]] == tri->bdMarkers[bd])))) { /* Get neighboring node number */ if (edges[edge*2] == bdNodes[bdNode]) { node = edges[edge*2+1]; } else { node = edges[edge*2]; } /* Find neighboring node in bdNodes[] */ for(nextBdNode = tri->bdBegin[bd]; nextBdNode < tri->bdBegin[bd+1]; nextBdNode++) if (bdNodes[nextBdNode] == node) break; #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " found connection along edge %d to bd node %d = node %d\n", edge, nextBdNode, node); #endif if (nextBdNode > bdNode) { /* Insert midnode */ if (numCorners == 6) { /* Move this node next to the connected one */ #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " moving node %d from %d to %d\n", bdNodes[nextBdNode], nextBdNode, bdNode+2); #endif tmp = bdNodes[bdNode+2]; bdNodes[bdNode+2] = bdNodes[nextBdNode]; bdNodes[nextBdNode] = tmp; nextBdNode = bdNode+2; /* Walk around the node looking for a boundary edge and reset containing element */ node = bdNodes[bdNode]; nextNode = bdNodes[nextBdNode]; ierr = MeshGetNodeSupport(mesh, node, cElem, °ree, &support); CHKERRQ(ierr); for(supElem = 0, midnode = -1; supElem < degree; supElem++) { for(supCorner = 0; supCorner < 3; supCorner++) { cElem = support[supElem]; if ((faces[cElem*numCorners+supCorner] == node) && (faces[cElem*numCorners+((supCorner+1)%3)] == nextNode)) { midnode = faces[cElem*numCorners+((((supCorner*2+1)%3)*2)%3 + 3)]; supElem = degree; break; } else if ((faces[cElem*numCorners+supCorner] == node) && (faces[cElem*numCorners+((supCorner+2)%3)] == nextNode)) { midnode = faces[cElem*numCorners+((((supCorner*2+2)%3)*2)%3 + 3)]; supElem = degree; break; } } } ierr = MeshRestoreNodeSupport(mesh, node, cElem, °ree, &support); CHKERRQ(ierr); /* Find midnode in bdNodes[] */ for(midBdNode = tri->bdBegin[bd]; midBdNode < tri->bdBegin[bd+1]; midBdNode++) if (bdNodes[midBdNode] == midnode) break; if (midBdNode == tri->bdBegin[bd+1]) SETERRQ(PETSC_ERR_PLIB, "Unable to locate midnode on boundary"); #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " moving midnode %d in elem %d from %d to %d\n", midnode, cElem, midBdNode, bdNode+1); #endif tmp = bdNodes[bdNode+1]; bdNodes[bdNode+1] = bdNodes[midBdNode]; bdNodes[midBdNode] = tmp; bdNode++; } else { /* numCorners == 3 */ /* Move this node next to the connected one */ #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " moving node %d from %d to %d\n", bdNodes[nextBdNode], nextBdNode, bdNode+1); #endif tmp = bdNodes[bdNode+1]; bdNodes[bdNode+1] = bdNodes[nextBdNode]; bdNodes[nextBdNode] = tmp; } /* Reorder bdEdges[] */ #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " moving edge %d from %d to %d\n", edge, bdEdge, numBdEdges); #endif tmp = bdEdges[numBdEdges]; bdEdges[numBdEdges] = bdEdges[bdEdge]; bdEdges[bdEdge] = tmp; numBdEdges++; break; } /* Check that first and last nodes are connected (closed boundary) */ else if ((nextBdNode == tri->bdBegin[bd]) && (bdNode != (numCorners == 6 ? nextBdNode+2 : nextBdNode+1))) { if (bdEdges[numBdEdges++] != edge) SETERRQ(PETSC_ERR_PLIB, "Invalid closing edge"); closed = PETSC_TRUE; /* End the loop since only a midnode is left */ bdNode = tri->bdBegin[bd+1]; #ifdef MESH_TRACE_BOUNDARY_CREATE PetscPrintf(mesh->comm, " adding edge %d total: %d\n", edge, numBdEdges); #endif break; } } } if (bdEdge == tri->bdEdgeBegin[bd+1]) SETERRQ(PETSC_ERR_PLIB, "No connection"); } /* Check for closed boundary */ if (closed == PETSC_FALSE) { PetscPrintf(PETSC_COMM_SELF, "Boundary %d with marker %d is not closed\n", bd, tri->bdMarkers[bd]); } /* Check size of boundary */ if (tri->bdBegin[bd+1] != numBdEdges*(numCorners == 3 ? 1 : 2)) { PetscPrintf(PETSC_COMM_SELF, "On boundary %d with marker %d, edge count %d x %d should equal %d (was %d)\n", bd, tri->bdMarkers[bd], numBdEdges, (numCorners == 3 ? 1 : 2), tri->bdBegin[bd+1], tri->bdBegin[bd]); SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge marker information"); } #ifdef MESH_TRACE_BOUNDARY_CREATE /* Check boundary nodes */ for(bdNode = tri->bdBegin[bd]; bdNode < tri->bdBegin[bd+1]; bdNode++) PetscPrintf(mesh->comm, "bd: %4d bdNode: %4d node: %5d\n", bd, bdNode, bdNodes[bdNode], bdEdges[bdNode], edges[bdEdges[bdNode]*2], edges[bdEdges[bdNode]*2+1]); /* Check boundary edges */ for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) PetscPrintf(mesh->comm, "bd: %4d bdEdge: %4d edge: %5d start: %5d end: %5d\n", bd, bdEdge, bdEdges[bdEdge], edges[bdEdges[bdEdge]*2], edges[bdEdges[bdEdge]*2+1]); #endif bdSeen[bd] = 1; } } } ierr = PetscFree(bdSeen); CHKERRQ(ierr); /* Set fields */ tri->bdNodes = bdNodes; tri->bdEdges = bdEdges; /* Diagnostics */ #ifdef MESH_TRACE_BOUNDARY_REORDER int minNode, minBdNode, minEdge, minBdEdge, dir, nextBdEdge; ierr = PetscMalloc(mesh->numBdNodes * sizeof(int), &bdNodes); CHKERRQ(ierr); ierr = PetscMalloc(mesh->numBdEdges * sizeof(int), &bdEdges); CHKERRQ(ierr); for(bd = 0; bd < numBd; bd++) { /* Find smallest node */ for(bdNode = tri->bdBegin[bd], minNode = tri->numNodes, minBdNode = -1; bdNode < tri->bdBegin[bd+1]; bdNode++) if (tri->bdNodes[bdNode] < minNode) { minNode = tri->bdNodes[bdNode]; minBdNode = bdNode; } /* Proceed in the direction of the smaller node */ if (minBdNode == tri->bdBegin[bd]) { if (tri->bdNodes[minBdNode+1] < tri->bdNodes[tri->bdBegin[bd+1]-1]) dir = 1; else dir = 0; } else if (minBdNode == tri->bdBegin[bd+1]-1) { if (tri->bdNodes[tri->bdBegin[bd]] < tri->bdNodes[minBdNode-1]) dir = 0; else dir = 1; } else if (tri->bdNodes[minBdNode+1] < tri->bdNodes[minBdNode-1]) { dir = 1; } else { dir = 0; } if (dir) { for(nextBdNode = tri->bdBegin[bd], bdNode = minBdNode; bdNode < tri->bdBegin[bd+1]; nextBdNode++, bdNode++) bdNodes[nextBdNode] = tri->bdNodes[bdNode]; for(bdNode = tri->bdBegin[bd]; bdNode < minBdNode; nextBdNode++, bdNode++) bdNodes[nextBdNode] = tri->bdNodes[bdNode]; } else { for(nextBdNode = tri->bdBegin[bd], bdNode = minBdNode; bdNode >= tri->bdBegin[bd]; nextBdNode++, bdNode--) bdNodes[nextBdNode] = tri->bdNodes[bdNode]; for(bdNode = tri->bdBegin[bd+1]-1; bdNode > minBdNode; nextBdNode++, bdNode--) bdNodes[nextBdNode] = tri->bdNodes[bdNode]; } } for(bd = 0; bd < numBd; bd++) { /* Find smallest edge */ for(bdEdge = tri->bdEdgeBegin[bd], minEdge = tri->numEdges, minBdEdge = -1; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) if (tri->bdEdges[bdEdge] < minEdge) { minEdge = tri->bdEdges[bdEdge]; minBdEdge = bdEdge; } /* Proceed in the direction of the smaller edge */ if (minBdEdge == tri->bdEdgeBegin[bd]) { if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[tri->bdEdgeBegin[bd+1]-1]) dir = 1; else dir = 0; } else if (minBdEdge == tri->bdEdgeBegin[bd+1]-1) { if (tri->bdEdges[tri->bdEdgeBegin[bd]] < tri->bdEdges[minBdEdge-1]) dir = 0; else dir = 1; } else if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[minBdEdge-1]) { dir = 1; } else { dir = 0; } if (dir) { for(nextBdEdge = tri->bdEdgeBegin[bd], bdEdge = minBdEdge; bdEdge < tri->bdEdgeBegin[bd+1]; nextBdEdge++, bdEdge++) bdEdges[nextBdEdge] = tri->bdEdges[bdEdge]; for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < minBdEdge; nextBdEdge++, bdEdge++) bdEdges[nextBdEdge] = tri->bdEdges[bdEdge]; } else { for(nextBdEdge = tri->bdEdgeBegin[bd], bdEdge = minBdEdge; bdEdge >= tri->bdEdgeBegin[bd]; nextBdEdge++, bdEdge--) bdEdges[nextBdEdge] = tri->bdEdges[bdEdge]; for(bdEdge = tri->bdEdgeBegin[bd+1]-1; bdEdge > minBdEdge; nextBdEdge++, bdEdge--) bdEdges[nextBdEdge] = tri->bdEdges[bdEdge]; } } ierr = PetscMemcpy(tri->bdNodes, bdNodes, mesh->numBdNodes * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(tri->bdEdges, bdEdges, mesh->numBdEdges * sizeof(int)); CHKERRQ(ierr); ierr = PetscFree(bdNodes); CHKERRQ(ierr); ierr = PetscFree(bdEdges); CHKERRQ(ierr); #endif #ifdef MESH_TRACE_BOUNDARY int minNode, minBdNode, minEdge, minBdEdge, dir; for(bd = 0; bd < numBd; bd++) { /* Find smallest node */ for(bdNode = tri->bdBegin[bd], minNode = tri->numNodes, minBdNode = -1; bdNode < tri->bdBegin[bd+1]; bdNode++) if (tri->bdNodes[bdNode] < minNode) { minNode = tri->bdNodes[bdNode]; minBdNode = bdNode; } /* Proceed in the direction of the smaller node */ if (minBdNode == tri->bdBegin[bd]) { if (tri->bdNodes[minBdNode+1] < tri->bdNodes[tri->bdBegin[bd+1]-1]) dir = 1; else dir = 0; } else if (minBdNode == tri->bdBegin[bd+1]-1) { if (tri->bdNodes[tri->bdBegin[bd]] < tri->bdNodes[minBdNode-1]) dir = 0; else dir = 1; } else if (tri->bdNodes[minBdNode+1] < tri->bdNodes[minBdNode-1]) { dir = 1; } else { dir = 0; } if (dir) { for(bdNode = minBdNode; bdNode < tri->bdBegin[bd+1]; bdNode++) PetscPrintf(mesh->comm, "bd: %4d node: %5d\n", bd, tri->bdNodes[bdNode]); for(bdNode = tri->bdBegin[bd]; bdNode < minBdNode; bdNode++) PetscPrintf(mesh->comm, "bd: %4d node: %5d\n", bd, tri->bdNodes[bdNode]); } else { for(bdNode = minBdNode; bdNode >= tri->bdBegin[bd]; bdNode--) PetscPrintf(mesh->comm, "bd: %4d node: %5d\n", bd, tri->bdNodes[bdNode]); for(bdNode = tri->bdBegin[bd+1]-1; bdNode > minBdNode; bdNode--) PetscPrintf(mesh->comm, "bd: %4d node: %5d\n", bd, tri->bdNodes[bdNode]); } } for(bd = 0; bd < numBd; bd++) { /* Find smallest edge */ for(bdEdge = tri->bdEdgeBegin[bd], minEdge = tri->numEdges, minBdEdge = -1; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) if (tri->bdEdges[bdEdge] < minEdge) { minEdge = tri->bdEdges[bdEdge]; minBdEdge = bdEdge; } /* Proceed in the direction of the smaller edge */ if (minBdEdge == tri->bdEdgeBegin[bd]) { if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[tri->bdEdgeBegin[bd+1]-1]) dir = 1; else dir = 0; } else if (minBdEdge == tri->bdEdgeBegin[bd+1]-1) { if (tri->bdEdges[tri->bdEdgeBegin[bd]] < tri->bdEdges[minBdEdge-1]) dir = 0; else dir = 1; } else if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[minBdEdge-1]) { dir = 1; } else { dir = 0; } if (dir) { for(bdEdge = minBdEdge; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5d\n", bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]); for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < minBdEdge; bdEdge++) PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5d\n", bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]); } else { for(bdEdge = minBdEdge; bdEdge >= tri->bdEdgeBegin[bd]; bdEdge--) PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5d\n", bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]); for(bdEdge = tri->bdEdgeBegin[bd+1]-1; bdEdge > minBdEdge; bdEdge--) PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5d\n", bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]); } } #endif ierr = PetscLogEventEnd(MESH_SetupBoundary, mesh, 0, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshAssemble_Private" /* MeshAssemble_Private - This function assembles the mesh entirely on the first processor. Collective on Mesh Input Parameters: . mesh - The mesh being refined Output Parameter: + nodes - The node coordinates . markers - The node markers . faces - The nodes for each element - edges - The nodes for each edge Level: developer .keywords mesh, assemble .seealso MeshInitRefineInput_Triangle() */ int MeshAssemble_Private(Mesh mesh, double **nodes, int **markers, int **faces, int **edges) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; Partition p = mesh->part; Partition_Triangular_2D *q = (Partition_Triangular_2D *) p->data; int numNodes = q->numNodes; int numFaces = p->numElements; int numEdges = q->numEdges; int numCorners = mesh->numCorners; int numProcs = p->numProcs; int rank = p->rank; int *locFaces = PETSC_NULL; int *locEdges = PETSC_NULL; int *numRecvNodes, *numRecvMarkers, *numRecvFaces, *numRecvEdges; int *nodeOffsets, *markerOffsets, *faceOffsets, *edgeOffsets; int proc, elem, edge, size; int ierr; PetscFunctionBegin; /* Allocate global arrays */ if (rank == 0) { ierr = PetscMalloc(numNodes*2 * sizeof(double), nodes); CHKERRQ(ierr); ierr = PetscMalloc(numNodes * sizeof(int), markers); CHKERRQ(ierr); ierr = PetscMalloc(numFaces*numCorners * sizeof(int), faces); CHKERRQ(ierr); ierr = PetscMalloc(numEdges*2 * sizeof(int), edges); CHKERRQ(ierr); } if (numProcs > 1) { if (mesh->numFaces > 0) { ierr = PetscMalloc(mesh->numFaces*numCorners * sizeof(int), &locFaces); CHKERRQ(ierr); } if (mesh->numEdges > 0) { ierr = PetscMalloc(mesh->numEdges*2 * sizeof(int), &locEdges); CHKERRQ(ierr); } /* Calculate offsets */ ierr = PetscMalloc(numProcs * sizeof(int), &numRecvNodes); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &numRecvMarkers); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &numRecvEdges); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &numRecvFaces); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &nodeOffsets); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &markerOffsets); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &edgeOffsets); CHKERRQ(ierr); ierr = PetscMalloc(numProcs * sizeof(int), &faceOffsets); CHKERRQ(ierr); for(proc = 0; proc < numProcs; proc++) { numRecvNodes[proc] = (q->firstNode[proc+1] - q->firstNode[proc])*2; numRecvMarkers[proc] = (q->firstNode[proc+1] - q->firstNode[proc]); numRecvEdges[proc] = (q->firstEdge[proc+1] - q->firstEdge[proc])*2; numRecvFaces[proc] = (p->firstElement[proc+1] - p->firstElement[proc])*numCorners; } nodeOffsets[0] = 0; markerOffsets[0] = 0; edgeOffsets[0] = 0; faceOffsets[0] = 0; for(proc = 1; proc < numProcs; proc++) { nodeOffsets[proc] = numRecvNodes[proc-1] + nodeOffsets[proc-1]; markerOffsets[proc] = numRecvMarkers[proc-1] + markerOffsets[proc-1]; edgeOffsets[proc] = numRecvEdges[proc-1] + edgeOffsets[proc-1]; faceOffsets[proc] = numRecvFaces[proc-1] + faceOffsets[proc-1]; } /* Local to global node number conversion */ for(elem = 0; elem < mesh->numFaces*numCorners; elem++) { ierr = PartitionLocalToGlobalNodeIndex(p, tri->faces[elem], &locFaces[elem]); CHKERRQ(ierr); } for(edge = 0; edge < mesh->numEdges*2; edge++) { ierr = PartitionLocalToGlobalNodeIndex(p, tri->edges[edge], &locEdges[edge]); CHKERRQ(ierr); } /* Collect global arrays */ size = mesh->numNodes*2; ierr = MPI_Gatherv(tri->nodes, size, MPI_DOUBLE, *nodes, numRecvNodes, nodeOffsets, MPI_DOUBLE, 0, p->comm); CHKERRQ(ierr); size = mesh->numNodes; ierr = MPI_Gatherv(tri->markers, size, MPI_INT, *markers, numRecvMarkers, markerOffsets, MPI_INT, 0, p->comm); CHKERRQ(ierr); size = mesh->numEdges*2; ierr = MPI_Gatherv(locEdges, size, MPI_INT, *edges, numRecvEdges, edgeOffsets, MPI_INT, 0, p->comm); CHKERRQ(ierr); size = mesh->numFaces*numCorners; ierr = MPI_Gatherv(locFaces, size, MPI_INT, *faces, numRecvFaces, faceOffsets, MPI_INT, 0, p->comm); CHKERRQ(ierr); /* Cleanup */ if (locFaces != PETSC_NULL) { ierr = PetscFree(locFaces); CHKERRQ(ierr); } if (locEdges != PETSC_NULL) { ierr = PetscFree(locEdges); CHKERRQ(ierr); } ierr = PetscFree(numRecvNodes); CHKERRQ(ierr); ierr = PetscFree(numRecvMarkers); CHKERRQ(ierr); ierr = PetscFree(numRecvEdges); CHKERRQ(ierr); ierr = PetscFree(numRecvFaces); CHKERRQ(ierr); ierr = PetscFree(nodeOffsets); CHKERRQ(ierr); ierr = PetscFree(markerOffsets); CHKERRQ(ierr); ierr = PetscFree(edgeOffsets); CHKERRQ(ierr); ierr = PetscFree(faceOffsets); CHKERRQ(ierr); } else { /* Uniprocessor case */ ierr = PetscMemcpy(*nodes, tri->nodes, numNodes*2 * sizeof(double)); CHKERRQ(ierr); ierr = PetscMemcpy(*markers, tri->markers, numNodes * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(*faces, tri->faces, numFaces*numCorners * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(*edges, tri->edges, numEdges*2 * sizeof(int)); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshCoarsen_Triangular_2D" /* MeshCoarsen_Triangular_2D - Coarsens a two dimensional unstructured mesh using area constraints Collective on Mesh Input Parameters: + mesh - The mesh begin coarsened - area - A function which gives an area constraint when evaluated inside an element Output Parameters: . newmesh - The coarse mesh Note: If PETSC_NULL is given as the 'area' argument, the mesh consisting only of vertices is returned. Level: developer .keywords unstructured grid .seealso GridCoarsen(), MeshRefine(), MeshReform() */ static int MeshCoarsen_Triangular_2D(Mesh mesh, PointFunction area, Mesh *newmesh) { ParameterDict dict; Mesh m; Mesh_Triangular *tri; Partition p; Partition_Triangular_2D *q; Mesh_Triangular *triOld = (Mesh_Triangular *) mesh->data; Partition pOld = mesh->part; Partition_Triangular_2D *qOld = (Partition_Triangular_2D *) pOld->data; int (*refine)(Mesh, PointFunction, Mesh); int *newNodes; /* newNodes[fine node #] = coarse node # or -1 */ int *FineMapping; /* FineMapping[n] = node n in the fine mesh numbering */ int *CoarseMapping; /* CoarseMapping[n] = node n in the coarse mesh numbering */ PetscTruth hasIndex; int *nodePerm, *temp; int numGhostElements, numGhostNodes; int proc, bd, elem, edge, corner, node, newNode, gNode, ghostNode, bdNode, count; int ierr; PetscFunctionBegin; if (area == PETSC_NULL) { if (mesh->numCorners == 3) { /* Return the mesh and increment the reference count instead of copying */ ierr = PetscObjectReference((PetscObject) mesh); CHKERRQ(ierr); *newmesh = mesh; PetscFunctionReturn(0); } else if (mesh->numCorners != 6) { SETERRQ1(PETSC_ERR_SUP, "Number of local nodes %d not supported", mesh->numCorners); } /* We would like to preserve the partition of vertices, although it may be somewhat unbalanced, since this makes mapping from the coarse mesh back to the original much easier. */ ierr = MeshCreate(mesh->comm, &m); CHKERRQ(ierr); ierr = ParameterDictCreate(mesh->comm, &dict); CHKERRQ(ierr); ierr = ParameterDictSetObject(dict, "bdCtx", mesh->bdCtx); CHKERRQ(ierr); ierr = ParameterDictSetInteger(dict, "dim", 2); CHKERRQ(ierr); ierr = ParameterDictSetInteger(dict, "numLocNodes", mesh->numCorners); CHKERRQ(ierr); ierr = PetscObjectSetParameterDict((PetscObject) m, dict); CHKERRQ(ierr); ierr = ParameterDictDestroy(dict); CHKERRQ(ierr); ierr = PetscNew(Mesh_Triangular, &tri); CHKERRQ(ierr); PetscLogObjectMemory(m, sizeof(Mesh_Triangular)); ierr = PetscMemcpy(m->ops, mesh->ops, sizeof(struct _MeshOps)); CHKERRQ(ierr); ierr = PetscStrallocpy(mesh->type_name, &m->type_name); CHKERRQ(ierr); ierr = PetscStrallocpy(mesh->serialize_name, &m->serialize_name); CHKERRQ(ierr); m->data = (void *) tri; m->dim = 2; m->startX = mesh->startX; m->endX = mesh->endX; m->sizeX = mesh->sizeX; m->startY = mesh->startY; m->endY = mesh->endY; m->sizeY = mesh->sizeY; m->isPeriodic = mesh->isPeriodic; m->isPeriodicDim[0] = mesh->isPeriodicDim[0]; m->isPeriodicDim[1] = mesh->isPeriodicDim[1]; m->isPeriodicDim[2] = mesh->isPeriodicDim[2]; m->nodeOrdering = PETSC_NULL; m->partitioned = 1; m->highlightElement = mesh->highlightElement; m->usr = mesh->usr; /* Copy function list */ ierr = PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", (void (**)(void)) &refine); CHKERRQ(ierr); if (refine != PETSC_NULL) { #ifdef PETSC_HAVE_TRIANGLE ierr = PetscObjectComposeFunction((PetscObject) m, "MeshTriangular2D_Refine_C", "MeshRefine_Triangle", (void (*)(void)) MeshRefine_Triangle); #else /* The query needs to return the name also */ ierr = PetscObjectComposeFunction((PetscObject) m, "MeshTriangular2D_Refine_C", "", (void (*)(void)) refine); #endif CHKERRQ(ierr); } /* Create the partition object */ ierr = PartitionCreate(m, &p); CHKERRQ(ierr); ierr = PetscNew(Partition_Triangular_2D, &q); CHKERRQ(ierr); PetscLogObjectMemory(p, sizeof(Partition_Triangular_2D)); ierr = PetscMemcpy(p->ops, pOld->ops, sizeof(struct _PartitionOps)); CHKERRQ(ierr); ierr = PetscObjectChangeTypeName((PetscObject) p, pOld->type_name); CHKERRQ(ierr); ierr = PetscObjectChangeSerializeName((PetscObject) p, pOld->serialize_name); CHKERRQ(ierr); p->data = (void *) q; m->part = p; PetscLogObjectParent(m, p); /* Construct element partition */ p->numProcs = pOld->numProcs; p->rank = pOld->rank; p->ordering = PETSC_NULL; p->numLocElements = pOld->numLocElements; p->numElements = pOld->numElements; p->numOverlapElements = pOld->numOverlapElements; ierr = PetscMalloc((p->numProcs+1) * sizeof(int), &p->firstElement); CHKERRQ(ierr); ierr = PetscMemcpy(p->firstElement, pOld->firstElement, (p->numProcs+1) * sizeof(int)); CHKERRQ(ierr); p->ghostElements = PETSC_NULL; p->ghostElementProcs = PETSC_NULL; numGhostElements = pOld->numOverlapElements - pOld->numLocElements; if (numGhostElements > 0) { ierr = PetscMalloc(numGhostElements * sizeof(int), &p->ghostElements); CHKERRQ(ierr); ierr = PetscMalloc(numGhostElements * sizeof(int), &p->ghostElementProcs); CHKERRQ(ierr); PetscLogObjectMemory(p, numGhostElements*2 * sizeof(int)); ierr = PetscMemcpy(p->ghostElements, pOld->ghostElements, numGhostElements * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(p->ghostElementProcs, pOld->ghostElementProcs, numGhostElements * sizeof(int)); CHKERRQ(ierr); } q->nodeOrdering = PETSC_NULL; q->edgeOrdering = PETSC_NULL; q->numLocEdges = qOld->numLocEdges; q->numEdges = qOld->numEdges; ierr = PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstEdge); CHKERRQ(ierr); ierr = PetscMemcpy(q->firstEdge, qOld->firstEdge, (p->numProcs+1) * sizeof(int)); CHKERRQ(ierr); PetscLogObjectMemory(p, (p->numProcs+1)*2 * sizeof(int)); /* We must handle ghost members since we manually construct the partition */ m->numBd = mesh->numBd; m->numFaces = mesh->numFaces; m->numCorners = 3; m->numEdges = mesh->numEdges; ierr = PetscMalloc(p->numOverlapElements*m->numCorners * sizeof(int), &tri->faces); CHKERRQ(ierr); ierr = PetscMalloc(p->numOverlapElements*3 * sizeof(int), &tri->neighbors); CHKERRQ(ierr); ierr = PetscMalloc(m->numEdges*2 * sizeof(int), &tri->edges); CHKERRQ(ierr); ierr = PetscMalloc(m->numEdges * sizeof(int), &tri->edgemarkers); CHKERRQ(ierr); ierr = PetscMalloc(qOld->numOverlapNodes * sizeof(int), &newNodes); CHKERRQ(ierr); PetscLogObjectMemory(m, (p->numOverlapElements*(m->numCorners + 3) + m->numEdges*3) * sizeof(int)); for(node = 0; node < qOld->numOverlapNodes; node++) newNodes[node] = -1; /* Renumber the vertices and construct the face list */ for(elem = 0, newNode = 0, numGhostNodes = 0; elem < p->numOverlapElements; elem++) { for(corner = 0; corner < 3; corner++) { node = triOld->faces[elem*mesh->numCorners+corner]; if (newNodes[node] == -1) { if (node >= mesh->numNodes) { /* Mark them as ghost nodes */ numGhostNodes++; newNodes[node] = -2; } else { newNodes[node] = newNode++; } } tri->faces[elem*m->numCorners+corner] = node; } } ierr = PetscMemcpy(tri->neighbors, triOld->neighbors, p->numOverlapElements*3 * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(tri->edges, triOld->edges, m->numEdges*2 * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(tri->edgemarkers, triOld->edgemarkers, m->numEdges * sizeof(int)); CHKERRQ(ierr); /* We may have extra ghost nodes from the edges */ for(edge = 0; edge < q->numLocEdges; edge++) { for(corner = 0; corner < 2; corner++) { node = tri->edges[edge*2+corner]; if (newNodes[node] == -1) { if (node >= mesh->numNodes) { /* Mark them as ghost nodes */ numGhostNodes++; newNodes[node] = -2; } } } } /* Compute vertex sizes */ m->numNodes = newNode; m->numVertices = m->numNodes; ierr = MPI_Allreduce(&m->numNodes, &q->numNodes, 1, MPI_INT, MPI_SUM, m->comm); CHKERRQ(ierr); q->numLocNodes = newNode; q->numOverlapNodes = newNode + numGhostNodes; /* Compute vertex partition */ ierr = PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstNode); CHKERRQ(ierr); PetscLogObjectMemory(p, (p->numProcs+1) * sizeof(int)); ierr = MPI_Allgather(&q->numLocNodes, 1, MPI_INT, &q->firstNode[1], 1, MPI_INT, p->comm); CHKERRQ(ierr); for(proc = 1, q->firstNode[0] = 0; proc <= p->numProcs; proc++) { q->firstNode[proc] += q->firstNode[proc-1]; } q->ghostNodes = PETSC_NULL; q->ghostNodeProcs = PETSC_NULL; if (numGhostNodes > 0) { ierr = PetscMalloc(numGhostNodes * sizeof(int), &q->ghostNodes); CHKERRQ(ierr); ierr = PetscMalloc(numGhostNodes * sizeof(int), &q->ghostNodeProcs); CHKERRQ(ierr); PetscLogObjectMemory(p, numGhostNodes*2 * sizeof(int)); } /* Create the mapping from coarse nodes to vertices in the original mesh */ ierr = PetscMalloc(m->numNodes * sizeof(int), &FineMapping); CHKERRQ(ierr); ierr = PetscMalloc(m->numNodes * sizeof(int), &CoarseMapping); CHKERRQ(ierr); for(node = 0, count = 0; node < mesh->numNodes; node++) { newNode = newNodes[node]; if (newNode >= 0) { /* These are all interior nodes */ ierr = PartitionLocalToGlobalNodeIndex(pOld, node, &FineMapping[count]); CHKERRQ(ierr); ierr = PartitionLocalToGlobalNodeIndex(p, newNode, &CoarseMapping[count]); CHKERRQ(ierr); count++; } } if (count != m->numNodes) SETERRQ2(PETSC_ERR_PLIB, "Invalid coarse map size %d should be %d", count, m->numNodes); ierr = AOCreateMapping(m->comm, m->numNodes, FineMapping, CoarseMapping, &m->coarseMap); CHKERRQ(ierr); ierr = PetscFree(FineMapping); CHKERRQ(ierr); ierr = PetscFree(CoarseMapping); CHKERRQ(ierr); /* Setup ghost nodes */ for(ghostNode = mesh->numNodes, count = 0; ghostNode < qOld->numOverlapNodes; ghostNode++) { if (newNodes[ghostNode] == -2) { ierr = PartitionLocalToGlobalNodeIndex(pOld, ghostNode, &gNode); CHKERRQ(ierr); ierr = AOApplicationToPetsc(m->coarseMap, 1, &gNode); CHKERRQ(ierr); q->ghostNodes[count] = gNode; q->ghostNodeProcs[count] = qOld->ghostNodeProcs[ghostNode-mesh->numNodes]; count++; } } if (count != numGhostNodes) SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghost nodes %d should be %d", count, numGhostNodes); /* Cleanup */ ierr = PetscFree(newNodes); CHKERRQ(ierr); /* Resort ghost nodes */ if (numGhostNodes > 0) { ierr = PetscMalloc(numGhostNodes*2 * sizeof(int), &nodePerm); CHKERRQ(ierr); temp = nodePerm + numGhostNodes; for(node = 0; node < numGhostNodes; node++) nodePerm[node] = node; ierr = PetscSortIntWithPermutation(numGhostNodes, q->ghostNodes, nodePerm); CHKERRQ(ierr); for(node = 0; node < numGhostNodes; node++) temp[node] = q->ghostNodes[nodePerm[node]]; for(node = 0; node < numGhostNodes; node++) q->ghostNodes[node] = temp[node]; for(node = 0; node < numGhostNodes; node++) temp[node] = q->ghostNodeProcs[nodePerm[node]]; for(node = 0; node < numGhostNodes; node++) q->ghostNodeProcs[node] = temp[node]; ierr = PetscFree(nodePerm); CHKERRQ(ierr); } /* Copy vertex information */ ierr = PetscMalloc(q->numOverlapNodes*2 * sizeof(double), &tri->nodes); CHKERRQ(ierr); ierr = PetscMalloc(q->numOverlapNodes * sizeof(int), &tri->markers); CHKERRQ(ierr); PetscLogObjectMemory(m, q->numOverlapNodes*2 * sizeof(double)); PetscLogObjectMemory(m, q->numOverlapNodes * sizeof(int)); for(newNode = 0; newNode < q->numOverlapNodes; newNode++) { ierr = PartitionLocalToGlobalNodeIndex(p, newNode, &gNode); CHKERRQ(ierr); ierr = AOPetscToApplication(m->coarseMap, 1, &gNode); CHKERRQ(ierr); ierr = PartitionGlobalToLocalNodeIndex(pOld, gNode, &node); CHKERRQ(ierr); tri->nodes[newNode*2] = triOld->nodes[node*2]; tri->nodes[newNode*2+1] = triOld->nodes[node*2+1]; tri->markers[newNode] = triOld->markers[node]; } /* Renumber nodes */ for(elem = 0; elem < p->numOverlapElements*m->numCorners; elem++) { ierr = PartitionLocalToGlobalNodeIndex(pOld, tri->faces[elem], &tri->faces[elem]); CHKERRQ(ierr); } for(edge = 0; edge < m->numEdges*2; edge++) { ierr = PartitionLocalToGlobalNodeIndex(pOld, tri->edges[edge], &tri->edges[edge]); CHKERRQ(ierr); } ierr = AOApplicationToPetsc(m->coarseMap, p->numOverlapElements*m->numCorners, tri->faces); CHKERRQ(ierr); ierr = AOApplicationToPetsc(m->coarseMap, m->numEdges*2, tri->edges); CHKERRQ(ierr); for(elem = 0; elem < p->numOverlapElements*m->numCorners; elem++) { ierr = PartitionGlobalToLocalNodeIndex(p, tri->faces[elem], &tri->faces[elem]); CHKERRQ(ierr); } for(edge = 0; edge < m->numEdges*2; edge++) { ierr = PartitionGlobalToLocalNodeIndex(p, tri->edges[edge], &tri->edges[edge]); CHKERRQ(ierr); } /* Construct derived and boundary information */ tri->areas = PETSC_NULL; tri->aspectRatios = PETSC_NULL; m->numBdEdges = mesh->numBdEdges; ierr = PetscMalloc(m->numBd * sizeof(int), &tri->bdMarkers); CHKERRQ(ierr); ierr = PetscMalloc((m->numBd+1) * sizeof(int), &tri->bdBegin); CHKERRQ(ierr); ierr = PetscMalloc((m->numBd+1) * sizeof(int), &tri->bdEdgeBegin); CHKERRQ(ierr); ierr = PetscMalloc(m->numBdEdges * sizeof(int), &tri->bdEdges); CHKERRQ(ierr); PetscLogObjectMemory(m, (m->numBd + (m->numBd+1)*2 + m->numBdEdges) * sizeof(int)); ierr = PetscMemcpy(tri->bdMarkers, triOld->bdMarkers, m->numBd * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(tri->bdEdgeBegin, triOld->bdEdgeBegin, (m->numBd+1) * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(tri->bdEdges, triOld->bdEdges, m->numBdEdges * sizeof(int)); CHKERRQ(ierr); tri->bdBegin[0] = 0; for(bd = 0; bd < m->numBd; bd++) { tri->bdBegin[bd+1] = tri->bdBegin[bd]; for(bdNode = triOld->bdBegin[bd]; bdNode < triOld->bdBegin[bd+1]; bdNode++) { node = triOld->bdNodes[bdNode]; ierr = AOMappingHasApplicationIndex(m->coarseMap, node, &hasIndex); CHKERRQ(ierr); if (hasIndex == PETSC_TRUE) tri->bdBegin[bd+1]++; } } m->numBdNodes = tri->bdBegin[m->numBd]; /* Assume closed boundaries */ if (m->numBdNodes != m->numBdEdges) SETERRQ(PETSC_ERR_PLIB, "Invalid mesh boundary"); #ifdef PETSC_USE_BOPT_g ierr = MPI_Scan(&m->numBdNodes, &node, 1, MPI_INT, MPI_SUM, m->comm); CHKERRQ(ierr); if (node != m->numBdNodes*(p->rank+1)) SETERRQ(PETSC_ERR_PLIB, "Invalid mesh boundary"); #endif ierr = PetscMalloc(m->numBdNodes * sizeof(int), &tri->bdNodes); CHKERRQ(ierr); PetscLogObjectMemory(m, m->numBdNodes * sizeof(int)); for(bd = 0, count = 0; bd < m->numBd; bd++) { for(bdNode = triOld->bdBegin[bd]; bdNode < triOld->bdBegin[bd+1]; bdNode++) { node = triOld->bdNodes[bdNode]; ierr = AOMappingHasApplicationIndex(m->coarseMap, node, &hasIndex); CHKERRQ(ierr); if (hasIndex == PETSC_TRUE) { ierr = AOApplicationToPetsc(m->coarseMap, 1, &node); CHKERRQ(ierr); tri->bdNodes[count++] = node; } } } if (count != m->numBdNodes) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node partition"); /* Partition boundary nodes */ ierr = PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstBdNode); CHKERRQ(ierr); PetscLogObjectMemory(p, (p->numProcs+1) * sizeof(int)); q->numLocBdNodes = 0; for(bdNode = 0; bdNode < m->numBdNodes; bdNode++) { newNode = tri->bdNodes[bdNode]; if ((newNode >= q->firstNode[p->rank]) && (newNode < q->firstNode[p->rank+1])) q->numLocBdNodes++; } ierr = MPI_Allgather(&q->numLocBdNodes, 1, MPI_INT, &q->firstBdNode[1], 1, MPI_INT, p->comm); CHKERRQ(ierr); for(proc = 1, q->firstBdNode[0] = 0; proc <= p->numProcs; proc++) { q->firstBdNode[proc] += q->firstBdNode[proc-1]; } q->numBdNodes = q->firstBdNode[p->numProcs]; /* Process ghost boundary nodes */ q->numOverlapBdNodes = q->numLocBdNodes; for(node = 0; node < numGhostNodes; node++) { if (tri->markers[m->numNodes+node] != 0) q->numOverlapBdNodes++; } q->ghostBdNodes = PETSC_NULL; if (q->numOverlapBdNodes > q->numLocBdNodes) { ierr = PetscMalloc((q->numOverlapBdNodes - q->numLocBdNodes) * sizeof(int), &q->ghostBdNodes); CHKERRQ(ierr); for(node = 0, bdNode = 0; node < numGhostNodes; node++) { if (tri->markers[m->numNodes+node] != 0) q->ghostBdNodes[bdNode++] = node; } if (bdNode != q->numOverlapBdNodes - q->numLocBdNodes) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node partition"); } /* Copy holes from previous mesh */ m->numHoles = mesh->numHoles; m->holes = PETSC_NULL; if (m->numHoles > 0) { ierr = PetscMalloc(m->numHoles*2 * sizeof(double), &m->holes); CHKERRQ(ierr); ierr = PetscMemcpy(m->holes, mesh->holes, m->numHoles*2 * sizeof(double)); CHKERRQ(ierr); PetscLogObjectMemory(m, m->numHoles*2 * sizeof(double)); } /* Calculate maximum degree of vertices */ m->maxDegree = mesh->maxDegree; ierr = PetscMalloc(m->numNodes * sizeof(int), &tri->degrees); CHKERRQ(ierr); ierr = PetscMalloc(m->maxDegree * sizeof(int), &m->support); CHKERRQ(ierr); for(newNode = 0; newNode < m->numNodes; newNode++) { ierr = PartitionLocalToGlobalNodeIndex(p, newNode, &gNode); CHKERRQ(ierr); ierr = AOPetscToApplication(m->coarseMap, 1, &gNode); CHKERRQ(ierr); ierr = PartitionGlobalToLocalNodeIndex(pOld, gNode, &node); CHKERRQ(ierr); tri->degrees[newNode] = triOld->degrees[node]; } m->supportTaken = PETSC_FALSE; #ifdef PETSC_USE_BOPT_g /* Check mesh integrity */ ierr = MeshDebug_Triangular_2D(m, PETSC_TRUE); CHKERRQ(ierr); #endif /* Initialize save space */ tri->nodesOld = PETSC_NULL; } else { SETERRQ(PETSC_ERR_SUP, "No coarsening strategies currently supported"); } ierr = PetscObjectSetName((PetscObject) m, "Coarse Mesh"); CHKERRQ(ierr); *newmesh = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRefine_Triangular_2D" /* MeshRefine_Triangular_2D - Refines a two dimensional unstructured mesh using area constraints Collective on Mesh Input Parameters: + mesh - The mesh begin refined - area - A function which gives an area constraint when evaluated inside an element Output Parameter: . newmesh - The refined mesh Level: developer .keywords unstructured grid .seealso GridRefine(), MeshCoarsen(), MeshReform() */ static int MeshRefine_Triangular_2D(Mesh mesh, PointFunction area, Mesh *newmesh) { ParameterDict dict; Mesh m; int (*f)(Mesh, PointFunction, Mesh); int ierr; PetscFunctionBegin; ierr = MeshCreate(mesh->comm, &m); CHKERRQ(ierr); ierr = ParameterDictCreate(mesh->comm, &dict); CHKERRQ(ierr); ierr = ParameterDictSetObject(dict, "bdCtx", mesh->bdCtx); CHKERRQ(ierr); ierr = ParameterDictSetInteger(dict, "dim", 2); CHKERRQ(ierr); ierr = ParameterDictSetInteger(dict, "numLocNodes", mesh->numCorners); CHKERRQ(ierr); ierr = PetscObjectSetParameterDict((PetscObject) m, dict); CHKERRQ(ierr); ierr = ParameterDictDestroy(dict); CHKERRQ(ierr); ierr = PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", (void (**)(void)) &f); CHKERRQ(ierr); if (f == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, "Mesh has no refinement function"); ierr = (*f)(mesh, area, m); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) m, "Refined Mesh"); CHKERRQ(ierr); *newmesh = m; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshResetNodes_Triangular_2D" static int MeshResetNodes_Triangular_2D(Mesh mesh, PetscTruth resetBd) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int *elements = tri->faces; int *markers = tri->markers; double *nodes = tri->nodes; int numCorners = mesh->numCorners; int elem, corner, node1, node2, midNode, midCorner; PetscFunctionBegin; if (numCorners == 3) PetscFunctionReturn(0); if (mesh->isPeriodic == PETSC_TRUE) { for(elem = 0; elem < mesh->part->numOverlapElements; elem++) { for(corner = 0; corner < 3; corner++) { node1 = elements[elem*numCorners+corner]; node2 = elements[elem*numCorners+((corner+1)%3)]; if ((markers[node1] != markers[node2]) || (markers[node1] == 0) || (resetBd == PETSC_TRUE)) { midCorner = (corner+2)%3 + 3; midNode = elements[elem*numCorners+midCorner]; nodes[midNode*2] = MeshPeriodicX(mesh, 0.5*(nodes[node1*2] + MeshPeriodicRelativeX(mesh, nodes[node2*2], nodes[node1*2]))); nodes[midNode*2+1] = MeshPeriodicY(mesh, 0.5*(nodes[node1*2+1] + MeshPeriodicRelativeY(mesh, nodes[node2*2+1], nodes[node1*2+1]))); } } } } else { for(elem = 0; elem < mesh->part->numOverlapElements; elem++) { for(corner = 0; corner < 3; corner++) { node1 = elements[elem*numCorners+corner]; node2 = elements[elem*numCorners+((corner+1)%3)]; if ((markers[node1] != markers[node2]) || (markers[node1] == 0) || (resetBd == PETSC_TRUE)) { midCorner = (corner+2)%3 + 3; midNode = elements[elem*numCorners+midCorner]; nodes[midNode*2] = 0.5*(nodes[node1*2] + nodes[node2*2]); nodes[midNode*2+1] = 0.5*(nodes[node1*2+1] + nodes[node2*2+1]); } } } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshPartition_Triangular_2D" int MeshPartition_Triangular_2D(Mesh mesh) { int ierr; PetscFunctionBegin; ierr = PartitionCreateTriangular2D(mesh, &mesh->part); mesh->partitioned = 1; PetscFunctionReturn(ierr); } #undef __FUNCT__ #define __FUNCT__ "MeshUpdateBoundingBox_Triangular_2D" int MeshUpdateBoundingBox_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; double *nodes = tri->nodes; Partition part; int dim, numOverlapNodes; int node; int ierr; /* Calculate the local bounding box */ PetscFunctionBegin; ierr = MeshGetPartition(mesh, &part); CHKERRQ(ierr); ierr = MeshGetDimension(mesh, &dim); CHKERRQ(ierr); ierr = PartitionGetNumOverlapNodes(part, &numOverlapNodes); CHKERRQ(ierr); if (numOverlapNodes > 0) { mesh->locStartX = mesh->locEndX = nodes[0]; mesh->locStartY = mesh->locEndY = nodes[1]; } else { mesh->locStartX = mesh->locEndX = 0.0; mesh->locStartY = mesh->locEndY = 0.0; } for(node = 0; node < numOverlapNodes; node++) { if (nodes[node*dim] < mesh->locStartX) mesh->locStartX = nodes[node*dim]; if (nodes[node*dim] > mesh->locEndX) mesh->locEndX = nodes[node*dim]; if (nodes[node*dim+1] < mesh->locStartY) mesh->locStartY = nodes[node*dim+1]; if (nodes[node*dim+1] > mesh->locEndY) mesh->locEndY = nodes[node*dim+1]; } mesh->locSizeX = mesh->locEndX - mesh->locStartX; mesh->locSizeY = mesh->locEndY - mesh->locStartY; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetElemNeighbor_Triangular_2D" int MeshGetElemNeighbor_Triangular_2D(Mesh mesh, int elem, int corner, int *neighbor) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; #ifdef PETSC_USE_BOPT_g int numOverlapElements; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapElements(mesh->part, &numOverlapElements); CHKERRQ(ierr); if ((elem < 0) || (elem >= numOverlapElements)) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid element %d should be in [0,%d)", elem, numOverlapElements); } if ((corner < 0) || (corner >= mesh->numCorners)) { SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid corner %d should be in [0,%d)", corner, mesh->numCorners); } #endif *neighbor = tri->neighbors[elem*3+corner]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeCoords_Triangular_2D" int MeshGetNodeCoords_Triangular_2D(Mesh mesh, int node, double *x, double *y, double *z) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int dim = mesh->dim; #ifdef PETSC_USE_BOPT_g int numOverlapNodes; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified"); #endif if (x != PETSC_NULL) *x = tri->nodes[node*dim]; if (y != PETSC_NULL) *y = tri->nodes[node*dim+1]; if (z != PETSC_NULL) *z = 0.0; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSetNodeCoords_Triangular_2D" int MeshSetNodeCoords_Triangular_2D(Mesh mesh, int node, double x, double y, double z) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int dim = mesh->dim; #ifdef PETSC_USE_BOPT_g int numOverlapNodes; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified"); #endif tri->nodes[node*dim] = MeshPeriodicX(mesh, x); tri->nodes[node*dim+1] = MeshPeriodicY(mesh, y); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeCoordsSaved_Triangular_2D" int MeshGetNodeCoordsSaved_Triangular_2D(Mesh mesh, int node, double *x, double *y, double *z) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int dim = mesh->dim; #ifdef PETSC_USE_BOPT_g int numOverlapNodes; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified"); #endif if (tri->nodesOld != PETSC_NULL) { if (x != PETSC_NULL) *x = tri->nodesOld[node*dim]; if (y != PETSC_NULL) *y = tri->nodesOld[node*dim+1]; if (z != PETSC_NULL) *z = 0.0; } else { if (x != PETSC_NULL) *x = tri->nodes[node*dim]; if (y != PETSC_NULL) *y = tri->nodes[node*dim+1]; if (z != PETSC_NULL) *z = 0.0; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeBoundary_Triangular_2D" int MeshGetNodeBoundary_Triangular_2D(Mesh mesh, int node, int *bd) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; #ifdef PETSC_USE_BOPT_g int numOverlapNodes; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if ((node < 0) || (node >= numOverlapNodes)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", node, numOverlapNodes); } #endif *bd = tri->markers[node]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeFromElement_Triangular_2D" int MeshGetNodeFromElement_Triangular_2D(Mesh mesh, int elem, int corner, int *node) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int numCorners = mesh->numCorners; PetscFunctionBegin; *node = tri->faces[elem*numCorners+corner]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeFromEdge_Triangular_2D" int MeshGetNodeFromEdge_Triangular_2D(Mesh mesh, int edge, int corner, int *node) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; PetscFunctionBegin; *node = tri->edges[edge*2+corner]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshNodeIsVertex_Triangular_2D" int MeshNodeIsVertex_Triangular_2D(Mesh mesh, int node, PetscTruth *isVertex) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; #ifdef PETSC_USE_BOPT_g int numOverlapNodes; int ierr; #endif PetscFunctionBegin; #ifdef PETSC_USE_BOPT_g ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if ((node < 0) || (node >= numOverlapNodes)) { SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", node, numOverlapNodes); } #endif if (tri->degrees[node] == 0) *isVertex = PETSC_FALSE; else *isVertex = PETSC_TRUE; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetBoundarySize_Triangular_2D" int MeshGetBoundarySize_Triangular_2D(Mesh mesh, int boundary, int *size) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int b; /* Canonical boundary number */ int ierr; PetscFunctionBegin; /* Find canonical boundary number */ ierr = MeshGetBoundaryIndex(mesh, boundary, &b); CHKERRQ(ierr); *size = tri->bdBegin[b+1] - tri->bdBegin[b]; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetBoundaryIndex_Triangular_2D" int MeshGetBoundaryIndex_Triangular_2D(Mesh mesh, int boundary, int *index) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int b; /* Canonical boundary number */ PetscFunctionBegin; /* Find canonical boundary number */ for(b = 0; b < mesh->numBd; b++) if (tri->bdMarkers[b] == boundary) break; if (b == mesh->numBd) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid boundary marker %d", boundary); *index = b; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetBoundaryStart_Triangular_2D" int MeshGetBoundaryStart_Triangular_2D(Mesh mesh, int boundary, PetscTruth ghost, int *node) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int b; /* Canonical boundary number */ int ierr; PetscFunctionBegin; /* Find canonical boundary number */ ierr = MeshGetBoundaryIndex(mesh, boundary, &b); CHKERRQ(ierr); if (mesh->activeBd != -1) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Already iterating over boundary %d", mesh->activeBd); /* Find first boundary node */ mesh->activeBd = b; mesh->activeBdOld = b; mesh->activeBdNode = tri->bdBegin[b] - 1; ierr = MeshGetBoundaryNext(mesh, boundary, ghost, node); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetBoundaryNext_Triangular_2D" int MeshGetBoundaryNext_Triangular_2D(Mesh mesh, int boundary, PetscTruth ghost, int *node) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; Partition p = mesh->part; Partition_Triangular_2D *q = (Partition_Triangular_2D *) p->data; int numNodes = mesh->numNodes; int b; /* Canonical boundary number */ int tempNode; /* Canonical local node number */ int ierr; PetscFunctionBegin; /* Find canonical boundary number */ ierr = MeshGetBoundaryIndex(mesh, boundary, &b); CHKERRQ(ierr); /* Find next boundary node */ if (mesh->activeBd != b) SETERRQ1(PETSC_ERR_ARG_WRONG, "Boundary %d not active", boundary); do { mesh->activeBdNode++; if (mesh->activeBdNode == tri->bdBegin[b+1]) { mesh->activeBd = mesh->activeBdNode = -1; *node = -1; PetscFunctionReturn(0); } tempNode = tri->bdNodes[mesh->activeBdNode] - q->firstNode[p->rank]; /* Translate ghost nodes */ if (ghost == PETSC_TRUE) { if ((tempNode < 0) || (tempNode >= numNodes)) { ierr = PartitionGhostNodeIndex_Private(p, tri->bdNodes[mesh->activeBdNode], &tempNode); if (ierr == 0) { tempNode += numNodes; break; } } } } while((tempNode < 0) || (tempNode >= numNodes)); *node = tempNode; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetActiveBoundary_Triangular_2D" int MeshGetActiveBoundary_Triangular_2D(Mesh mesh, int *boundary) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; PetscFunctionBegin; if (mesh->activeBdOld >= 0) { *boundary = tri->bdMarkers[mesh->activeBdOld]; } else { *boundary = 0; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshGetNodeSupport_Triangular_2D" int MeshGetNodeSupport_Triangular_2D(Mesh mesh, int node, int startElem, int *degree, int **support) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int numOverlapElements = mesh->numFaces; int numCorners = mesh->numCorners; int *elements = tri->faces; int *neighbors = tri->neighbors; int deg = -1; PetscTruth found; int firstElem; int prevElem, nextElem; int elem, neighbor, corner; PetscFunctionBegin; /* First check suggested element */ firstElem = -1; for(corner = 0; corner < numCorners; corner++) { if (elements[startElem*numCorners+corner] == node) { firstElem = startElem; if (corner >= 3) { deg = 1; } else { deg = 0; } break; } } /* Locate an element containing the node */ if (mesh->part != PETSC_NULL) { numOverlapElements = mesh->part->numOverlapElements; } for(elem = 0; (elem < numOverlapElements) && (firstElem < 0); elem++) { for(corner = 0; corner < numCorners; corner++) { if (elements[elem*numCorners+corner] == node) { firstElem = elem; if (corner >= 3) { deg = 1; } else { deg = 0; } break; } } } if (firstElem < 0) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Node %d not found in mesh", node); } /* Check for midnode */ if (deg != 0) { mesh->support[0] = firstElem; mesh->support[1] = -1; /* Find other element containing node */ for(neighbor = 0; neighbor < 3; neighbor++) { elem = neighbors[firstElem*3+neighbor]; if (elem < 0) continue; for(corner = 3; corner < numCorners; corner++) { if (elements[elem*numCorners+corner] == node) { mesh->support[1] = elem; deg++; break; } } if (corner < numCorners) break; } #ifdef PETSC_USE_BOPT_g if (mesh->support[1] == -1) { if (tri->markers[node] == 0) SETERRQ1(PETSC_ERR_PLIB, "Interior node %d with incomplete support", node); } #endif *support = mesh->support; *degree = deg; PetscFunctionReturn(0); } /* Walk around node */ prevElem = -1; nextElem = -1; found = PETSC_TRUE; while((nextElem != firstElem) && (found == PETSC_TRUE)) { /* First time through */ if (prevElem == -1) nextElem = firstElem; /* Add element to the list */ mesh->support[deg++] = nextElem; /* Look for a neighbor containing the node */ found = PETSC_FALSE; for(neighbor = 0; neighbor < 3; neighbor++) { /* Reject the element we just came from (we could reject the opposite face also) */ elem = neighbors[nextElem*3+neighbor]; if ((elem == prevElem) || (elem < 0)) continue; for(corner = 0; corner < 3; corner++) { if (elements[elem*numCorners+corner] == node) { prevElem = nextElem; nextElem = elem; found = PETSC_TRUE; break; } } if (corner < 3) break; } if (found == PETSC_FALSE) { #ifdef PETSC_USE_BOPT_g /* Make sure that this is a boundary node */ if (tri->markers[node] == 0) SETERRQ(PETSC_ERR_PLIB, "Interior node with incomplete support"); #endif /* Also walk the other way for a boundary node -- We could be nice and reorder the elements afterwards to all be adjacent */ if (deg > 1) { prevElem = mesh->support[1]; nextElem = mesh->support[0]; found = PETSC_TRUE; while(found == PETSC_TRUE) { /* Look for a neighbor containing the node */ found = PETSC_FALSE; for(neighbor = 0; neighbor < 3; neighbor++) { /* Reject the element we just came from (we could reject the opposite face also) */ elem = neighbors[nextElem*3+neighbor]; if ((elem == prevElem) || (elem < 0)) continue; for(corner = 0; corner < 3; corner++) { if (elements[elem*numCorners+corner] == node) { prevElem = nextElem; nextElem = elem; found = PETSC_TRUE; /* Add element to the list */ mesh->support[deg++] = nextElem; break; } } if (corner < 3) break; } } } } #ifdef PETSC_USE_BOPT_g /* Check for overflow */ if (deg > mesh->maxDegree) SETERRQ1(PETSC_ERR_MEM, "Maximum degree %d exceeded", mesh->maxDegree); #endif } *support = mesh->support; *degree = deg; PetscFunctionReturn(0); } #if 0 #undef __FUNCT__ #define __FUNCT__ "MeshLocatePoint_Triangular_2D_Directed" int MeshLocatePoint_Triangular_2D_Directed(Mesh mesh, double xx, double yy, double zz, int *containingElem) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; double start_x, start_y; /* Coordinates of a vertex */ double v_x, v_y; /* Coordinates of an edge vector */ double p_x, p_y; /* Coordinates of the vector from a vertex (start) to a point */ double cross; /* The magnitude of the cross product of two vectors */ PetscTruth inside[3]; /* Are we in the correct halfplane? */ PetscTruth found; /* We have found the element containing a point */ int nextElem; /* The next element to search */ int prevElem; /* The last element searched */ int numCorners = mesh->numCorners; int numElements = mesh->numFaces; double *nodes = tri->nodes; int *elements = tri->faces; int *neighbors = tri->neighbors; double x = xx; double y = yy; double coords[6]; int neighbor; int elem, edge; PetscFunctionBegin; *containingElem = -1; /* Quick bounding rectangle check */ if ((mesh->isPeriodic == PETSC_FALSE) && ((mesh->locStartX > x) || (mesh->locEndX < x) || (mesh->locStartY > y) || (mesh->locEndY < y))) PetscFunctionReturn(0); /* Simple search */ for(elem = 0, found = PETSC_FALSE; elem < numElements; elem++) { /* A point lies within a triangle is it is on the inner side of every edge -- We take the cross product of the edge oriented counterclockwise with the vector from the edge start to the point: 2 |\ | \ | \ 0---1----> v_0 \ \ \theta \ P If |v_i \cross P| \equiv |v_i| |P| \sin\theta > 0 then the point lies in the correct half-plane. This has the advantage of only using elementary arithmetic. Robust predicates could be used to catch points near a boundary. */ if (mesh->isPeriodic == PETSC_TRUE) { coords[0] = nodes[elements[elem*numCorners]*2]; coords[1] = nodes[elements[elem*numCorners]*2+1]; coords[2] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+1]*2], nodes[elements[elem*numCorners]*2]); coords[3] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+1]*2+1], nodes[elements[elem*numCorners]*2+1]); coords[4] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+2]*2], nodes[elements[elem*numCorners]*2]); coords[5] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+2]*2+1], nodes[elements[elem*numCorners]*2+1]); x = MeshPeriodicRelativeX(mesh, xx, nodes[elements[elem*numCorners]*2]); y = MeshPeriodicRelativeY(mesh, yy, nodes[elements[elem*numCorners]*2+1]); } for(edge = 0; edge < 3; edge++) { /* First edge 0--1 */ if (mesh->isPeriodic == PETSC_TRUE) { start_x = coords[edge*2]; start_y = coords[edge*2+1]; v_x = coords[((edge+1)%3)*2] - start_x; v_y = coords[((edge+1)%3)*2+1] - start_y; } else { start_x = nodes[elements[elem*numCorners+edge]*2]; start_y = nodes[elements[elem*numCorners+edge]*2+1]; v_x = nodes[elements[elem*numCorners+((edge+1)%3)]*2] - start_x; v_y = nodes[elements[elem*numCorners+((edge+1)%3)]*2+1] - start_y; } p_x = x - start_x; p_y = y - start_y; cross = v_x*p_y - v_y*p_x; /* We will use this routine to find integration points which lie on the boundary so we must give a little slack to the test. I do not think this will bias our results. */ if ((cross >= 0) || (PetscAbsScalar(cross) < PETSC_MACHINE_EPSILON)) { inside[edge] = PETSC_TRUE; } else { inside[edge] = PETSC_FALSE; } } /* This information could be used for a directed search -- \ 4 T| \ F| \ F| \ | \| 2 T |\ T 5 T | \ F 3 F |T \ T ----0---1------ F | F \ F 6 T | T 1 \ F 2 F | T \ T and clearly all F's means that the area of the triangle is 0. The regions are ordered so that all the Hamming distances are 1 (Gray code). */ if (inside[2] == PETSC_TRUE) { if (inside[1] == PETSC_TRUE) { if (inside[0] == PETSC_TRUE) { /* Region 0 -- TTT */ found = PETSC_TRUE; } else { /* Region 1 -- FTT */ nextElem = neighbors[elem*3+2]; } } else { if (inside[0] == PETSC_TRUE) { /* Region 3 -- TFT */ nextElem = neighbors[elem*3]; } else { /* Region 2 -- FFT */ if (neighbors[elem*3] >= 0) { nextElem = neighbors[elem*3]; } else { nextElem = neighbors[elem*3+2]; } } } } else { if (inside[1] == PETSC_TRUE) { if (inside[0] == PETSC_TRUE) { /* Region 5 -- TTF */ nextElem = neighbors[elem*3+1]; } else { /* Region 6 -- FTF */ if (neighbors[elem*3+1] >= 0) { nextElem = neighbors[elem*3+1]; } else { nextElem = neighbors[elem*3+2]; } } } else { if (inside[0] == PETSC_TRUE) { /* Region 4 -- TFF */ if (neighbors[elem*3] >= 0) { nextElem = neighbors[elem*3]; } else { nextElem = neighbors[elem*3+1]; } } else { /* Region 7 -- FFF */ PetscPrintf(PETSC_COMM_SELF, "Element %d: (%g, %g)-(%g, %g)-(%g, %g) and (%g, %g)\n", elem, nodes[elements[elem*numCorners]*2], nodes[elements[elem*numCorners]*2+1], nodes[elements[elem*numCorners+1]*2], nodes[elements[elem*numCorners+1]*2+1], nodes[elements[elem*numCorners+2]*2], nodes[elements[elem*numCorners+2]*2+1], x, y); SETERRQ1(PETSC_ERR_MESH_NULL_ELEM, "Element %d has no area", elem); } } } /* We found the point */ if (found) { *containingElem = elem; break; } /* Choose next direction */ if (nextElem < 0) { /* Choose randomly */ for(neighbor = 0; neighbor < 3; neighbor++) { tempElem = neighbors[elem*3+neighbor]; if ((tempElem >= 0) && (tempElem != prevElem)) { nextElem = tempElem; break; } } if (nextElem < 0) { SETERRQ1(PETSC_ERR_MAX_ITER, "Element search hit dead end in element %d", elem); } } prevElem = elem; } PetscFunctionReturn(0); } #endif #undef __FUNCT__ #define __FUNCT__ "MeshLocatePoint_Triangular_2D" int MeshLocatePoint_Triangular_2D(Mesh mesh, double xx, double yy, double zz, int *containingElem) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; double start_x, start_y; /* Coordinates of a vertex */ double v_x, v_y; /* Coordinates of an edge vector */ double p_x, p_y; /* Coordinates of the vector from a vertex (start) to a point */ double cross; /* The magnitude of the cross product of two vectors */ PetscTruth inside[3]; /* Are we in the correct halfplane? */ PetscTruth found; /* We have found the element containing a point */ int numCorners = mesh->numCorners; int numElements = mesh->numFaces; double *nodes = tri->nodes; int *markers = tri->markers; int *elements = tri->faces; int *neighbors = tri->neighbors; double x = xx; double y = yy; double coords[6]; double pNormSq; int elem, edge; PetscFunctionBegin; *containingElem = -1; /* Quick bounding rectangle check */ if ((mesh->isPeriodic == PETSC_FALSE) && ((mesh->locStartX > x) || (mesh->locEndX < x) || (mesh->locStartY > y) || (mesh->locEndY < y))) PetscFunctionReturn(0); /* Simple search */ for(elem = 0, found = PETSC_FALSE; elem < numElements; elem++) { /* A point lies within a triangle is it is on the inner side of every edge -- We take the cross product of the edge oriented counterclockwise with the vector from the edge start to the point: 2 |\ | \ | \ 0---1----> v_0 \ \ \theta \ P If |v_i \cross P| \equiv |v_i| |P| \sin\theta > 0 then the point lies in the correct half-plane. This has the advantage of only using elementary arithmetic. Robust predicates could be used to catch points near a boundary. */ if (mesh->isPeriodic == PETSC_TRUE) { coords[0] = nodes[elements[elem*numCorners]*2]; coords[1] = nodes[elements[elem*numCorners]*2+1]; coords[2] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+1]*2], nodes[elements[elem*numCorners]*2]); coords[3] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+1]*2+1], nodes[elements[elem*numCorners]*2+1]); coords[4] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+2]*2], nodes[elements[elem*numCorners]*2]); coords[5] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+2]*2+1], nodes[elements[elem*numCorners]*2+1]); x = MeshPeriodicRelativeX(mesh, xx, nodes[elements[elem*numCorners]*2]); y = MeshPeriodicRelativeY(mesh, yy, nodes[elements[elem*numCorners]*2+1]); } for(edge = 0; edge < 3; edge++) { /* First edge 0--1 */ if (mesh->isPeriodic == PETSC_TRUE) { start_x = coords[edge*2]; start_y = coords[edge*2+1]; v_x = coords[((edge+1)%3)*2] - start_x; v_y = coords[((edge+1)%3)*2+1] - start_y; } else { start_x = nodes[elements[elem*numCorners+edge]*2]; start_y = nodes[elements[elem*numCorners+edge]*2+1]; v_x = nodes[elements[elem*numCorners+((edge+1)%3)]*2] - start_x; v_y = nodes[elements[elem*numCorners+((edge+1)%3)]*2+1] - start_y; } p_x = x - start_x; p_y = y - start_y; cross = v_x*p_y - v_y*p_x; pNormSq = p_x*p_x+p_y*p_y; /* Check for identical points, being looser at the boundary */ if (pNormSq < PETSC_MACHINE_EPSILON) { *containingElem = elem; PetscFunctionReturn(0); } else if (markers[elements[elem*numCorners+edge]] != 0) { if (pNormSq < 0.001*PetscMax(v_x*v_x, v_y*v_y)) { *containingElem = elem; PetscFunctionReturn(0); } } if (cross >= 0) { inside[edge] = PETSC_TRUE; } else if ((neighbors[elem*3+((edge+2)%3)] < 0) && (cross*cross < 0.087155743*(v_x*v_x+v_y*v_y)*pNormSq)) { /* We are losing precision here (I think) for large problems (since points are turning up inside particles) and therefore I am putting a 5 degree limit on this check (sin 5 = 0.087155743), which is only active for boundary edges. */ inside[edge] = PETSC_TRUE; } else { inside[edge] = PETSC_FALSE; } } /* This information could be used for a directed search -- \ 4 T| \ F| \ F| \ | \| 2 T |\ T 5 T | \ F 3 F |T \ T ----0---1------ F | F \ F 6 T | T 1 \ F 2 F | T \ T and clearly all F's means that the area of the triangle is 0. The regions are ordered so that all the Hamming distances are 1 (Gray code). */ if ((inside[2] == PETSC_TRUE) && (inside[1] == PETSC_TRUE) && (inside[0] == PETSC_TRUE)) { double minX = coords[0]; double maxX = coords[0]; double minY = coords[1]; double maxY = coords[1]; double offX, offY; if (coords[2] < minX) minX = coords[2]; if (coords[2] > maxX) maxX = coords[2]; if (coords[3] < minY) minY = coords[3]; if (coords[3] > maxY) maxY = coords[3]; if (coords[4] < minX) minX = coords[4]; if (coords[4] > maxX) maxX = coords[4]; if (coords[5] < minY) minY = coords[5]; if (coords[5] > maxY) maxY = coords[5]; offX = 0.01*(maxX - minX); offY = 0.01*(maxY - minY); if ((x >= minX-offX) && (x <= maxX+offX) && (y >= minY-offY) && (y <= maxY+offY)) { /* Region 0 -- TTT */ found = PETSC_TRUE; } } else if ((inside[2] == PETSC_FALSE) && (inside[1] == PETSC_FALSE) && (inside[0] == PETSC_FALSE)) { /* Region 7 -- FFF */ PetscPrintf(PETSC_COMM_SELF, "Element %d: (%g, %g)-(%g, %g)-(%g, %g) and (%g, %g)\n", elem, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5], x, y); SETERRQ1(PETSC_ERR_MESH_NULL_ELEM, "Element %d has no area or is inverted", elem); } /* We found the point */ if (found) { *containingElem = elem; break; } } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshNearestNode_Triangular_2D" int MeshNearestNode_Triangular_2D(Mesh mesh, double x, double y, double z, PetscTruth outsideDomain, int *node) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; Partition p = mesh->part; Partition_Triangular_2D *q = (Partition_Triangular_2D *) p->data; int numCorners = mesh->numCorners; int *elements = tri->faces; int numNodes = mesh->numNodes; double *nodes = tri->nodes; int *firstNode = q->firstNode; int rank = p->rank; int elem; double minDist, dist, allMinDist; int minNode, globalMinNode, allMinNode; int corner, n; int ierr; PetscFunctionBegin; if (outsideDomain == PETSC_FALSE) { ierr = MeshLocatePoint(mesh, x, y, z, &elem); CHKERRQ(ierr); if (elem >= 0) { minDist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners]*2] - x)) + PetscSqr(MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners]*2+1] - y)); minNode = elements[elem*numCorners]; for(corner = 1; corner < numCorners; corner++) { dist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners+corner]*2] - x)) + PetscSqr(MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners+corner]*2+1] - y)); if (dist < minDist) { minDist = dist; minNode = elements[elem*numCorners+corner]; } } ierr = PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode); CHKERRQ(ierr); } else { minNode = -1; globalMinNode = -1; } /* The minimum node might still be a ghost node */ ierr = MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, p->comm); CHKERRQ(ierr); if ((allMinNode >= firstNode[rank]) && (allMinNode < firstNode[rank+1])) { *node = allMinNode - firstNode[rank]; } else { *node = -1; } if (allMinNode < 0) PetscFunctionReturn(1); } else { /* Brute force check */ minNode = 0; minDist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[0*2] - x)) + PetscSqr(MeshPeriodicDiffY(mesh, nodes[0*2+1] - y)); for(n = 1; n < numNodes; n++) { dist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[n*2] - x)) + PetscSqr(MeshPeriodicDiffY(mesh, nodes[n*2+1] - y)); if (dist < minDist) { minDist = dist; minNode = n; } } /* Find the minimum distance */ ierr = MPI_Allreduce(&minDist, &allMinDist, 1, MPI_DOUBLE, MPI_MIN, p->comm); CHKERRQ(ierr); if (minDist == allMinDist) { ierr = PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode); CHKERRQ(ierr); } else { globalMinNode = -1; } /* We might still have ties */ ierr = MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, p->comm); CHKERRQ(ierr); if ((allMinNode >= firstNode[rank]) && (allMinNode < firstNode[rank+1])) { *node = allMinNode - firstNode[rank]; } else { *node = -1; } if (allMinNode < 0) PetscFunctionReturn(1); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshSaveMesh_Triangular_2D" int MeshSaveMesh_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int dim = mesh->dim; int numOverlapNodes; int ierr; PetscFunctionBegin; ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if (tri->nodesOld == PETSC_NULL) { ierr = PetscMalloc(numOverlapNodes*dim * sizeof(double), &tri->nodesOld); CHKERRQ(ierr); PetscLogObjectMemory(mesh, numOverlapNodes*dim * sizeof(double)); } ierr = PetscMemcpy(tri->nodesOld, tri->nodes, numOverlapNodes*dim * sizeof(double)); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshRestoreMesh_Triangular_2D" int MeshRestoreMesh_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; int dim = mesh->dim; int numOverlapNodes; int ierr; PetscFunctionBegin; ierr = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes); CHKERRQ(ierr); if (tri->nodesOld != PETSC_NULL) { ierr = PetscMemcpy(tri->nodes, tri->nodesOld, numOverlapNodes*dim * sizeof(double)); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "MeshIsDistorted_Triangular_2D" int MeshIsDistorted_Triangular_2D(Mesh mesh, PetscTruth *flag) { Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data; double *nodes = tri->nodes; int *faces = tri->faces; int numCorners = mesh->numCorners; int numElements = mesh->numFaces; double maxAspectRatio = mesh->maxAspectRatio; double area; /* Twice the area of the triangle */ double aspectRatio; /* (Longest edge)^2/(Twice the area of the triangle) */ double max; PetscTruth locFlag; double x21, y21, x31, y31, x32, y32; int elem; int locRetVal; int retVal; int ierr; PetscFunctionBegin; for(elem = 0, locFlag = PETSC_FALSE, locRetVal = 0, max = 0.0; elem < numElements; elem++) { /* Find the area */ if (mesh->isPeriodic == PETSC_TRUE) { x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]); y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]); x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]); y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]); x32 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners+1]*2]); y32 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]); } else { x21 = nodes[faces[elem*numCorners+1]*2] - nodes[faces[elem*numCorners]*2]; y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]; x31 = nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners]*2]; y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]; x32 = nodes[faces[elem*numCorners+2]*2] - nodes[faces[elem*numCorners+1]*2]; y32 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]; } area = x21*y31 - x31*y21; /* Check for inverted elements */ if (area < 0.0) { locFlag = PETSC_TRUE; locRetVal = 1; break; } else if (area == 0.0) { SETERRQ1(PETSC_ERR_ARG_CORRUPT, "Element %d has no area", elem); } /* Find the aspect ratio = (longest edge length)^2 / (2 area) */ aspectRatio = PetscMax(x21*x21 + y21*y21, x31*x31 + y31*y31); aspectRatio = PetscMax(x32*x32 + y32*y32, aspectRatio); aspectRatio /= area; /* We cannot bail here since we must check for inverted elements */ max = PetscMax(max, aspectRatio); } if (max > maxAspectRatio) { PetscLogInfo(mesh, "Aspect ratio too large in element %d: %g > %g\n", elem, max, maxAspectRatio); locFlag = PETSC_TRUE; } else { PetscLogInfo(mesh, "Maximum aspect ratio in element %d: %g\n", elem, max); } PetscLogFlops(elem*19); ierr = MPI_Allreduce(&locFlag, flag, 1, MPI_INT, MPI_LOR, mesh->comm); CHKERRQ(ierr); ierr = MPI_Allreduce(&locRetVal, &retVal, 1, MPI_INT, MPI_LOR, mesh->comm); CHKERRQ(ierr); PetscFunctionReturn(retVal); } static struct _MeshOps MeshOps = {/* Generic Operations */ PETSC_NULL/* MeshSetup */, PETSC_NULL/* MeshSetFromOptions */, MeshView_Triangular_2D, PETSC_NULL/* MeshCopy */, PETSC_NULL/* MeshDuplicate */, MeshDestroy_Triangular_2D, /* Mesh-Specific Operations */ MeshPartition_Triangular_2D, MeshCoarsen_Triangular_2D, MeshRefine_Triangular_2D, MeshResetNodes_Triangular_2D, MeshSaveMesh_Triangular_2D, MeshRestoreMesh_Triangular_2D, /* Mesh Query Functions */ MeshUpdateBoundingBox_Triangular_2D, MeshIsDistorted_Triangular_2D, /* Mesh Boundary Query Functions */ MeshGetBoundarySize_Triangular_2D, MeshGetBoundaryIndex_Triangular_2D, MeshGetBoundaryStart_Triangular_2D, MeshGetBoundaryNext_Triangular_2D, MeshGetActiveBoundary_Triangular_2D, /* Mesh Node Query Functions */ MeshGetNodeBoundary_Triangular_2D, MeshNodeIsVertex_Triangular_2D, MeshGetNodeCoords_Triangular_2D, MeshSetNodeCoords_Triangular_2D, MeshGetNodeCoordsSaved_Triangular_2D, MeshNearestNode_Triangular_2D, MeshGetNodeSupport_Triangular_2D, /* Mesh Element Query Functions */ MeshGetElemNeighbor_Triangular_2D, MeshLocatePoint_Triangular_2D, /* Mesh Embedding Query Functions */ MeshGetNodeFromElement_Triangular_2D, MeshGetNodeFromEdge_Triangular_2D, /* CSR Support Functions */ MeshCreateLocalCSR_Triangular_2D, MeshCreateFullCSR_Triangular_2D, MeshCreateDualCSR_Triangular_2D}; #undef __FUNCT__ #define __FUNCT__ "MeshCreateTriangular2DCSR" /*@ MeshCreateTriangular2DCSR - Creates a basic two dimensional unstructured grid from an existing CSR representation. Collective on MPI_Comm Input Parameters: + comm - The communicator that shares the grid . numNodes - The number of mesh nodes, here identical to vertices . numFaces - The number of elements in the mesh . nodes - The coordinates for each node . offsets - The offset of each node's adjacency list . adj - The adjacency list for the mesh - faces - The list of nodes for each element Output Parameter: . mesh - The mesh Level: beginner .keywords unstructured mesh .seealso MeshCreateTriangular2D(), MeshCreateTriangular2DPeriodic() @*/ int MeshCreateTriangular2DCSR(MPI_Comm comm, int numNodes, int numFaces, double *nodes, int *offsets, int *adj, int *faces, Mesh *mesh) { ParameterDict dict; MeshBoundary2D *bdCtx; int ierr; PetscFunctionBegin; /* REWORK: We are really stretching this interface, but I need this yesterday */ ierr = PetscMalloc(sizeof(MeshBoundary2D), &bdCtx); CHKERRQ(ierr); bdCtx->numBd = 0; bdCtx->numVertices = numNodes; bdCtx->vertices = nodes; bdCtx->markers = offsets; bdCtx->segMarkers = adj; bdCtx->numSegments = numFaces; bdCtx->segments = faces; bdCtx->numHoles = 0; bdCtx->holes = PETSC_NULL; ierr = MeshCreate(comm, mesh); CHKERRQ(ierr); ierr = ParameterDictCreate(comm, &dict); CHKERRQ(ierr); ierr = ParameterDictSetObject(dict, "bdCtx", bdCtx); CHKERRQ(ierr); ierr = ParameterDictSetInteger(dict, "numLocNodes", 3); CHKERRQ(ierr); ierr = PetscObjectSetParameterDict((PetscObject) *mesh, dict); CHKERRQ(ierr); ierr = ParameterDictDestroy(dict); CHKERRQ(ierr); ierr = MeshSetDimension(*mesh, 2); CHKERRQ(ierr); /* Setup special mechanisms */ ierr = PetscObjectComposeFunction((PetscObject) *mesh, "MeshTriangular2D_Create_C", "MeshCreate_CSR", (void (*)(void)) MeshCreate_CSR); CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject) *mesh, "MeshTriangular2D_Refine_C", "MeshRefine_CSR", (void (*)(void)) MeshRefine_CSR); CHKERRQ(ierr); ierr = MeshSetType(*mesh, MESH_TRIANGULAR_2D); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) *mesh, "Constructed Mesh"); CHKERRQ(ierr); ierr = PetscFree(bdCtx); CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "MeshSerialize_Triangular_2D" int MeshSerialize_Triangular_2D(MPI_Comm comm, Mesh *mesh, PetscViewer viewer, PetscTruth store) { Mesh m; Mesh_Triangular *tri; int fd; int zero = 0; int one = 1; int numNodes, numEdges, numFaces, numCorners, numBd, numBdNodes, numBdEdges, old; int ierr; PetscFunctionBegin; ierr = PetscViewerBinaryGetDescriptor(viewer, &fd); CHKERRQ(ierr); if (store) { m = *mesh; ierr = PetscBinaryWrite(fd, &m->highlightElement, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->maxDegree, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->startX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->endX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->sizeX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->startY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->endY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->sizeY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->startZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->endZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->sizeZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locStartX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locEndX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locSizeX, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locStartY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locEndY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locSizeY, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locStartZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locEndZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->locSizeZ, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->isPeriodic, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->isPeriodicDim, 3, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numBd, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numVertices, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numNodes, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numBdNodes, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numEdges, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numBdEdges, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numFaces, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numBdFaces, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numCells, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numCorners, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numCellCorners, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->numHoles, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PartitionSerialize(m, &m->part, viewer, store); CHKERRQ(ierr); ierr = PartitionGetNumOverlapElements(m->part, &numFaces); CHKERRQ(ierr); ierr = PartitionGetNumOverlapNodes(m->part, &numNodes); CHKERRQ(ierr); numEdges = m->numEdges; numCorners = m->numCorners; numBd = m->numBd; numBdNodes = m->numBdNodes; numBdEdges = m->numBdEdges; tri = (Mesh_Triangular *) (*mesh)->data; ierr = PetscBinaryWrite(fd, tri->nodes, numNodes*2, PETSC_DOUBLE, 0); CHKERRQ(ierr); if (tri->nodesOld == PETSC_NULL) { ierr = PetscBinaryWrite(fd, &zero, 1, PETSC_INT, 0); CHKERRQ(ierr); } else { ierr = PetscBinaryWrite(fd, &one, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->nodesOld, numNodes*2, PETSC_DOUBLE, 0); CHKERRQ(ierr); } ierr = PetscBinaryWrite(fd, tri->markers, numNodes, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->degrees, numNodes, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->edges, numEdges*2, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->edgemarkers, numEdges, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->faces, numFaces*numCorners, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->neighbors, numFaces*3, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->bdNodes, numBdNodes, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->bdEdges, numBdEdges, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->bdMarkers, numBd, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->bdBegin, numBd+1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, tri->bdEdgeBegin, numBd+1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, m->holes, m->numHoles*2, PETSC_DOUBLE, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &m->maxAspectRatio, 1, PETSC_DOUBLE, 0); CHKERRQ(ierr); } else { /* Create the mesh context */ ierr = MeshCreate(comm, &m); CHKERRQ(ierr); ierr = PetscNew(Mesh_Triangular, &tri); CHKERRQ(ierr); PetscLogObjectMemory(m, sizeof(Mesh_Triangular)); ierr = PetscMemcpy(m->ops, &MeshOps, sizeof(struct _MeshOps)); CHKERRQ(ierr); ierr = PetscStrallocpy(MESH_TRIANGULAR_2D, &m->type_name); CHKERRQ(ierr); m->data = (void *) tri; m->dim = 2; ierr = PetscBinaryRead(fd, &m->highlightElement, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->maxDegree, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscMalloc(m->maxDegree * sizeof(int), &m->support); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->startX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->endX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->sizeX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->startY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->endY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->sizeY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->startZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->endZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->sizeZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locStartX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locEndX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locSizeX, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locStartY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locEndY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locSizeY, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locStartZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locEndZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->locSizeZ, 1, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->isPeriodic, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->isPeriodicDim, 3, PETSC_INT); CHKERRQ(ierr); m->activeBd = -1; m->activeBdOld = -1; m->activeBdNode = -1; tri->areas = PETSC_NULL; tri->aspectRatios = PETSC_NULL; ierr = PetscBinaryRead(fd, &m->numBd, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numVertices, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numNodes, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numBdNodes, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numEdges, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numBdEdges, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numFaces, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numBdFaces, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numCells, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numCorners, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numCellCorners, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->numHoles, 1, PETSC_INT); CHKERRQ(ierr); ierr = PartitionSerialize(m, &m->part, viewer, store); CHKERRQ(ierr); PetscLogObjectParent(m, m->part); ierr = PartitionGetNumOverlapElements(m->part, &numFaces); CHKERRQ(ierr); ierr = PartitionGetNumOverlapNodes(m->part, &numNodes); CHKERRQ(ierr); numEdges = m->numEdges; numCorners = m->numCorners; numBd = m->numBd; numBdNodes = m->numBdNodes; numBdEdges = m->numBdEdges; ierr = PetscMalloc(numNodes*2 * sizeof(double), &tri->nodes); CHKERRQ(ierr); ierr = PetscMalloc(numNodes * sizeof(int), &tri->markers); CHKERRQ(ierr); ierr = PetscMalloc(numNodes * sizeof(int), &tri->degrees); CHKERRQ(ierr); ierr = PetscMalloc(numEdges*2 * sizeof(int), &tri->edges); CHKERRQ(ierr); ierr = PetscMalloc(numEdges * sizeof(int), &tri->edgemarkers); CHKERRQ(ierr); ierr = PetscMalloc(numFaces*numCorners * sizeof(int), &tri->faces); CHKERRQ(ierr); ierr = PetscMalloc(numFaces*3 * sizeof(int), &tri->neighbors); CHKERRQ(ierr); ierr = PetscMalloc(numBdNodes * sizeof(int), &tri->bdNodes); CHKERRQ(ierr); ierr = PetscMalloc(numBdEdges * sizeof(int), &tri->bdEdges); CHKERRQ(ierr); ierr = PetscMalloc(numBd * sizeof(int), &tri->bdMarkers); CHKERRQ(ierr); ierr = PetscMalloc((numBd+1) * sizeof(int), &tri->bdBegin); CHKERRQ(ierr); ierr = PetscMalloc((numBd+1) * sizeof(int), &tri->bdEdgeBegin); CHKERRQ(ierr); if (m->numHoles > 0) { ierr = PetscMalloc(m->numHoles*2 * sizeof(double), &m->holes); CHKERRQ(ierr); } else { m->holes = PETSC_NULL; } PetscLogObjectMemory(m, (numNodes*2 + m->numHoles*2) * sizeof(double) + (numNodes*2 + numEdges*3 + numFaces*(numCorners+1) + numBdNodes + numBdEdges + numBd*3+2)*sizeof(int)); ierr = PetscBinaryRead(fd, tri->nodes, numNodes*2, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &old, 1, PETSC_INT); CHKERRQ(ierr); if (old) { ierr = PetscMalloc(numNodes*2 * sizeof(double), &tri->nodesOld); CHKERRQ(ierr); PetscLogObjectMemory(m, numNodes*2 * sizeof(double)); ierr = PetscBinaryRead(fd, tri->nodesOld, numNodes*2, PETSC_DOUBLE); CHKERRQ(ierr); } else { tri->nodesOld = PETSC_NULL; } ierr = PetscBinaryRead(fd, tri->markers, numNodes, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->degrees, numNodes, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->edges, numEdges*2, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->edgemarkers, numEdges, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->faces, numFaces*numCorners, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->neighbors, numFaces*3, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->bdNodes, numBdNodes, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->bdEdges, numBdEdges, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->bdMarkers, numBd, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->bdBegin, numBd+1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, tri->bdEdgeBegin, numBd+1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, m->holes, m->numHoles*m->dim, PETSC_DOUBLE); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &m->maxAspectRatio, 1, PETSC_DOUBLE); CHKERRQ(ierr); #ifdef PETSC_USE_BOPT_g /* Check mesh integrity */ ierr = MeshDebug_Triangular_2D(m, PETSC_TRUE); CHKERRQ(ierr); #endif *mesh = m; } PetscFunctionReturn(0); } EXTERN_C_END EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "MeshCreate_Triangular_2D" int MeshCreate_Triangular_2D(Mesh mesh) { Mesh_Triangular *tri; int (*f)(MeshBoundary2D *, int, Mesh); MPI_Comm comm; PetscTruth opt; int ierr; PetscFunctionBegin; ierr = PetscNew(Mesh_Triangular, &tri); CHKERRQ(ierr); ierr = PetscMemcpy(mesh->ops, &MeshOps, sizeof(struct _MeshOps)); CHKERRQ(ierr); mesh->data = (void *) tri; ierr = PetscStrallocpy(MESH_SER_TRIANGULAR_2D_BINARY, &mesh->serialize_name); CHKERRQ(ierr); mesh->dim = 2; /* Setup the periodic domain */ if (mesh->isPeriodic == PETSC_TRUE) { ierr = MeshBoundary2DSetBoundingBox(mesh, mesh->bdCtx); CHKERRQ(ierr); } /* Setup default mechanisms */ ierr = PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", (PetscVoidFunction) &f); CHKERRQ(ierr); if (f == PETSC_NULL) { #ifdef PETSC_HAVE_TRIANGLE if (mesh->isPeriodic == PETSC_TRUE) { ierr = PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", "MeshCreate_Periodic", (void (*)(void)) MeshCreate_Periodic); CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", "MeshRefine_Periodic", (void (*)(void)) MeshRefine_Periodic); CHKERRQ(ierr); } else { ierr = PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", "MeshCreate_Triangle", (void (*)(void)) MeshCreate_Triangle); CHKERRQ(ierr); ierr = PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", "MeshRefine_Triangle", (void (*)(void)) MeshRefine_Triangle); } #endif } ierr = MeshCheckBoundary_Triangular_2D(mesh); CHKERRQ(ierr); ierr = PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", (PetscVoidFunction) &f); CHKERRQ(ierr); if (!f) SETERRQ(1,"Unable to find mesh generation function"); ierr = (*f)(mesh->bdCtx, mesh->numCorners, mesh); CHKERRQ(ierr); /* Calculate maximum degree of vertices */ ierr = MeshSetupSupport_Triangular_2D(mesh); CHKERRQ(ierr); /* Construct derived and boundary information */ ierr = MeshSetupBoundary_Triangular_2D(mesh); CHKERRQ(ierr); /* Reorder nodes before distributing mesh */ ierr = PetscOptionsHasName(mesh->prefix, "-mesh_reorder", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { /* MUST FIX: Since we have duplicated the whole mesh, we may impersonate a serial mesh */ comm = mesh->comm; mesh->comm = PETSC_COMM_SELF; ierr = MeshGetOrdering(mesh, MESH_ORDER_TRIANGULAR_2D_RCM, &mesh->nodeOrdering); CHKERRQ(ierr); mesh->comm = comm; /* Permute arrays implicitly numbered by node numbers */ ierr = AOApplicationToPetscPermuteReal(mesh->nodeOrdering, mesh->dim, tri->nodes); CHKERRQ(ierr); ierr = AOApplicationToPetscPermuteInt(mesh->nodeOrdering, 1, tri->markers); CHKERRQ(ierr); ierr = AOApplicationToPetscPermuteInt(mesh->nodeOrdering, 1, tri->degrees); CHKERRQ(ierr); /* Renumber arrays dependent on the canonical node numbering */ ierr = AOApplicationToPetsc(mesh->nodeOrdering, mesh->numEdges*2, tri->edges); CHKERRQ(ierr); ierr = AOApplicationToPetsc(mesh->nodeOrdering, mesh->numFaces*mesh->numCorners, tri->faces); CHKERRQ(ierr); ierr = AOApplicationToPetsc(mesh->nodeOrdering, mesh->numBdNodes, tri->bdNodes); CHKERRQ(ierr); } /* Initialize partition */ ierr = MeshPartition(mesh); CHKERRQ(ierr); /* Initialize save space */ tri->nodesOld = PETSC_NULL; /* Reset the midnodes */ if (mesh->isPeriodic == PETSC_TRUE) { ierr = MeshResetNodes(mesh, PETSC_TRUE); CHKERRQ(ierr); } PetscFunctionReturn(0); } EXTERN_C_END