import gmsh import numpy as np import itertools gmsh.initialize() class TrackerUI: def __init__(self, initial_dt=0.05): gmsh.option.set_number("General.Antialiasing",1) gmsh.option.set_number("General.SmallAxes",0) gmsh.option.set_number("Print.Background",1) gmsh.option.set_number("Mesh.LineWidth",2) gmsh.option.set_color("Mesh.One",0,0,0) gmsh.option.set_number("Mesh.Algorithm",2) self.view = gmsh.view.add("front") gmsh.option.set_number(f"View[{self.view}].ShowScale",0) gmsh.option.set_number(f"View[{self.view}].PointType",1) gmsh.option.set_number(f"View[{self.view}].PointSize",5) gmsh.option.set_number(f"View[{self.view}].ColormapBias",-0.15) gmsh.onelab.set(f"""[ {{ "type":"number", "name":"Pause", "values":[0], "choices":[0, 1] }}, {{ "type":"number", "name":"Record", "values":[0], "choices":[0, 1] }}, {{ "type":"number", "name":"dt", "values":[{initial_dt}]}} ]""") self.recid = 0 self.dt = gmsh.onelab.get_number("dt")[0] gmsh.fltk.initialize() def update(self, x, front): for tn in range(x.shape[0]): gmsh.model.mesh.set_node(tn+1, x[tn], [0,0]) self.dt = gmsh.onelab.get_number("dt")[0] self.pause = gmsh.onelab.get_number("Pause") == 1 frontnodes = np.where(front)[0] vs = np.zeros((frontnodes.size, 4)) vs[:,:3] = x[front] gmsh.view.addListData(self.view, "SP", vs.shape[0], vs.reshape([-1])) gmsh.graphics.draw() gmsh.model.set_current(gmsh.model.get_current()) #force redraw if gmsh.onelab.get_number("Record")[0] == 1 and not self.pause: gmsh.write(f"fig-{self.recid:05}.png") self.recid += 1 class FieldUI: def __init__(self, name, **opts): self.view = gmsh.view.add(name) gmsh.option.set_number(f"View[{self.view}].IntervalsType",3) for key, value in opts.items(): gmsh.option.set_number(f"View[{self.view}].{key}",value) def set(self, sol): nodes = np.arange(1, sol.shape[0]+1) gmsh.view.add_model_data(self.view, 0, "", "NodeData", nodes, sol[:,None]) class Mesh(): """ Topological mesh (no coordinates) """ @staticmethod def _create_sub_simplices(elements, subdim): n = elements.shape[1] nsub = subdim+1 comb = list(itertools.combinations(range(n), nsub)) sube = elements[:,comb].reshape(-1,nsub) asort = sube.argsort(axis=1) sube,emap = np.unique(np.take_along_axis(sube,asort,1),axis=0,return_inverse=True) emap = emap.reshape(elements.shape[0],-1) return sube, emap def _gen_neighbours(self): t = self.tri hedges = np.hstack([ [t[:,0],t[:,1]],[t[:,1],t[:,0]], [t[:,1],t[:,2]],[t[:,2],t[:,1]], [t[:,2],t[:,0]],[t[:,0],t[:,2]]]).T-1 hedges = np.unique(hedges,axis=0).astype(np.int64) self.neigh_start = np.cumsum(np.hstack([[0],np.bincount(hedges[:,0])])) self.neigh = hedges[:,1].copy() @classmethod def import_from_gmsh(cls): self = Mesh() gmsh.model.mesh.renumber_nodes() model = gmsh.model.get_current() self.tritags, tri = gmsh.model.mesh.get_elements_by_type(2) self.tri = tri.reshape([-1,3]) tags, x, _ = gmsh.model.mesh.get_nodes() order = np.argsort(tags) self.nodetags = tags[order] x = x.reshape([-1,3])[order] self.num_nodes = x.shape[0] assert(np.all(self.nodetags == np.arange(1, tags.size+1))) self._gen_neighbours() self.boundary_nodes = {} for gdim, gtag in gmsh.model.getPhysicalGroups(1): gname = gmsh.model.get_physical_name(gdim, gtag) if gname is not None: nodes = [] for gent in gmsh.model.get_entities_for_physical_group(gdim, gtag): nodes.append(gmsh.model.mesh.getNodes(gdim, gent, True)[0]) nodes = np.unique(np.hstack(nodes)) self.boundary_nodes[gname] = nodes self.edges, self.tri_edges = Mesh._create_sub_simplices(self.tri, 1) return self, x @classmethod def import_from_file(cls, filename): self = Mesh() gmsh.open(filename) gmsh.model.mesh.renumber_nodes() model = gmsh.model.get_current() self.tritags, tri = gmsh.model.mesh.get_elements_by_type(2) self.tri = tri.reshape([-1,3]) tags, x, _ = gmsh.model.mesh.get_nodes(includeBoundary=False) order = np.argsort(tags) self.nodetags = tags[order] x = x.reshape([-1,3])[order] self.num_nodes = x.shape[0] assert(np.all(self.nodetags == np.arange(1, tags.size+1))) self._gen_neighbours() self.boundary_nodes = {} for gdim, gtag in gmsh.model.getPhysicalGroups(1): gname = gmsh.model.get_physical_name(gdim, gtag) if gname is not None: nodes = [] for gent in gmsh.model.get_entities_for_physical_group(gdim, gtag): nodes.append(gmsh.model.mesh.getNodes(gdim, gent, True)[0]) nodes = np.unique(np.hstack(nodes)) self.boundary_nodes[gname] = nodes self.edges, self.tri_edges = Mesh._create_sub_simplices(self.tri, 1) return self, x class FrontTracker(): def __init__(self, mesh,x): self.mesh = mesh self.front = np.full(x.shape[0], False) self.x0 = x.copy() self.previous_sign = np.full(x.shape[0], 1) def move_front(self, x, levelset): print("in xmesh") def move_node(tn) : all_neighbours = neigh[neigh_start[tn]:neigh_start[tn+1]] neighbours = all_neighbours[ np.logical_and( sign[all_neighbours] == -sign[tn], np.logical_not(front[all_neighbours]))] if len(neighbours) != 0 : dist = np.linalg.norm(x[neighbours]- x[tn], axis=1) grad = (levelset[neighbours]-levelset[tn])/dist target = neighbours[np.argmin(sign[tn]*grad)] alpha = np.clip((levelset[tn])/(levelset[tn]-levelset[target]),0,1) xn[tn] = alpha*x[target]+(1-alpha)*x[tn] # if alpha>0.95: # print("HERE ") # return True # else: # return False return alpha == 1 else: return False neigh = self.mesh.neigh boundary_nodes = self.mesh.boundary_nodes neigh_start = self.mesh.neigh_start front = self.front sign = np.where(levelset > 0, 1, -1) xn = x.copy() # remove from the front all nodes connected to a single side for tn in range(x.shape[0]): if not front[tn]: continue all_neighbours = neigh[neigh_start[tn]:neigh_start[tn+1]] nsign = sign[all_neighbours] nsign[front[all_neighbours]] = 0 nsign = np.r_[sign[tn],nsign] print(nsign) if nsign.max() - nsign.min() != 2 : front[tn] = False # move the nodes already in the front torm = [] for tn in range(x.shape[0]): if front[tn]: if move_node(tn): torm.append(tn) # add to the front the nodes whose sign changed for tn in range(x.shape[0]): if sign[tn] != self.previous_sign[tn] and not(front[tn]): if move_node(tn): torm.append(tn) front[tn] = True # remove from the front the nodes that hit another node # in an attempt to avoid jumping to another front nodes front[torm] = False self.previous_sign[:] = sign return xn def relax(self, x, alpha): xn = x.copy() for tn in range(x.shape[0]): if (not self.front[tn]): xn[tn] = x[tn]*(1-alpha)+self.x0[tn]*alpha return xn