Demonstration of using shapely to create FiPy mesh and check point containment in different domains

LICENSE.md

Terms of Use

This software was developed by employees of the National Institute of Standards and Technology NIST, an agency of the Federal Government and is being made available as a public service. Pursuant to title 17 United States Code Section 105, works of NIST employees are not subject to copyright protection in the United States. This software may be subject to foreign copyright. Permission in the United States and in foreign countries, to the extent that NIST may hold copyright, to use, copy, modify, create derivative works, and distribute this software and its documentation without fee is hereby granted on a non-exclusive basis, provided that this notice and disclaimer of warranty appears in all copies.

THE SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY THAT THE SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND FREEDOM FROM INFRINGEMENT, AND ANY WARRANTY THAT THE DOCUMENTATION WILL CONFORM TO THE SOFTWARE, OR ANY WARRANTY THAT THE SOFTWARE WILL BE ERROR FREE. IN NO EVENT SHALL NIST BE LIABLE FOR ANY DAMAGES, INCLUDING, BUT NOT LIMITED TO, DIRECT, INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES, ARISING OUT OF, RESULTING FROM, OR IN ANY WAY CONNECTED WITH THIS SOFTWARE, WHETHER OR NOT BASED UPON WARRANTY, CONTRACT, TORT, OR OTHERWISE, WHETHER OR NOT INJURY WAS SUSTAINED BY PERSONS OR PROPERTY OR OTHERWISE, AND WHETHER OR NOT LOSS WAS SUSTAINED FROM, OR AROSE OUT OF THE RESULTS OF, OR USE OF, THE SOFTWARE OR SERVICES PROVIDED HEREUNDER.

MeshContainment.ipynb

shapely2geo.py

from shapely.geometry import LineString, Point, MultiPoint, MultiLineString

class Shapely(dict):
    def __init__(self, arg, cellsize=1e-7):
        dict.__init__(self, arg)
        
        self.points = MultiPoint()
        self.lines = MultiLineString()
        self.loops = []
        self.surfaces = []
        self.physicalSurfaces = {}
        self.physicalLines = {}
        self.cellsize = cellsize
        
    def _lineString2lines(self, line):
        lines = []
        coords = list(line.coords)
        for p0, p1 in zip(coords[:-1], coords[1:]):
            l0 = LineString((p0, p1))
            if self.lines.contains(l0):
                for i, l1 in enumerate(self.lines):
                    if l0.equals(l1):
                        # is there no easier way to tell if two lines are anti-parallel?
                        if l0.parallel_offset(l0.length, 'left').equals(l1.parallel_offset(l1.length, 'left')):
                            lines.append(i+1)
                        else:
                            lines.append(-(i+1))
                        break
            else:
                self.lines = MultiLineString(list(self.lines) + [l0])
                lines.append(len(self.lines))
                
        return lines

    def _polygon2geo(self, polygon):
        surface = []

        self.loops.append(self._lineString2lines(polygon.exterior))
        surface.append(len(self.loops))
        for hole in polygon.interiors:
            self.loops.append(self._lineString2lines(polygon.exterior))
            surface.append(len(self.loops))
            
        return surface
        
    @property
    def geo(self):
        for name, patch in self.items():
            if patch.geom_type in ['Polygon']:
                surface = self._polygon2geo(patch)
                            
                self.surfaces.append(surface)
                self.physicalSurfaces[name] = [len(self.surfaces)]
            elif patch.geom_type in ['MultiPolygon']:
                physical = []
                for polygon in patch.geoms:
                    surface = self._polygon2geo(polygon)
                                
                    self.surfaces.append(surface)
                    physical.append(len(self.surfaces))
                    
                self.physicalSurfaces[name] = physical
            elif patch.geom_type in ['LineString']:
                self.physicalLines[name] = self._lineString2lines(patch)
            elif patch.geom_type in ['MultiLineString']:
                physical = []
                for line in patch.geoms:
                    physical += self._lineString2lines(line)
                    
                self.physicalLines[name] = physical
                

                
        lines = []
        for line in self.lines:
            indices = []
            for p0 in line.boundary:
                if self.points.contains(p0):
                    for i, p1 in enumerate(self.points):
                        if p0.equals(p1):
                            indices.append(i+1)
                else:
                    self.points = MultiPoint(list(self.points) + [p0])
                    indices.append(len(self.points))

            lines.append(indices)
                
        geo = ["\n".join(["Point(%d) = {%g, %g, 0.0, %g};" % (i+1, p.x, p.y, self.cellsize) for i, p in enumerate(list(self.points))]),
               "\n".join(["Line(%d) = {%d, %d};" % (i+1, p0, p1) for i, (p0, p1) in enumerate(lines)]),
               "\n".join(["Line Loop(%d) = {%s};" % (i+1, ", ".join([str(l) for l in loop])) for i, loop in enumerate(self.loops)]),
               "\n".join(["Plane Surface(%d) = {%s};" % (i+1, ", ".join([str(l) for l in loop])) for i, loop in enumerate(self.surfaces)]),
               "\n".join(["Physical Surface(\"%s\") = {%s};" % (key, ", ".join([str(s) for s in surfaces])) for key, surfaces in self.physicalSurfaces.items()]),
               "\n".join(["Physical Line(\"%s\") = {%s};" % (key, ", ".join([str(l) for l in lines])) for key, lines in self.physicalLines.items()])]
               
        return "\n\n".join(geo) + "\n"