From 33312a75949fc63c88ab8eb50a49f7a408eab689 Mon Sep 17 00:00:00 2001 From: "Jonathan E. Guyer" Date: Thu, 9 Aug 2012 20:30:38 +0000 Subject: [PATCH 1/4] creating branch for ticket452 git-svn-id: svn+ssh://code.matforge.org/fipy/branches/ticket452@5283 d80e17d7-ff13-0410-a124-85740d801063 From cce069526101a39d10e736fa18c00e8af4cc70f5 Mon Sep 17 00:00:00 2001 From: "Jonathan E. Guyer" Date: Fri, 10 Aug 2012 17:48:18 +0000 Subject: [PATCH 2/4] `sumAll` is (and always was) silly git-svn-id: svn+ssh://code.matforge.org/fipy/branches/ticket452@5284 d80e17d7-ff13-0410-a124-85740d801063 --- fipy/meshes/gmshImport.py | 4 ++-- fipy/tools/comms/commWrapper.py | 5 +---- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/fipy/meshes/gmshImport.py b/fipy/meshes/gmshImport.py index 5abb18e42b..02693dc130 100755 --- a/fipy/meshes/gmshImport.py +++ b/fipy/meshes/gmshImport.py @@ -1647,7 +1647,7 @@ def __init__(self, self.mshFile.close() if communicator.Nproc > 1: - self.globalNumberOfCells = communicator.sumAll(len(self.cellGlobalIDs)) + self.globalNumberOfCells = communicator.sum(len(self.cellGlobalIDs)) parprint(" I'm solving with %d cells total." % self.globalNumberOfCells) parprint(" Got global number of cells") @@ -1982,7 +1982,7 @@ def __init__(self, arg, communicator=parallel, order=1, background=None): _TopologyClass=_GmshTopology) if self.communicator.Nproc > 1: - self.globalNumberOfCells = self.communicator.sumAll(len(self.cellGlobalIDs)) + self.globalNumberOfCells = self.communicator.sum(len(self.cellGlobalIDs)) (self.physicalCellMap, self.geometricalCellMap, diff --git a/fipy/tools/comms/commWrapper.py b/fipy/tools/comms/commWrapper.py index a8d102e053..fad360113b 100644 --- a/fipy/tools/comms/commWrapper.py +++ b/fipy/tools/comms/commWrapper.py @@ -87,11 +87,8 @@ def allgather(self, sendobj=None, recvobj=None): return recvobj - def sumAll(self, a): - return self.epetra_comm.SumAll(numerix.array(a)) - def sum(self, a, axis=None): - return self.epetra_comm.SumAll(numerix.array(a.sum(axis=axis))) + return self.epetra_comm.SumAll(numerix.array(a).sum(axis=axis)) def __getstate__(self): return {'dummy': 0} From 8055286868523e4624ead7c5a7fea6e2ef349cf3 Mon Sep 17 00:00:00 2001 From: "Jonathan E. Guyer" Date: Mon, 13 Aug 2012 13:31:15 +0000 Subject: [PATCH 3/4] * write POSFile in parallel * remove cruft git-svn-id: svn+ssh://code.matforge.org/fipy/branches/ticket452@5285 d80e17d7-ff13-0410-a124-85740d801063 --- fipy/meshes/gmshImport.py | 202 +++++++++++++------------------------- 1 file changed, 67 insertions(+), 135 deletions(-) diff --git a/fipy/meshes/gmshImport.py b/fipy/meshes/gmshImport.py index 02693dc130..bbbc5ee80e 100755 --- a/fipy/meshes/gmshImport.py +++ b/fipy/meshes/gmshImport.py @@ -56,6 +56,8 @@ from fipy.meshes.mesh2D import Mesh2D from fipy.meshes.topologies.meshTopology import _MeshTopology +from fipy.tools.debug import PRINT + __all__ = ["openMSHFile", "openPOSFile", "Gmsh2D", "Gmsh2DIn3DSpace", "Gmsh3D", "GmshGrid2D", "GmshGrid3D", @@ -203,16 +205,20 @@ def openMSHFile(name, dimensions=None, coordDimensions=None, communicator=parall gmshFlags += ["-format", "msh"] - if communicator.procID == 0: - if background is not None: + if background is not None: + if communicator.procID == 0: f, bgmf = tempfile.mkstemp(suffix=".pos") os.close(f) - f = openPOSFile(name=bgmf, mode='w') - f.write(background) - f.close() - - gmshFlags += ["-bgm", bgmf] + else: + bgmf = None + bgmf = communicator.bcast(bgmf) + f = openPOSFile(name=bgmf, mode='w') + f.write(background) + f.close() + + gmshFlags += ["-bgm", bgmf] + if communicator.procID == 0: (f, mshFile) = tempfile.mkstemp('.msh') os.close(f) fileIsTemporary = True @@ -380,30 +386,33 @@ def _writeValues(self, var, dimensions, time=0.0, timeindex=0): nb_scalar_pyramids2 = nb_vector_pyramids2 = nb_tensor_pyramids2 = 0 nb_text2d = nb_text2d_chars = nb_text3d = nb_text3d_chars = 0 + nonOverlapping = nx.zeros((var.mesh.numberOfCells,), dtype=bool) + nonOverlapping[..., var.mesh._localNonOverlappingCellIDs] = True + cellTopology = var.mesh._cellTopology t = var.mesh.topology._elementTopology if var.rank == 0: - nb_scalar_triangles = (cellTopology == t["triangle"]).sum() - nb_scalar_quadrangles = (cellTopology == t["quadrangle"]).sum() - nb_scalar_tetrahedra = (cellTopology == t["tetrahedron"]).sum() - nb_scalar_hexahedra = (cellTopology == t["hexahedron"]).sum() - nb_scalar_prisms = (cellTopology == t["prism"]).sum() - nb_scalar_pyramids = (cellTopology == t["pyramid"]).sum() + nb_scalar_triangles = self.communicator.sum((cellTopology == t["triangle"]) & nonOverlapping) + nb_scalar_quadrangles = self.communicator.sum((cellTopology == t["quadrangle"]) & nonOverlapping) + nb_scalar_tetrahedra = self.communicator.sum((cellTopology == t["tetrahedron"]) & nonOverlapping) + nb_scalar_hexahedra = self.communicator.sum((cellTopology == t["hexahedron"]) & nonOverlapping) + nb_scalar_prisms = self.communicator.sum((cellTopology == t["prism"]) & nonOverlapping) + nb_scalar_pyramids = self.communicator.sum((cellTopology == t["pyramid"]) & nonOverlapping) elif var.rank == 1: - nb_vector_triangles = (cellTopology == t["triangle"]).sum() - nb_vector_quadrangles = (cellTopology == t["quadrangle"]).sum() - nb_vector_tetrahedra = (cellTopology == t["tetrahedron"]).sum() - nb_vector_hexahedra = (cellTopology == t["hexahedron"]).sum() - nb_vector_prisms = (cellTopology == t["prism"]).sum() - nb_vector_pyramids = (cellTopology == t["pyramid"]).sum() + nb_vector_triangles = self.communicator.sum((cellTopology == t["triangle"]) & nonOverlapping) + nb_vector_quadrangles = self.communicator.sum((cellTopology == t["quadrangle"]) & nonOverlapping) + nb_vector_tetrahedra = self.communicator.sum((cellTopology == t["tetrahedron"]) & nonOverlapping) + nb_vector_hexahedra = self.communicator.sum((cellTopology == t["hexahedron"]) & nonOverlapping) + nb_vector_prisms = self.communicator.sum((cellTopology == t["prism"]) & nonOverlapping) + nb_vector_pyramids = self.communicator.sum((cellTopology == t["pyramid"]) & nonOverlapping) elif var.rank == 2: - nb_tensor_triangles = (cellTopology == t["triangle"]).sum() - nb_tensor_quadrangles = (cellTopology == t["quadrangle"]).sum() - nb_tensor_tetrahedra = (cellTopology == t["tetrahedron"]).sum() - nb_tensor_hexahedra = (cellTopology == t["hexahedron"]).sum() - nb_tensor_prisms = (cellTopology == t["prism"]).sum() - nb_tensor_pyramids = (cellTopology == t["pyramid"]).sum() + nb_tensor_triangles = self.communicator.sum((cellTopology == t["triangle"]) & nonOverlapping) + nb_tensor_quadrangles = self.communicator.sum((cellTopology == t["quadrangle"]) & nonOverlapping) + nb_tensor_tetrahedra = self.communicator.sum((cellTopology == t["tetrahedron"]) & nonOverlapping) + nb_tensor_hexahedra = self.communicator.sum((cellTopology == t["hexahedron"]) & nonOverlapping) + nb_tensor_prisms = self.communicator.sum((cellTopology == t["prism"]) & nonOverlapping) + nb_tensor_pyramids = self.communicator.sum((cellTopology == t["pyramid"]) & nonOverlapping) else: raise ValueError("rank must be 0, 1, or 2") @@ -436,31 +445,23 @@ def _writeValues(self, var, dimensions, time=0.0, timeindex=0): cellVertexIDs = var.mesh._orderedCellVertexIDs for shape in ("triangle", "quadrangle", "tetrahedron", "hexahedron", "prism", "pyramid"): - for i in (cellTopology == t[shape]).nonzero()[0]: - nodes = cellVertexIDs[..., i] - self._writeNodesAndValues(vertexCoords=vertexCoords, - nodes=nodes.compressed(), - value=value[..., i]) - -# cellFaceVertices = nx.take(var.mesh.faceVertexIDs, var.mesh.cellFaceIDs, axis=1) -# faceOrder = nx.argsort(var.mesh._unsortedNodesPerFace, axis=0)[::-1] -# faceOrientations = nx.take(var.mesh._rightHandOrientation, var.mesh.cellFaceIDs, axis=0) -# -# for shape, writeShape in (("triangle", self._writeTriangle), -# ("quadrangle", self._writeQuadrangle), -# ("tetrahedron", self._writeTetrahedron), -# ("hexahedron", self._writeHexahedron), -# ("prism", self._writePrism), -# ("pyramid", self._writePyramid)): -# for i in (cellTopology == t[shape]).nonzero()[0]: -# writeShape(vertexCoords=vertexCoords, -# cell=cellFaceVertices[..., faceOrder[..., i], i], -# orientations=faceOrientations[faceOrder[..., i], i], -# value=value[..., i]) - + for proc in range(self.communicator.Nproc): + # write the elements in sequence to avoid pulling everything + # to one processor. + if proc == self.communicator.procID: +# print self.communicator.procID, (cellTopology == t[shape]), (cellTopology == t[shape]) & nonOverlapping + for i in ((cellTopology == t[shape]) & nonOverlapping).nonzero()[0]: + nodes = cellVertexIDs[..., i] + self._writeNodesAndValues(vertexCoords=vertexCoords, + nodes=nodes.compressed(), + value=value[..., i]) + self.communicator.Barrier() + # need to synchronize all processes at the end of the file + # this is ridiculously difficult to achieve + self.fileobj.seek(self.communicator.MaxAll(self.fileobj.tell())) + self.fileobj.write("$EndView\n") - # lists of double precision numbers giving the node coordinates and the # values associated with the nodes of the nb-scalar-points scalar points, # nb-vector-points vector points, ..., for each of the time-step-values. @@ -476,86 +477,13 @@ def _writeValues(self, var, dimensions, time=0.0, timeindex=0): # comp1-node2-time2 comp2-node2-time2 comp3-node2-time2 # comp1-node3-time2 comp2-node3-time2 comp3-node3-time2 - def _writeTriangle(self, vertexCoords, cell, orientations, value): - # triangle is defined by one face and the remaining point - # from either of the other faces - nodes = cell[..., 0] - nodes = nx.concatenate((nodes, cell[~nx.in1d(cell[..., 1], nodes), 1])) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - - def _writeQuadrangle(self, vertexCoords, cell, orientations, value): - # quadrangle is defined by one face and the opposite face - face0 = cell[..., 0] - for face1 in cell[..., 1:].swapaxes(0, 1): - if not nx.in1d(face1, face0).any(): - break - - # need to ensure face1 is oriented same way as face0 - cross01 = nx.cross((vertexCoords[..., face1[0]] - - vertexCoords[..., face0[0]]), - (vertexCoords[..., face1[1]] - - vertexCoords[..., face0[0]])) - - cross10 = nx.cross((vertexCoords[..., face0[0]] - - vertexCoords[..., face1[0]]), - (vertexCoords[..., face0[1]] - - vertexCoords[..., face1[0]])) - - dim = vertexCoords.shape[0] - if ((dim == 2 and cross01 * cross10 < 0) - or (dim == 3 and nx.dot(cross01, cross10) < 0)): - face1 = face1[::-1] - - nodes = nx.concatenate((face0, face1)) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - - def _writeTetrahedron(self, vertexCoords, cell, orientations, value): - # tetrahedron is defined by one face and the remaining point - # from any of the other faces - nodes = cell[..., 0][::orientations[0]] - nodes = nx.concatenate((nodes, cell[~nx.in1d(cell[..., 1], nodes), 1])) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - - def _writeHexahedron(self, vertexCoords, cell, orientations, value): - # hexahedron is defined by one face and the opposite face - face0 = cell[..., 0][::orientations[0]] - for i, face1 in enumerate(cell[..., 1:].swapaxes(0, 1)): - if not nx.any([node in face1 for node in face0]): - face1 = face1[::orientations[i+1]] - break - - nodes = nx.concatenate((face0, face1)) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - - def _writePrism(self, vertexCoords, cell, orientations, value): - # prism is defined by the two three-sided faces - face0 = cell[..., 3][::orientations[3]] - face1 = cell[..., 4][::orientations[4]] - - nodes = nx.concatenate((face0, face1)) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - - - def _writePyramid(self, vertexCoords, cell, orientations, value): - # pyramid is defined by four-sided face and the remaining point - # from any of the other faces - nodes = cell[..., 0][::orientations[0]] - nodes = nx.concatenate((nodes, cell[~nx.in1d(cell[..., 1], nodes), 1])) - - self._writeNodesAndValues(vertexCoords=vertexCoords, nodes=nodes, value=value) - def _writeNodesAndValues(self, vertexCoords, nodes, value): # strip out masked values nodes = nodes[nodes != -1] numNodes = len(nodes) data = [] dim = vertexCoords.shape[0] - data = [[str(vertexCoords[..., j, nodes[i]]) for i in range(numNodes)] for j in range(dim)] + data = [[str(vertexCoords[..., j, node]) for node in nodes] for j in range(dim)] if dim == 2: data += [["0.0"] * numNodes] data += [[str(value)] * numNodes] @@ -1754,12 +1682,10 @@ def _test(self): True >>> (ftmp, mshFile) = tempfile.mkstemp('.msh') + >>> os.close(ftmp) >>> f = openMSHFile(name=mshFile, mode='w') # doctest: +GMSH >>> f.write(circle) # doctest: +GMSH >>> f.close() # doctest: +GMSH - >>> import sys - >>> if sys.platform == 'win32': - ... os.close(ftmp) >>> from fipy import doctest_raw_input >>> if __name__ == "__main__": @@ -1824,13 +1750,15 @@ def _test(self): >>> from fipy import CellVariable >>> vol = CellVariable(mesh=sqrTri, value=sqrTri.cellVolumes) # doctest: +GMSH - >>> (ftmp, posFile) = tempfile.mkstemp('.pos') + >>> if parallel.procID == 0: + ... (ftmp, posFile) = tempfile.mkstemp('.pos') + ... os.close(ftmp) + ... else: + ... posFile = None + >>> posFile = parallel.bcast(posFile) >>> f = openPOSFile(posFile, mode='w') # doctest: +GMSH >>> f.write(vol) # doctest: +GMSH >>> f.close() # doctest: +GMSH - >>> import sys - >>> if sys.platform == 'win32': - ... os.close(ftmp) >>> f = open(posFile, mode='r') >>> print "".join(f.readlines()) @@ -1869,7 +1797,8 @@ def _test(self): >>> f.close() - >>> os.remove(posFile) + >>> if parallel.procID == 0: + ... os.remove(posFile) """ class Gmsh2DIn3DSpace(Gmsh2D): @@ -2113,13 +2042,15 @@ def _test(self): >>> from fipy import CellVariable >>> vol = CellVariable(mesh=tetPriPyr, value=tetPriPyr.cellVolumes, name="volume") # doctest: +GMSH - >>> (ftmp, posFile) = tempfile.mkstemp('.pos') + >>> if parallel.procID == 0: + ... (ftmp, posFile) = tempfile.mkstemp('.pos') + ... os.close(ftmp) + ... else: + ... posFile = None + >>> posFile = parallel.bcast(posFile) >>> f = openPOSFile(posFile, mode='w') # doctest: +GMSH >>> f.write(vol) # doctest: +GMSH >>> f.close() # doctest: +GMSH - >>> import sys - >>> if sys.platform == 'win32': - ... os.close(ftmp) >>> f = open(posFile, mode='r') # doctest: +GMSH >>> l = f.readlines() # doctest: +GMSH @@ -2175,7 +2106,8 @@ def _test(self): >>> print numerix.allclose(a1, a2) # doctest: +GMSH True - >>> os.remove(posFile) + >>> if parallel.procID == 0: + ... os.remove(posFile) """ class GmshGrid2D(Gmsh2D): From 738e369b29288d8d6754ea9bed6fed1026d6ac0e Mon Sep 17 00:00:00 2001 From: "Jonathan E. Guyer" Date: Mon, 13 Aug 2012 13:57:18 +0000 Subject: [PATCH 4/4] rename `gmshImport.py` to `gmshMesh.py` git-svn-id: svn+ssh://code.matforge.org/fipy/branches/ticket452@5286 d80e17d7-ff13-0410-a124-85740d801063 --- examples/diffusion/circle.py | 2 +- examples/diffusion/circleQuad.py | 2 +- examples/diffusion/steadyState/mesh20x20/modifiedMeshInput.py | 3 ++- fipy/meshes/__init__.py | 4 ++-- fipy/meshes/{gmshImport.py => gmshMesh.py} | 2 +- fipy/meshes/mesh2D.py | 2 +- fipy/meshes/numMesh/gmshImport.py | 2 +- fipy/meshes/test.py | 2 +- fipy/tests/testClass.py | 2 +- 9 files changed, 11 insertions(+), 10 deletions(-) rename fipy/meshes/{gmshImport.py => gmshMesh.py} (99%) diff --git a/examples/diffusion/circle.py b/examples/diffusion/circle.py index 9dd0ebd22b..db73a9888d 100755 --- a/examples/diffusion/circle.py +++ b/examples/diffusion/circle.py @@ -50,7 +50,7 @@ .. _gmsh manual: http://www.geuz.org/gmsh/doc/texinfo/gmsh.html The mesh created by :term:`Gmsh` is then imported into :term:`FiPy` using the -:class:`~fipy.meshes.gmshImport.Gmsh2D` object. +:class:`~fipy.meshes.gmshMesh.Gmsh2D` object. >>> from fipy import * >>> mesh = Gmsh2D(''' diff --git a/examples/diffusion/circleQuad.py b/examples/diffusion/circleQuad.py index fb19a308b2..dd221a88bd 100755 --- a/examples/diffusion/circleQuad.py +++ b/examples/diffusion/circleQuad.py @@ -50,7 +50,7 @@ .. _gmsh manual: http://www.geuz.org/gmsh/doc/texinfo/gmsh.html The mesh created by :term:`Gmsh` is then imported into :term:`FiPy` using the -:class:`~fipy.meshes.gmshImporter.Gmsh2D` object. +:class:`~fipy.meshes.gmshMesh.Gmsh2D` object. >>> from fipy import * >>> mesh = Gmsh2D(''' diff --git a/examples/diffusion/steadyState/mesh20x20/modifiedMeshInput.py b/examples/diffusion/steadyState/mesh20x20/modifiedMeshInput.py index 500f38082c..19880e5473 100755 --- a/examples/diffusion/steadyState/mesh20x20/modifiedMeshInput.py +++ b/examples/diffusion/steadyState/mesh20x20/modifiedMeshInput.py @@ -39,7 +39,8 @@ This input file again solves a 1D diffusion problem as in `./examples/diffusion/steadyState/mesh1D/input.py`. The difference -being that it uses a triangular mesh loaded in using gmshImport. +being that it uses a triangular mesh loaded in using the +:class:`~fipy.meshes.gmshMesh.Gmsh2D` object. The result is again tested in the same way: diff --git a/fipy/meshes/__init__.py b/fipy/meshes/__init__.py index 25d012c7a4..9b639d3c4b 100644 --- a/fipy/meshes/__init__.py +++ b/fipy/meshes/__init__.py @@ -3,7 +3,7 @@ from fipy.meshes.periodicGrid2D import * from fipy.meshes.skewedGrid2D import * from fipy.meshes.tri2D import * -from fipy.meshes.gmshImport import * +from fipy.meshes.gmshMesh import * __all__ = [] __all__.extend(factoryMeshes.__all__) @@ -11,5 +11,5 @@ __all__.extend(periodicGrid2D.__all__) __all__.extend(skewedGrid2D.__all__) __all__.extend(tri2D.__all__) -__all__.extend(gmshImport.__all__) +__all__.extend(gmshMesh.__all__) diff --git a/fipy/meshes/gmshImport.py b/fipy/meshes/gmshMesh.py similarity index 99% rename from fipy/meshes/gmshImport.py rename to fipy/meshes/gmshMesh.py index bbbc5ee80e..4ea38c037b 100755 --- a/fipy/meshes/gmshImport.py +++ b/fipy/meshes/gmshMesh.py @@ -4,7 +4,7 @@ # ################################################################### # FiPy - a finite volume PDE solver in Python # -# FILE: "gmshImport.py" +# FILE: "gmshMesh.py" # # Author: James O'Beirne # Author: Alexander Mont diff --git a/fipy/meshes/mesh2D.py b/fipy/meshes/mesh2D.py index a5b2155c44..1b05cc4e23 100644 --- a/fipy/meshes/mesh2D.py +++ b/fipy/meshes/mesh2D.py @@ -198,7 +198,7 @@ def _extrude(self, mesh, extrudeFunc, layers): ## should extrude cnahe self rather than creating a new mesh? ## the following allows the 2D mesh to be in 3D space, this can be the case for a - ## GmshImporter2DIn3DSpace which would then be extruded. + ## Gmsh2DIn3DSpace which would then be extruded. oldVertices = mesh.vertexCoords if oldVertices.shape[0] == 2: oldVertices = numerix.resize(oldVertices, (3, len(oldVertices[0]))) diff --git a/fipy/meshes/numMesh/gmshImport.py b/fipy/meshes/numMesh/gmshImport.py index dfbe16c4a6..aaf5937ee0 100644 --- a/fipy/meshes/numMesh/gmshImport.py +++ b/fipy/meshes/numMesh/gmshImport.py @@ -1,6 +1,6 @@ from fipy.meshes.numMesh.deprecatedWarning import numMeshDeprecated -from fipy.meshes.gmshImport import * +from fipy.meshes.gmshMesh import * numMeshDeprecated() diff --git a/fipy/meshes/test.py b/fipy/meshes/test.py index 1629481e7a..02e1e8a52e 100755 --- a/fipy/meshes/test.py +++ b/fipy/meshes/test.py @@ -47,7 +47,7 @@ def _suite(): 'fipy.meshes.grid2D', 'fipy.meshes.grid3D', 'fipy.meshes.tri2D', - 'fipy.meshes.gmshImport', + 'fipy.meshes.gmshMesh', 'fipy.meshes.periodicGrid1D', 'fipy.meshes.periodicGrid2D', 'fipy.meshes.uniformGrid1D', diff --git a/fipy/tests/testClass.py b/fipy/tests/testClass.py index 28e62ca01d..0bfb512189 100644 --- a/fipy/tests/testClass.py +++ b/fipy/tests/testClass.py @@ -165,7 +165,7 @@ def printPackageInfo(self): ## Gmsh version try: - from fipy.meshes.gmshImport import gmshVersion + from fipy.meshes.gmshMesh import gmshVersion gmshversion = gmshVersion() if gmshversion is None: print 'gmsh is not installed'