xmesh2d.py 8.04 KB
Newer Older
dazerki's avatar
dazerki committed
1
2
3
4
5
6
7
8
import gmsh
import numpy as np
import itertools
gmsh.initialize()


class TrackerUI:

dazerki's avatar
dazerki committed
9
    def  __init__(self, initial_dt=0.05):
dazerki's avatar
dazerki committed
10
11
12
13
14
15
16
17
18
19
20
21
        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"""[
dazerki's avatar
dazerki committed
22
                {{ "type":"number", "name":"Pause", "values":[0], "choices":[0, 1] }},
dazerki's avatar
dazerki committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
                {{ "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

dazerki's avatar
dazerki committed
112
113
114
115
116
117
118
119
    @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])
tleyssens6's avatar
tleyssens6 committed
120
        tags, x, _ = gmsh.model.mesh.get_nodes(includeBoundary=False)
dazerki's avatar
dazerki committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
        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


dazerki's avatar
dazerki committed
140
141
142
143
144
145
146
147
148
149

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):
tleyssens6's avatar
tleyssens6 committed
150
        print("in xmesh")
dazerki's avatar
dazerki committed
151
152
153
154
155
156
157
158
159
160
161
162
163

        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]
dazerki's avatar
dazerki committed
164
165
166
167
168
                # if alpha>0.95:
                #     print("HERE ")
                #     return True
                # else:
                #     return False
dazerki's avatar
dazerki committed
169
170
171
172
173
                return alpha == 1
            else:
                return False

        neigh = self.mesh.neigh
tleyssens6's avatar
tleyssens6 committed
174
        boundary_nodes = self.mesh.boundary_nodes
dazerki's avatar
dazerki committed
175
176
177
178
179
180
181
182
183
184
185
        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]
tleyssens6's avatar
tleyssens6 committed
186
            print(nsign)
dazerki's avatar
dazerki committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
            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