Source code for UM2N.generator.polymesh

import random

import firedrake as fd
import gmsh
import numpy as np

__all__ = ["RandPolyMesh"]


[docs] class RandPolyMesh: """ Create a random polygonal mesh by spliting the edge of a square randomly. """ def __init__(self, scale=1.0, mesh_type=2): # params setup self.mesh_type = mesh_type self.scale = scale self.start = 0 self.end = self.scale self.split_threshold = 0.3 self.mid = (self.start + self.end) / 2 self.quater = (self.start + self.mid) / 2 self.three_quater = (self.mid + self.end) / 2 self.mid_interval = (self.end - self.start) / 3 self.quater_interval = (self.mid - self.start) / 4 # temp vars self.points = [] self.lines = [] # generate mesh self.get_rand_points() return
[docs] def get_mesh(self, res=1e-1, file_path="./temp.msh"): gmsh.initialize() gmsh.model.add("t1") # params setup self.lc = res self.start = 0 self.end = self.scale self.mid = (self.start + self.end) / 2 self.quater = (self.start + self.mid) / 2 self.three_quater = (self.mid + self.end) / 2 self.mid_interval = (self.end - self.start) / 3 self.quater_interval = (self.mid - self.start) / 4 self.file_path = file_path # temp vars self.points = [] self.lines = [] # generate mesh self.get_points() self.get_line() self.get_curve() self.get_plane() gmsh.model.geo.synchronize() gmsh.option.setNumber("Mesh.Algorithm", self.mesh_type) self.get_boundaries() gmsh.model.addPhysicalGroup(2, [1], name="My surface") gmsh.model.mesh.generate(2) gmsh.write(self.file_path) gmsh.finalize() self.num_boundary = len(self.lines) return fd.Mesh(self.file_path)
[docs] def get_rand(self, mean, interval): return random.uniform(mean - interval, mean + interval)
[docs] def get_rand_points(self): points = [] split_p = np.random.uniform(0, 1, 4) # edge 1 if split_p[0] < self.split_threshold: points.append([self.get_rand(self.quater, self.quater_interval), 0]) points.append([self.get_rand(self.three_quater, self.quater_interval), 0]) else: points.append([self.get_rand(self.mid, self.mid_interval), 0]) # edge 2 if split_p[1] < self.split_threshold: points.append( [self.scale, self.get_rand(self.quater, self.quater_interval)] ) points.append( [self.scale, self.get_rand(self.three_quater, self.quater_interval)] ) else: points.append([self.scale, self.get_rand(self.mid, self.mid_interval)]) # edge 3 if split_p[2] < self.split_threshold: points.append( [self.get_rand(self.three_quater, self.quater_interval), self.scale] ) points.append( [self.get_rand(self.quater, self.quater_interval), self.scale] ) else: points.append([self.get_rand(self.mid, self.mid_interval), self.scale]) # edge 4 if split_p[3] < self.split_threshold: points.append([0, self.get_rand(self.three_quater, self.quater_interval)]) points.append([0, self.get_rand(self.quater, self.quater_interval)]) else: points.append([0, self.get_rand(self.mid, self.mid_interval)]) # points.append(p1) self.raw_points = points return
[docs] def get_points(self): temp = [] for i in range(len(self.raw_points)): temp.append( gmsh.model.geo.addPoint( self.raw_points[i][0], self.raw_points[i][1], 0, self.lc ) ) self.points = temp
[docs] def get_line(self): for i in range(len(self.points)): if i < len(self.points) - 1: line = gmsh.model.geo.addLine(self.points[i], self.points[i + 1]) self.lines.append(line) else: line = gmsh.model.geo.addLine(self.points[i], self.points[0]) self.lines.append(line) return
[docs] def get_boundaries(self): print("in get_boundaries lines:", self.lines) for i, line_tag in enumerate(self.lines): gmsh.model.addPhysicalGroup(1, [line_tag], i + 1) gmsh.model.setPhysicalName(1, i + 1, "Boundary " + str(i + 1))
[docs] def get_curve(self): gmsh.model.geo.addCurveLoop([i for i in range(1, len(self.points) + 1)], 1)
[docs] def get_plane(self): gmsh.model.geo.addPlaneSurface([1], 1)
[docs] def show(self, file_path): mesh = fd.Mesh(file_path) fig = fd.triplot(mesh) return fig
if __name__ == "__main__": import matplotlib.pyplot as plt mesh_gen = RandPolyMesh(mesh_type=2) mesh_coarse = mesh_gen.get_mesh(res=5e-2, file_path="./temp1.msh") mesh_fine = mesh_gen.get_mesh(res=4e-2, file_path="./temp2.msh") mesh_gen.show("./temp1.msh") mesh_gen.show("./temp2.msh") plt.show()