netgen.py 54.1 KB
Newer Older
1
2
3
4
5
# generation and planar cutting of periodic pore networks
#
# Daniel W. Meyer
# Institute of Fluid Dynamics, ETH Zurich
# January 2019
sfritschi's avatar
sfritschi committed
6
import multiprocessing as mp
7
8
9
10
11
12
13
14
from typing import List, Dict, Set, Tuple # for type hints in argument lists

LABELS = ('', 'in', 'out', 'cut') # don't change the order

class Pore:
    """Class of a pore connected to other pores via throats."""
    id = 0
    def __init__(self, pos: List[float], r: float, label: str = LABELS[0], \
15
        throats: Set = None, id: int = -1, index: int = -1, originIndex: int = -1):
16
17
18
19
20
21
22
23
24
        if id == -1: # no id given
            self.id = Pore.id; Pore.id = Pore.id+1
        else: # id provided
            self.id = id; Pore.id = max(id+1, Pore.id)
        self.pos = pos.copy() # position vector
        if throats is None: throats = set()
        self.throats = throats # throats set
        self.r = r # radius
        self.label = label # label string like '', 'in', or 'out'
sfritschi's avatar
sfritschi committed
25
        self.index = index
26
        self.originIndex = originIndex  # points to original pore inside domain
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    def __repr__(self): return str(self.__class__) + ': ' + \
        str({k:self.__dict__[k] for k in self.__dict__ if k != 'throats'}) + \
        ' {0:d} throats'.format(len(self.throats)) # dont flush throat objects

class Throat:
    """Class of a throat connecting two pores.
    
    Periodic throats have a label 'X1 X2 X3' with Xc being element of {-1,0,1}.
    For Xc = 1, pore1 is at the right or upper domain bound in the c-direction
    and pore2 is at the left or lower bound. Vice versa for Xc = -1. With
    Xc = 0, the throat does not cross the bound periodically in c-direction.
    """
    id = 0
    def __init__(self, pore1: Pore, pore2: Pore, r: float, \
41
        label: str = LABELS[0], id: int = -1, index: int = -1):
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
        if id == -1: # no id given
            self.id = Throat.id; Throat.id = Throat.id+1
        else: # id provided
            self.id = id; Throat.id = max(id+1, Throat.id)
        # pore objects
        self.pore1 = pore1
        self.pore2 = pore2
        # radius and label
        self.r = r
        self.label = label # label string like '' or periodicity label '1 0 -1'
    def __repr__(self): return str(self.__class__) + ': ' + str(self.__dict__)

class Network:
    """Network class with pore and throat lists.
    
    Lmax must be smaller than the smallest side length of the network ub-lb.
    """
    def __init__(self, lb: List[float] = [], ub: List[float] = [], \
        pores: List[Pore] = None, throats: List[Throat] = None, \
sfritschi's avatar
sfritschi committed
61
        Lmax: float = 0.0, label: str = None):
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
        # lower/upper bounds for network volume
        self.lb = lb.copy(); self.ub = ub.copy()
        # sets with pores and throats
        if pores is None: pores = set()
        if throats is None: throats = set()
        self.pores = pores.copy(); self.throats = throats.copy()
        # length of longest throat (distance between connected pores)
        self.Lmax = Lmax
        # network label or name
        from datetime import datetime
        if label is None:
            self.label = datetime.today().isoformat(sep=' ',timespec='seconds')
        else:
            self.label = label
    def __repr__(self):
        # count in/out resp. labelled pores
        k = sum([int(pore.label != LABELS[0]) for pore in self.pores])
        # report
        return 'Network \'' + self.label + '\' from ' + str(self.lb) + \
            ' to ' + str(self.ub) + \
            ' with {0:d}({1:d}) pores, {2:d} throats, Lmax = {3:e}'. \
            format(len(self.pores), k, len(self.throats), self.Lmax)
    def add_pore(self, pore: Pore):
        self.pores.add(pore)
    def remove_pore(self, pore: Pore):
        """Remove pore and all throats connected to it."""
        # disconnect throats from pores connected to pore
        for throat in pore.throats.copy():
            connected_pore = throat.pore1
            if (connected_pore == pore): connected_pore = throat.pore2
            connected_pore.throats.remove(throat)
            pore.throats.remove(throat)
            self.throats.remove(throat)
        # remove pore
        self.pores.remove(pore)
    def connect_pores(self, pore1: Pore, pore2: Pore, r: float, \
        throat_id: int = -1, label: str = LABELS[0]) -> Throat:
        """Establish a throat connection between two pores."""
        throat = Throat(pore1, pore2, r, label, throat_id)
        self.throats.add(throat)
        pore1.throats.add(throat); pore2.throats.add(throat)
        return throat

105
class CellList:
sfritschi's avatar
sfritschi committed
106
107
108
109
    """CellList class dividing physical domain into equally-sized cells
    
    find_candidates() assumes pore is located in interior of domain (non-periodic)
    """
110
    def __init__(self, pores: List[Pore], domainSize: List[float], cellSize: float, \
sfritschi's avatar
sfritschi committed
111
112
                        basenet: Network, throatL: List[float], throatR: List[float], \
                        nInterior: int):
sfritschi's avatar
sfritschi committed
113
        from math import ceil
114
        self.dim = len(domainSize)
115
        self.Lmax = basenet.Lmax
116
        # Number of cells in each dimension
sfritschi's avatar
sfritschi committed
117
        self.nCells = [max(1, ceil(domainSize[i] / cellSize)) for i in range(self.dim)]
118
119
120
121
        # Inverse of cell sizes in domain
        self.invCellSizes = [self.nCells[i] / domainSize[i] for i in range(self.dim)]
        # Total number of cells
        self.totalCells = prod(self.nCells)
122
123
        # Offsets in throat arrays
        self.throatOffsets = [0]
sfritschi's avatar
sfritschi committed
124
125
        # Sort pores according to cell membership
        self.poresSorted = [set() for _ in range(self.totalCells)]
126
127
        # For each pore keep track of nbor pores connected to it
        self.connPores = [set() for _ in range(nInterior)]
sfritschi's avatar
sfritschi committed
128
129
130
131
        # List of (target) throat lengths from base network
        self.throatL = throatL
        # List of throat radii from base network
        self.throatR = throatR
sfritschi's avatar
sfritschi committed
132
        
133
134
        throatTotal = 0
        for i, pore in enumerate(pores):
135
            cellIdx = self.pore_to_index(pore)
136
            self.poresSorted[cellIdx].add(pore)
137
            
138
139
140
            if i < nInterior:
                throatTotal += len(pore.throats)
                self.throatOffsets.append(throatTotal)
141
                
sfritschi's avatar
sfritschi committed
142
        self.poreMatchTable = mp.RawArray('i', throatTotal)
143
144
145
146
147
        
    # Helper for converting 3D index triple into 1D (flattened) index
    def flatten(self, i: int, j: int, k: int) -> int:
        return i + self.nCells[0] * (j + self.nCells[1] * k)
    
sfritschi's avatar
sfritschi committed
148
    def pore_to_triplet(self, pore: Pore) -> List[int]:
149
150
        from math import floor
        # Shift pore position to be positive first (buffer layer)
sfritschi's avatar
sfritschi committed
151
        return [floor( (pore.pos[i] + self.Lmax) * self.invCellSizes[i]) \
152
153
154
155
156
157
158
159
160
                        for i in range(self.dim)]
    # Compute index of given pore in cell list based on position
    def pore_to_index(self, pore: Pore) -> int:
        from math import floor
        # Shift pore position to be positive first (buffer layer)
        cellIdx = self.flatten(*[floor( (pore.pos[i] + self.Lmax) * self.invCellSizes[i]) \
                        for i in range(self.dim)])
        return cellIdx
    
sfritschi's avatar
sfritschi committed
161
162
    # Called in parallel. Find all best neighbor pore candidates and
    # store their indices in poreMatchTable
163
    def match_maker(self, pores: List[Pore], tid: int):
164
        n = len(pores)
165
        for i, pore in enumerate(pores):
166
            if (tid == 0):
167
                print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
168
            
169
            candidates = self.find_candidates(pore)
170
            # Find best candidate for each throat of pore
171
            for idx, throatIdx in enumerate(pore.throats):
sfritschi's avatar
sfritschi committed
172
                lt = self.throatL[throatIdx]
173
                
174
175
176
                if (len(candidates) == 0):
                    self.set_throat(pore.index, -1, idx)
                    break  # no possible candidates left
177
                # Find best match among candidates
sfritschi's avatar
sfritschi committed
178
                matchIdx, _ = min(candidates.items(), key=lambda il: abs(il[1] - lt))
179
180
                # Remove match from candidates
                del candidates[matchIdx]
181
                # Put best match index in poreMatchTable
182
                self.set_throat(pore.index, matchIdx, idx)
183
    
sfritschi's avatar
sfritschi committed
184
    # Compute dictionary of all nbor candidates of given pore
185
    def find_candidates(self, pore: Pore) -> Dict[int, float]:    
186
        candidates  = {}
187
        
188
189
        porePos = pore.pos
        neighborhood = self.nbor_indices(pore)
sfritschi's avatar
sfritschi committed
190
        nearestCellIdx = neighborhood[0]
191
        
192
193
        # Set of pore indices already connected to pore
        pConn = self.connPores[pore.index]
sfritschi's avatar
sfritschi committed
194
        # Search current cell
195
        for nbor in self.poresSorted[nearestCellIdx]:
196
197
198
            nborIdx = nbor.index
            # No self connection; Check if already connected to this pore
            if (nborIdx == pore.index or nbor.originIndex in pConn): continue
199
            l = distance(nbor.pos, porePos)
sfritschi's avatar
sfritschi committed
200
            if (l < self.Lmax):
201
                candidates[nborIdx] = l
sfritschi's avatar
sfritschi committed
202
203
        # Search neighboring cells
        for nborCellIdx in neighborhood[1:]:
204
            for nbor in self.poresSorted[nborCellIdx]:
205
206
207
                nborIdx = nbor.index
                # Check if already connected to this pore
                if (nbor.originIndex in pConn): continue
208
                l = distance(nbor.pos, porePos)
sfritschi's avatar
sfritschi committed
209
                if (l < self.Lmax):
210
                    candidates[nborIdx] = l
sfritschi's avatar
sfritschi committed
211
        
212
        return candidates
sfritschi's avatar
sfritschi committed
213
    
sfritschi's avatar
sfritschi committed
214
    # Set appropriate field in poreMatchTable to index of best candidate
215
216
217
    def set_throat(self, poreIdx: int, matchIdx: int, idx: int):
        index = self.throatOffsets[poreIdx] + idx
        # Set matchIdx
sfritschi's avatar
sfritschi committed
218
        self.poreMatchTable[index] = matchIdx
219
    
sfritschi's avatar
sfritschi committed
220
221
222
223
224
225
    # Retrieve indices of best matches for given pore after they have
    # been computed
    def fetch_matches(self, pore: Pore) -> List[int]:
        lb = self.throatOffsets[pore.index]
        ub = lb + len(pore.throats)
        return self.poreMatchTable[lb:ub]
sfritschi's avatar
sfritschi committed
226
    
227
    # Consider all 27 neighboring cells
228
229
    def nbor_indices(self, pore: Pore) -> List[int]:
        base   = self.pore_to_index(pore)
sfritschi's avatar
sfritschi committed
230
231
        widthx = self.nCells[0]
        widthy = widthx * self.nCells[1]
232
233
234
235
236
237
238
239
240
241
242
243
244
        return [base, base - 1, base + 1,
                base - widthx, base + widthx,
                base - widthy, base + widthy,
                base + widthx - 1, base + widthx + 1,
                base - widthx - 1, base - widthx + 1,
                base + widthy - 1, base + widthy + 1,
                base - widthy - 1, base - widthy + 1,
                base + widthx + widthy, base + widthx + widthy - 1,
                base + widthx + widthy + 1, base + widthx - widthy,
                base + widthx - widthy - 1, base + widthx - widthy + 1,
                base - widthx + widthy, base - widthx + widthy - 1,
                base - widthx + widthy + 1, base - widthx - widthy,
                base - widthx - widthy - 1, base - widthx - widthy + 1]
245
246
247
248
    
    def update_connection(self, poreIdx: int, nborIdx: int):
        self.connPores[poreIdx].add(nborIdx)
        self.connPores[nborIdx].add(poreIdx)
sfritschi's avatar
sfritschi committed
249
250
251
252
        
    def nbor_cell_counts(self, neighborhood: List[int]) -> List[int]:
        return [len(self.poresSorted[i]) for i in neighborhood]
        
sfritschi's avatar
sfritschi committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
    def cell_stats(self):
        cellCounts = [len(self.poresSorted[i]) for i in range(self.totalCells)]
        n    = sum(cellCounts)
        max_ = max(cellCounts)
        min_ = min(cellCounts)
        avg  = n / self.totalCells
        
        """
        print("Cell statistics:")
        print("Number of cells (%d, %d, %d)" % (self.nCells[0], self.nCells[1], self.nCells[2]))
        print("Max. pores per cell: %d (%f%%)" % (max_, max_ / n * 100.))
        print("Min. pores per cell: %d (%f%%)" % (min_, min_ / n * 100.))
        print("Avg. pores per cell: %.1f (%f%%)" % (avg, avg / n * 100.))
        """
        return (max_, min_, avg)

sfritschi's avatar
sfritschi committed
269
270
    # Remove pore from sorted pores
    def expel(self, pore: Pore):
sfritschi's avatar
sfritschi committed
271
        poreCellIdx = self.pore_to_index(pore)
272
        self.poresSorted[poreCellIdx].remove(pore)
273
274
275
276
277
    
    def is_contained(self, pore: Pore) -> bool:
        poreCellIdx = self.pore_to_index(pore)
        return pore in self.poresSorted[poreCellIdx]
        
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
def distance(p1: List[float], p2: List[float]):
    """Distance between points p1 and p2."""
    mag = 0.0
    for c1, c2 in zip(p1,p2): mag += (c1-c2)**2
    return mag**0.5


def prod(v):
    """Cumulative product of vector components."""
    from functools import reduce
    return(reduce(lambda a,b: a*b, v))


def random_direction(d: int) -> List[float]:
    """Unity vector with uniformly distributed orientation."""
    from math import pi, sin, cos, acos
    from random import random
    phi = 2*pi * random()
    if (d == 2):
        return([cos(phi), sin(phi)])
    elif (d == 3):
        theta = acos(2*random() - 1)
        return([sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)])
    else:
        raise ValueError("random_direction supports d = 2 or 3")


def throat_ends(throat: Throat, L: List[float]) -> Tuple[List[float]]:
    """Provide the correct relative end points of a possibly periodic throat."""
    x1 = throat.pore1.pos.copy(); x2 = throat.pore2.pos.copy()
    # account for periodic throats
    if len(throat.label) != 0:
        direct = [int(k) for k in throat.label.split()] # periodicity
        x2 = [x2[k] + direct[k] * L[k] for k in range(len(L))]
    return (x1, x2)



def imperial_read(net_file_pref: str) -> Network:
    """Read network data in Imperial College format.
    
    For format specifics see PhD thesis of Taha Sochi from 2007 at Imperial.
    The names of the network files starts with the string net_file_pref.
    A Network object is returned.
    """
    # initialize network
    from os import sep
    network = Network(label = net_file_pref.split(sep)[-1])
    
    # read pore data
    with open(net_file_pref + '_node1.dat', 'r') as file:
        line = file.readline()
        size = [float(n) for n in line.replace('\t',' ').split()][1:]
        d = len(size) # dimensionality
        f = lambda words: (
            int(words[0]), # id
            [float(n) for n in words[1:(d+1)]], # pos
            [int(throat) for throat in words[(d+4+int(words[d+1])):]], # throats
            LABELS[int(words[d+2+int(words[d+1])]) + \
            2*int(words[d+3+int(words[d+1])])]) # label
        pores = [f(line.replace('\t',' ').split()) for line in file]
    pores.sort(key=lambda pore: pore[0])
    pores = [Pore(pore[1], 0.0, pore[3]) for pore in pores]
    # read pore radii
    with open(net_file_pref + '_node2.dat', 'r') as file:
        f = lambda words: (int(words[0]), float(words[2]))
        radii = [f(line.replace('\t',' ').split()) for line in file]
    radii.sort(key=lambda pore: pore[0])
    for k, pore in enumerate(pores): pore.r = radii[k][1]
    # add pores to network
    for pore in pores: network.add_pore(pore)

    # read throat data
    with open(net_file_pref + '_link1.dat', 'r') as file:
        n = int(file.readline().strip('\n ')) # number of throats
        f = lambda words: (
            int(words[0]), int(words[1])-1, int(words[2])-1, float(words[3]))
        throats = [f(line.replace('\t',' ').split()) for line in file]
    throats.sort(key=lambda throat: throat[0])
    # connect pores
    for throat in throats:
        if (throat[1] >= 0) and (throat[2] >= 0): # no in-/outflow throats
            network.connect_pores(pores[throat[1]], pores[throat[2]], throat[3])

    # determine maximum throat length
    Lmax = 0.0
    for throat in network.throats:
        L = distance(throat.pore1.pos, throat.pore2.pos)
        Lmax = max(L,Lmax)
    network.Lmax = Lmax

    network.lb = [0 for i in range(len(size))]
    network.ub = size
    return network


def plot_network(network: Network, labels: bool = False):
    """Plot pore network.
    
    Returns a figure object. If labels == True, pore and throat ids are shown.
    Pores with label == '' are plotted in black. Pores with labels starting
    with 'in', 'out', 'cut' are plotted in red, blue, magenta, respectively.
    """
    # setup plot
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    plt.ion()
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    # plot throats
    for throat in network.throats:
        if throat.label == LABELS[0]: col = 'black'
        else: col = 'grey'
        p1 = throat.pore1.pos; p2 = throat.pore2.pos
        ax.plot([p1[0], p2[0]], [p1[1], p2[1]], [p1[2], p2[2]],
            linewidth=0.5, color=col)
        if labels:
            ax.text3D((p1[0]+p2[0])/2, (p1[1]+p2[1])/2, (p1[2]+p2[2])/2,
                str(throat.id), fontsize=5)

    # plotting function for pores
    def plot_pores(label, size, color):
        pnts = [pore.pos for pore in network.pores \
            if (pore.label[:max(1,len(label))] == label)]
        if len(pnts) == 0: return # in case of empty point list
        pnts = list(map(list, zip(*pnts))) # transpose double list
        ax.scatter(*pnts, s=size, marker='.', c=color)
        if labels:
            for pore in network.pores:
                if (pore.label[:max(1,len(label))] == label):
                    ax.text3D(*pore.pos, str(pore.id), fontsize=7, color=color)
    # plot ordinary pores
    plot_pores(LABELS[0], 4, 'black')
    # plot inflow pores
    plot_pores(LABELS[1], 7, 'red')
    # plot outflow pores
    plot_pores(LABELS[2], 7, 'blue')
    # plot cut pores
    plot_pores(LABELS[3], 7, 'magenta')

    #ax.axis('equal')
    return fig


def generate_simple_net(n_pores: int, targetsize: List[float], \
    r_pore: float, r_throat: float, coordinatnumb: int, Lmax: float, \
    sd: int = None) -> Network:
    """Generate a spatially periodic network with uniform pore distribution."""
    d = len(targetsize) # number of spatial dimensions
    density = n_pores / prod(targetsize)
    origin = [0.0 for k in range(d)]
    Lb = (1/density)**(1/d)
    basenet = Network(origin, [Lb for k in range(d)], \
        [Pore(origin, r_pore, throats = {k for k in range(coordinatnumb)})], \
        [Throat(None, None, r_throat)], Lmax)
    basenet.label = 'simplenet_' + basenet.label
    return generate_imperial(basenet, targetsize, sd, False)


def generate_imperial(basenet: Network, targetsize: List[float], \
    sd: int = None, correlated: bool = True) -> Network:
    """Generate a new spatially periodic network.

    Based on an existing network basenet, generate a new network of size
    targetsize by using the algorithm outlined on p.54 of the PhD thesis of
    Nasiru Abiodun Idowu from 2009 at Imperial College.
    Pores are uniformly distribute with pore number density as in basenet.
    Pore radii and number of throat connections are sampled from basenet.
    Closest pores are connected. Throat radii are taken from basenet such
    that throat radius and radius-sum of the connected pores are correlated.
    Periodic throats are marked with a periodicity label, e.g., '1 0 -1' in 3d.
    """
    from random import seed, random, randint
    seed(sd)
    d = len(targetsize) # number of spatial dimensions
    # check target size (must be > basenet.Lmax)
    if any([L < basenet.Lmax for L in targetsize]):
        raise NameError('targetsize must be > Lmax of basenet!')

    # pores, uniformly distributed
    # number of pores in new network
    import math
    basesize = list(map(lambda a,b: a-b, basenet.ub, basenet.lb))
    n = math.ceil(prod(targetsize) * len(basenet.pores) / prod(basesize))
    # distribute pores uniformly
    basepores = [pore for pore in basenet.pores] # make indexable
    pores = [Pore(pos=[random()*L for L in targetsize], # pos
        # radius, label, number of throat connections
        r=basepores[idx].r, label=LABELS[0],
        throats=len(basepores[idx].throats)) for j, idx in \
        enumerate([randint(0,len(basepores)-1) for k in range(n)])]

    # add pore buffer layers for spatial periodicity
    n = len(pores)
    _ = __add_buffer_layers(pores, targetsize, basenet.Lmax)

    # throats, connect pores
    # triangulation
    import numpy as np
    from scipy.spatial import Delaunay
    points = np.array([pore.pos for pore in pores])
    tri = Delaunay(points) # possibility for missing (coplanar) points!
    # return neighbors of point pnt based on triangulation tri
    def neighbors_of(pnt, tri):
        return tri.vertex_neighbor_vertices[1][ \
            tri.vertex_neighbor_vertices[0][pnt]: \
            tri.vertex_neighbor_vertices[0][pnt+1]]
    # connect pores
    throats = []
    for k, pore in enumerate(pores[:n]):
        # initialize neighbor lists
        nbors = neighbors_of(k, tri)
        nbors = dict(zip(nbors,
            [distance(pore.pos, pores[j].pos) for j in nbors]))
        oldies = {k:0.0}
        # establish throat connections between pore and its neighbors
        while (pore.throats > 0) \
            and (len(nbors) > 0) and (min(nbors.values()) <= basenet.Lmax):
            # find closest neighbor
            for nbor in nbors:
                if nbors[nbor] == min(nbors.values()): break
            # if nbor is periodic copy, get original pore and periodicity label
            if (pores[nbor].label != LABELS[0]):
                nbor_o, lbl = pores[nbor].label.split(' ',1)
                nbor_o = int(nbor_o)
            else:
                nbor_o = nbor; lbl = LABELS[0]
            # is nbor a valid pore to connect to?
            if not ((nbor_o <= k) or (nbor in oldies) or \
                (pores[nbor_o].throats == 0)):
                # yes -> add connection to pore nbor_o
                throats.append([k, nbor_o, lbl])
                pore.throats = pore.throats - 1
                pores[nbor_o].throats = pores[nbor_o].throats - 1
            # nbor has been dealt with
            oldies[nbor] = nbors[nbor]
            # expand search neighborhood
            for j in neighbors_of(nbor, tri):
                nbors[j] = distance(pores[k].pos, pores[j].pos)
            for oldie in oldies.keys():
                if oldie in nbors: del nbors[oldie]
    
    # initialize throat sets of pores
    k = 0 # count of unrealised throats
    for pore in pores[:n]:
        k = k + pore.throats
        pore.throats = set() # initialize throats set

    # set throat radii
    # sort branch weights (branch weight is prop. to sum of pore radii)
    bw = [pores[throat[0]].r + pores[throat[1]].r for throat in throats]
    bw = list(enumerate(bw))
    bw.sort(key=lambda w: w[1])
    # sample throat radii from basenet
    basethroats = [throat for throat in basenet.throats] # make indexable
    r = [randint(0,len(basethroats)-1) for k in range(len(bw))]
    r = [basethroats[k].r for k in r]
    if correlated: r = sorted(r)
    # assign radii
    for k in range(len(r)): throats[bw[k][0]].append(r[k])

    # assemble and return network
    network = Network(lb=[0.0 for k in range(d)],
        ub=targetsize, Lmax=basenet.Lmax, label='from_' + basenet.label)
    for pore in pores[:n]: network.add_pore(pore)
    for throat in throats: network.connect_pores(pore1=pores[throat[0]],
        pore2=pores[throat[1]], label=throat[2], r=throat[3])
    return network
sfritschi's avatar
sfritschi committed
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638

def cell_list_scaling(basenet: Network, targetsize: List[int], \
    cutoff: float = float('inf'), sd: int = None, mute: bool = False):
    from random import seed, random, randint
    from itertools import product
    seed(sd)
    d = len(targetsize) # number of spatial dimensions
    # check target size multiplicator (must be >= 1)
    if any([i < 1 for i in targetsize]):
        raise NameError('targetsize must be >= 1!')
    # size of new network
    L = list(map(lambda lb,ub,i: (ub-lb)*i, basenet.lb, basenet.ub, targetsize))
    # check target size (must be > basenet.Lmax)
    if any([Lc < basenet.Lmax for Lc in L]):
        raise NameError('targetsize leads network < Lmax of basenet!')

    # pores, distributed based on dendrogram of basenet
    if (not mute): print("distributing pores...", end="", flush=True)
    # make indexable and discard in-/outflow pores
    basepores = [pore for pore in basenet.pores \
        if ((pore.label[:len(LABELS[1])] != LABELS[1]) \
        and (pore.label[:len(LABELS[2])] != LABELS[2]))]
    # dendrogram-based or uniform pore distribution
    pores = []
    if (cutoff != cutoff): # uniform
        print("\b"*21 + "uniform pore distribution")
        # loop over network sections & add uniformly distributed pores drawn from basenet
        for s in product(*[range(k) for k in targetsize]):
            for k, pore in enumerate(basepores):
                pos = [random()*l for l in L]
                pores.append(Pore(pos=pos, r=pore.r, label=LABELS[0], 
                    throats=pore.throats.copy()))
    else: # dendrogram-based
        # extract cluster hierarchy from basenet
        centroids = [pore.pos for pore in basepores]
        from scipy.cluster import hierarchy
        clustree = hierarchy.linkage(centroids, method = 'centroid')
        # determine cluster centroids and weights
        weights = len(basepores)*[1] # pores have weight = 1
        for cluster in clustree:
            p1 = int(cluster[0]); p2 = int(cluster[1])
            w1 = weights[p1]; w2 = weights[p2]
            centroids.append(list(map(lambda x1,x2: (w1*x1 + w2*x2)/(w1 + w2),
                centroids[p1], centroids[p2])))
            weights.append(w1 + w2)
        # loop over network sections (twisted basenet copies)
        for s in product(*[range(k) for k in targetsize]):
            # twist centers of sufficiently small cluster
            touched = [False]*len(centroids)
            # positions of pores based on rotations of linked cluster-pairs
            for k in range(len(centroids)-1,len(basepores)-1,-1):
                i = k-len(basepores) # index of cluster in clustree
                ctr = centroids[k] # center of rotation
                dist = clustree[i,2] # cluster distance
                # check if cluster/pore needs to be twisted
                if ((dist < cutoff) or touched[k]):
                    p1 = int(clustree[i,0]); p2 = int(clustree[i,1])
                    w1 = weights[p1]; w2 = weights[p2]
                    vec = random_direction(d)
                    centroids[p1] = \
                        [ctr[j] + w2*dist/(w1+w2)*vec[j] for j in range(d)]
                    centroids[p2] = \
                        [ctr[j] - w1*dist/(w1+w2)*vec[j] for j in range(d)]
                    touched[p1] = touched[p2] = True
            # add section pores in random order to pores of new network
            for k, pore in enumerate(basepores):
                pos = [c-lb + si*(ub-lb) \
                    for c,si,lb,ub in zip(centroids[k],s,basenet.lb,basenet.ub)]
                pores.insert(randint(0,len(pores)), Pore(pos=pos,
                    r=pore.r, label=LABELS[0], throats=pore.throats.copy()))
        print("\b"*21 + "left {0:d} of {1:d} clusters (incl. {2:d} pores) untouched".\
            format(sum([int(not j) for j in touched]), len(touched), len(basepores)))
        # flip pores outside back into domain
        for pore in pores:
            for k in range(d):
                if pore.pos[k] < 0:
                    pore.pos[k] = L[k] + pore.pos[k]
                elif pore.pos[k] >= L[k]:
                    pore.pos[k] = pore.pos[k] - L[k]

    # add pore buffer layers for spatial periodicity
    n = len(pores)
    copies = __add_buffer_layers(pores, L, basenet.Lmax)
    
    # Domain size including periodic buffer layers on all sides
    trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
    # Initialize cell-list; Place each pore in respective cell-set
    cellList = CellList(pores, trueDomainSize, basenet.Lmax, basenet.Lmax)
    # Print cell stats to console
    max_, min_, avg = cellList.cell_stats()
    return (max_, min_, avg)
    
639
    
640
def generate_dendrogram(basenet: Network, targetsize: List[int], \
641
    cutoff: float = float('inf'), sd: int = None, nthreads: int = mp.cpu_count(), \
642
    mute: bool = False) -> Network:
643
644
645
646
647
648
649
650
651
652
653
    """Generate a new spatially periodic network.

    Based on an existing network basenet, generate a new network of size
    targetsize by using a dendrogram-based algorithm that accounts for the
    spatial distribution/clustering of pores.
    The pore-clustering dendrogram is taken from basenet. Pore radii and number
    of throats are sampled from basenet. Pores are connected based on the
    throat lengths from the basenet data. Throat radii are taken from basenet.
    Periodic throats are marked with a periodicity label, e.g., '1 0 -1' in 3d.
    If cutoff == float('nan'), pores are uniformly distributed.
    """
sfritschi's avatar
sfritschi committed
654
    
sfritschi's avatar
sfritschi committed
655
656
    import multiprocessing as mp
    from random import seed, random, randint
657
    from itertools import product, accumulate
658
659
660
661
662
663
664
665
666
667
    seed(sd)
    d = len(targetsize) # number of spatial dimensions
    # check target size multiplicator (must be >= 1)
    if any([i < 1 for i in targetsize]):
        raise NameError('targetsize must be >= 1!')
    # size of new network
    L = list(map(lambda lb,ub,i: (ub-lb)*i, basenet.lb, basenet.ub, targetsize))
    # check target size (must be > basenet.Lmax)
    if any([Lc < basenet.Lmax for Lc in L]):
        raise NameError('targetsize leads network < Lmax of basenet!')
sfritschi's avatar
sfritschi committed
668
669
670
    # check valid number of threads
    if (nthreads < 1):
        raise ValueError('Number of threads must be >= 1!')
sfritschi's avatar
sfritschi committed
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
    
    # Initialize throat lengths and radii of base network to be used
    # later
    throatIdxMap = {}  # {int (throat id): int (throat index)}
    throatBaseIdx = 0
    throatL = []
    throatR = []
    baseL = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
    for pore in basenet.pores:
        for throat in pore.throats:
            lt = distance(*throat_ends(throat, baseL))
            throatL.append(lt)
            throatR.append(throat.r)
            throatIdxMap[throat] = throatBaseIdx
            throatBaseIdx += 1
sfritschi's avatar
sfritschi committed
686
        
sfritschi's avatar
sfritschi committed
687
    del baseL  # no longer needed
688
    # pores, distributed based on dendrogram of basenet
sfritschi's avatar
sfritschi committed
689
    if (not mute): print("distributing pores...")
690
691
692
693
694
695
696
    # make indexable and discard in-/outflow pores
    basepores = [pore for pore in basenet.pores \
        if ((pore.label[:len(LABELS[1])] != LABELS[1]) \
        and (pore.label[:len(LABELS[2])] != LABELS[2]))]
    # dendrogram-based or uniform pore distribution
    pores = []
    if (cutoff != cutoff): # uniform
697
        runningIndex = 0
698
        print("uniform pore distribution")
699
700
701
702
        # loop over network sections & add uniformly distributed pores drawn from basenet
        for s in product(*[range(k) for k in targetsize]):
            for k, pore in enumerate(basepores):
                pos = [random()*l for l in L]
sfritschi's avatar
sfritschi committed
703
                poreThroatIdx = set(throatIdxMap[throat] for throat in pore.throats)
704
                pores.append(Pore(pos=pos, r=pore.r, label=LABELS[0], 
sfritschi's avatar
sfritschi committed
705
                    throats=poreThroatIdx, index=runningIndex))
sfritschi's avatar
sfritschi committed
706
                runningIndex += 1
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
    else: # dendrogram-based
        # extract cluster hierarchy from basenet
        centroids = [pore.pos for pore in basepores]
        from scipy.cluster import hierarchy
        clustree = hierarchy.linkage(centroids, method = 'centroid')
        # determine cluster centroids and weights
        weights = len(basepores)*[1] # pores have weight = 1
        for cluster in clustree:
            p1 = int(cluster[0]); p2 = int(cluster[1])
            w1 = weights[p1]; w2 = weights[p2]
            centroids.append(list(map(lambda x1,x2: (w1*x1 + w2*x2)/(w1 + w2),
                centroids[p1], centroids[p2])))
            weights.append(w1 + w2)
        # loop over network sections (twisted basenet copies)
        for s in product(*[range(k) for k in targetsize]):
            # twist centers of sufficiently small cluster
            touched = [False]*len(centroids)
            # positions of pores based on rotations of linked cluster-pairs
            for k in range(len(centroids)-1,len(basepores)-1,-1):
                i = k-len(basepores) # index of cluster in clustree
                ctr = centroids[k] # center of rotation
                dist = clustree[i,2] # cluster distance
                # check if cluster/pore needs to be twisted
                if ((dist < cutoff) or touched[k]):
                    p1 = int(clustree[i,0]); p2 = int(clustree[i,1])
                    w1 = weights[p1]; w2 = weights[p2]
                    vec = random_direction(d)
                    centroids[p1] = \
                        [ctr[j] + w2*dist/(w1+w2)*vec[j] for j in range(d)]
                    centroids[p2] = \
                        [ctr[j] - w1*dist/(w1+w2)*vec[j] for j in range(d)]
                    touched[p1] = touched[p2] = True
            # add section pores in random order to pores of new network
sfritschi's avatar
sfritschi committed
740
            baseThroatIdx = 0
741
            for k, pore in enumerate(basepores):
sfritschi's avatar
sfritschi committed
742
                [throat in pore.throats]
743
744
                pos = [c-lb + si*(ub-lb) \
                    for c,si,lb,ub in zip(centroids[k],s,basenet.lb,basenet.ub)]
sfritschi's avatar
sfritschi committed
745
                poreThroatIdx = set(throatIdxMap[throat] for throat in pore.throats)
sfritschi's avatar
sfritschi committed
746
                pores.insert(randint(0,len(pores)), Pore(pos=pos,
sfritschi's avatar
sfritschi committed
747
748
749
750
751
                    r=pore.r, label=LABELS[0], throats=poreThroatIdx))
        
        # Throat index map no longer needed
        del throatIdxMap
        
752
        print("left {0:d} of {1:d} clusters (incl. {2:d} pores) untouched".\
753
            format(sum([int(not j) for j in touched]), len(touched), len(basepores)))
sfritschi's avatar
sfritschi committed
754
755
756
        # flip pores outside back into domain and set respective index of all pores
        for i, pore in enumerate(pores):
            pore.index = i
757
            pore.originIndex = i  # pore is already the original
758
759
760
761
762
763
764
765
766
            for k in range(d):
                if pore.pos[k] < 0:
                    pore.pos[k] = L[k] + pore.pos[k]
                elif pore.pos[k] >= L[k]:
                    pore.pos[k] = pore.pos[k] - L[k]

    # add pore buffer layers for spatial periodicity
    n = len(pores)
    copies = __add_buffer_layers(pores, L, basenet.Lmax)
767
            
768
769
    # Domain size including periodic buffer layers on all sides
    trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
sfritschi's avatar
sfritschi committed
770
    # Initialize cell-list; Place each pore in respective cell-set
sfritschi's avatar
sfritschi committed
771
772
    cellList = CellList(pores, trueDomainSize, basenet.Lmax, basenet, \
                          throatL, throatR, n)
sfritschi's avatar
sfritschi committed
773
    
774
775
776
    # Divide by two (since we are counting every throat twice)
    totalThroats = len(cellList.poreMatchTable) // 2
    
777
    # throats, connect pores
sfritschi's avatar
sfritschi committed
778
    if (not mute): print("\nconnecting")
779
780
    throats = []
    
781
    # DEBUG
782
783
784
785
786
787
788
    """
    for i, pore in enumerate(pores):
        if (pore.label != LABELS[0]):
            j, _ = pore.label.split(' ',1)
            assert(pore.originIndex == int(j))
        assert(pore.index == i)
    """
789
790
    # coordination number for each pore
    #coords = [len(pore.throats) for pore in pores[:n]]
791
792
    avg_throat_diff = 0.
    count = 0
793
    
sfritschi's avatar
sfritschi committed
794
    # Throats-To-Be-Realized
sfritschi's avatar
sfritschi committed
795
796
    # Copy throats of pores to keep track of the ones already realized
    ttbrs = [pore.throats.copy() for pore in pores[:n]]
sfritschi's avatar
sfritschi committed
797
    
798
799
800
    if (not mute): print("Total throats: %d" % totalThroats)
    
    throatsLeft = totalThroats
801
    throatsUnrealized = totalThroats
802
    
803
804
805
    # Threshold of remaining pores at which we compute the best
    # matches serially instead of using threads
    serialThresh = nthreads * 4
806
    # Pores in interior that are not fully-connected yet
807
    poresRemain = pores[:n]
sfritschi's avatar
sfritschi committed
808
    while (throatsUnrealized > 0):
809
        
sfritschi's avatar
sfritschi committed
810
        if (not mute):
811
            print("Throats left: %d Throats Unrealized: %d (%.2f%%)" % \
sfritschi's avatar
sfritschi committed
812
                 (throatsLeft, throatsUnrealized, throatsLeft / totalThroats * 100.))
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
        
        nRemain = len(poresRemain)
        if (nRemain > serialThresh):
            # Compute loads and displacements for all threads
            remainder = nRemain % nthreads
            loads = [nRemain // nthreads + (i < remainder) for i in range(nthreads)]
        
            # Compute inclusive-scan yielding offsets in poresRemain threads
            displs = [0]
            displs += accumulate(loads)
            
            workers = []
            if (not mute): print("computing best matches")
            # Compute best matches for each pore in parallel using available threads
            for tid in range(nthreads):
                poreList = poresRemain[displs[tid] : displs[tid+1]]
                worker = mp.Process(target=cellList.match_maker, args=(poreList, tid))
                worker.start()
                workers.append(worker)
            
            for worker in workers:
                worker.join()
        else:  # serial
            if (not mute): print("\nSearching serially...")
            cellList.match_maker(poresRemain, 0)
sfritschi's avatar
sfritschi committed
838
        
839
        # Iterate over remaining (reduced) pores
840
        for it, pore in enumerate(poresRemain):
841
            
842
            if (not mute):
843
                print(f"progress {((it+1) / n)*100.:.1f}%", end="\r", flush=True)
844
            # Set of throats to be realized
sfritschi's avatar
sfritschi committed
845
            poreTtbr = ttbrs[pore.index]
846
            if (len(poreTtbr) == 0): continue  # Nothing left to do
847
848
            # Set of indices of connected pores
            pConn = cellList.connPores[pore.index]
sfritschi's avatar
sfritschi committed
849
            matches = cellList.fetch_matches(pore)
850
            
851
            for throatIdx, match in zip(pore.throats, matches):
852
                
853
854
                if throatIdx not in poreTtbr:
                    continue  # Throat already realized
855
                
856
857
858
                if match == -1:  # give up
                    poreTtbr.clear()
                    break
859
                    
860
861
862
863
864
865
                nbor = pores[match]
                nborc = nbor # memorize potential periodic copy
                # find original of buffer layer pore
                if (nbor.label != LABELS[0]):
                    j, lbl = nbor.label.split(' ',1)
                    nbor = pores[int(j)]
866
                    assert(nborc.originIndex == int(j))
867
868
869
870
871
                else:
                    lbl = LABELS[0]
                
                nborTtbr = ttbrs[nbor.index]
                
872
873
                # Check that nbor is not fully connected and that
                # this pore is not already connected to it
874
                if (len(nborTtbr) == 0 or nbor.index in pConn):
875
                    continue
876
877
                # Update sets of connected pores
                cellList.update_connection(pore.index, nbor.index)
878
                
sfritschi's avatar
sfritschi committed
879
880
                r = throatR[throatIdx]
                targetLen = throatL[throatIdx]
881
882
883
884
885
886
                # Compute actual throat length between pores
                lt = distance(pore.pos, nborc.pos)
                
                # connect pore to nbor
                throats.append([pore, nbor, lbl, r])
                # remove most similar throat of nbor
sfritschi's avatar
sfritschi committed
887
                ranking = [(nborThroatIdx, abs(throatL[nborThroatIdx] - lt)) \
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
                            for nborThroatIdx in nborTtbr]
                nborThroatIdx, _ = min(ranking, key=lambda il: il[1])
                
                diff = abs(targetLen - lt)
                avg_throat_diff += diff
                count += 1
                
                # Successfully connected this throat
                poreTtbr.remove(throatIdx)
                # Remove throat from neighbor
                nborTtbr.remove(nborThroatIdx)
                
                # remove fully-connected nbor from search neighborhood
                if (len(nborTtbr) == 0):
                    cellList.expel(nbor)
                    
                    for cpore in copies[nbor]:
                        cellList.expel(cpore)
906
            
907
908
909
            # remove fully-connected pore from search neighborhood
            if (len(poreTtbr) == 0):
                cellList.expel(pore)
910
            
911
                for cpore in copies[pore]:
sfritschi's avatar
sfritschi committed
912
                    cellList.expel(cpore)
913
        
914
915
916
917
918
919
920
921
922
        # DEBUG
        """
        for pore in poresRemain:
            matches = cellList.fetch_matches(pore)
            for match in matches:
                if (match == -1): break  # No match found
                mPore = pores[pores[match].originIndex]  # original pore
                assert(mPore in poresRemain)
        """
923
        # Copy throats back to original pores and update #throats left
924
        # and compute remaining (reduced) pores
925
        throatCount = 0
926
927
        nextPores = []
        for pore in poresRemain:
928
            poreTtbr = ttbrs[pore.index]
929
            pore.throats = poreTtbr.copy()
930
931
            
            if len(poreTtbr) > 0:
932
                throatCount += len(poreTtbr)
933
                # Insert in random order for improved load balancing
934
                nextPores.insert(randint(0,len(nextPores)), pore)
935
936
937
938
            else:
                # Free memory associated with conn-pores of fully
                # connected pores
                cellList.connPores[pore.index].clear()
939
        
940
        poresRemain = nextPores
sfritschi's avatar
sfritschi committed
941
942
        throatsUnrealized = throatCount // 2  # Counted all throats twice
        throatsLeft = totalThroats - len(throats)
943
        
944
    
945
    # Unrealized throats
sfritschi's avatar
sfritschi committed
946
    nUnrealized = throatsLeft
947
    
948
    # Free memory (not needed anymore)
949
    del cellList
950
    del ttbrs
951
    
sfritschi's avatar
sfritschi committed
952
    # DEBUG
953
    percent = nUnrealized  / totalThroats * 100.
954
    print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.3f}%)".\
955
        format(len(throats), nUnrealized, percent))
956
957
958
    avg_throat_diff /= count
    print("Avg. throat length difference: %e" % avg_throat_diff)
    print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))
959

960
961
962
    # assemble and return network
    network = Network(lb=[0.0 for k in range(d)],
        ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label)
963
    for pore in pores[:n]:
964
        pore.throats.clear()  # remove (potentially) remaining pores
965
        network.add_pore(pore)
966
967
    for throat in throats: network.connect_pores(pore1=throat[0],
        pore2=throat[1], label=throat[2], r=throat[3])
968
    
969
    # DEBUG
970
    """
971
972
973
974
975
976
977
978
979
980
981
982
    # Check for correctness of generated network
    checkUnrealized = 0
    for pore in network.pores:
        nthroats = len(pore.throats)
        # Make sure pore coordination number is respected
        coordNum = coords[pore.index]
        assert(nthroats <= coordNum)
        checkUnrealized += coordNum - nthroats
        
        nborInd = set()
        for throat in pore.throats:
            assert(pore.index == throat.pore1.index or pore.index == throat.pore2.index)
983
984
985
986
            nbor = throat.pore2 if (throat.pore1.index == pore.index) else throat.pore1
            if (nbor.label != LABELS[0]):
                nbor = pores[int(j)]
            nborInd.add(nbor.index)
987
988
989
990
991
        # Make sure pore isn't connected to the same nbor twice
        assert(len(nborInd) == nthroats)
    
    checkUnrealized //= 2  # counted every throat twice
    # Make sure missing throats are accounted for
992
    assert(checkUnrealized == nUnrealized)
993
    """
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
    return network


def __add_buffer_layers(pores: List[Pore], targetsize: List[float],
    Lbuffer: float) -> Dict[Pore,Set[Pore]]:
    """Add pore buffer layers for spatial periodicity.
    
    Returns a dict where each key is a pore and the corresponding value
    is a set containing the copies of that pore.
    """
    n = len(pores) # number of pores inside domain
    d = len(targetsize) # number of spatial dimensions
    # label = pore index + periodicity label (-1,0,1) in each dim.
    for k, pore in enumerate(pores): pore.label = str(k) + d * ' 0'
    # add pore buffer layers
    copies = {pore:set() for pore in pores} # {pore:{set of periodic copies}}
sfritschi's avatar
sfritschi committed
1010
    runningIndex = n
1011
1012
1013
1014
1015
    for k in range(d):
        # left bound
        pores_layer = [pore for pore in pores \
            if targetsize[k]-Lbuffer <= pore.pos[k] < targetsize[k]]
        for pore in pores_layer:
1016
1017
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, \
                            id=pore.id, index=runningIndex)
1018
1019
1020
            pcopy.pos[k] = pcopy.pos[k]-targetsize[k]
            lbl = pore.label.split(); lbl[k + 1] = '-1'
            pcopy.label = ' '.join(lbl) # periodicity label
1021
1022
            originIndex = int(lbl[0])
            pcopy.originIndex = originIndex
1023
            pores.append(pcopy)
1024
            copies[pores[originIndex]].add(pcopy)
sfritschi's avatar
sfritschi committed
1025
            runningIndex += 1
1026
1027
1028
1029
        # right bound
        pores_layer = [pore for pore in pores \
            if 0 <= pore.pos[k] < Lbuffer]
        for pore in pores_layer:
1030
1031
            pcopy = Pore(pore.pos, pore.r, throats=pore.throats, \
                            id=pore.id, index=runningIndex)
1032
1033
1034
            pcopy.pos[k] = pcopy.pos[k]+targetsize[k]
            lbl = pore.label.split(); lbl[k + 1] = '1'
            pcopy.label = ' '.join(lbl) # periodicity label
1035
1036
            originIndex = int(lbl[0])
            pcopy.originIndex = originIndex
1037
            pores.append(pcopy)
1038
            copies[pores[originIndex]].add(pcopy)
sfritschi's avatar
sfritschi committed
1039
            runningIndex += 1
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
    # reset labels of pores inside domain (without buffer layers)
    for pore in pores[:n]: pore.label = LABELS[0]
    return copies


def open_periodic_network(network: Network, c: int, inouts: Dict = {}):
    """Open periodic throats in the spatial direction c>0 of a network.
    
    To avoid identical in/out pores resulting from different periodic throats
    connecting to the same original pore, one can keep track of newly generated
    in/out pores through the inouts pore dict {orig1:{"0 0 1":copy1, ...}, ...}.
    """
    # loop over (periodic) throats
    for throat in network.throats.copy(): # shallow copy (iteration changes set)
        # is throat a periodic connection in direction c?
        if len(throat.label) == 0: continue
        direct = [int(c) for c in throat.label.split()]
        if direct[c-1] == 0: continue
        # open periodic throat connection
        __open_periodic_throat(network, throat, direct, c, inouts)


def __open_periodic_throat(network: Network,
    throat: Throat, direct: List[int], c: int, inouts: Set[Pore]) \
    -> Set[Throat]:
    """Open a periodic throat connection.
    
    direct is the encoded periodicity label of the throat and c>0 is the
    spatial direction used for the in-/outflow labeling. The inouts dict
    is used to avoid the formation of pores with identical positions.
    """
    d = len(network.lb) # number of spatial dimensions
    p1 = throat.pore1; p2 = throat.pore2; r = throat.r
    # identify existing in/out pore or construct new one
    # needed in case of many periodic throats connecting to the same pore
    def get_inout_pore(pore_o, direct, pos, label, inouts) -> Pore:
        dirkey = str([abs(d) for d in direct]) # periodicity direction key
        # find existing inout pore copy (or periodic original)
        if pore_o in inouts: # is there an inout copy of the original pore?
            if dirkey in inouts[pore_o]: # is there copy in right direct?
                return inouts[pore_o][dirkey] # yes, return inout copy
        # if no inout copy exists, create a new one
        pore = Pore(pos, pore_o.r, label)
        if pore_o in inouts: # add new direct
            inout = inouts[pore_o] # dict dirkey:inout-pore
            inout[dirkey] = pore
        else: inouts[pore_o] = {dirkey:pore} # add original to inout dict
        return pore
    # pore2
    pos = [0.0 for j in range(d)] # init new position object
    # copy and shift pore2
    for j in range(d):
        pos[j] = p2.pos[j] + direct[j] * (network.ub[j]-network.lb[j])
    p2c = get_inout_pore(p2, direct, pos,
        LABELS[1+round((direct[c-1]+1)/2)] + str(c), inouts)
    network.add_pore(p2c)
    t1 = network.connect_pores(p1, p2c, throat.r)
    # pore1
    pos = [0.0 for j in range(d)] # init new position object
    # copy and shift pore1
    for j in range(d):
        pos[j] = p1.pos[j] - direct[j] * (network.ub[j]-network.lb[j])
    p1c = get_inout_pore(p1, direct, pos,
        LABELS[1+round((-direct[c-1]+1)/2)] + str(c), inouts)
    network.add_pore(p1c)
    t2 = network.connect_pores(p2, p1c, throat.r)
    # remove periodic throat
    network.throats.remove(throat)
    p1.throats.remove(throat); p2.throats.remove(throat)
    return {t1, t2}


def cut_network(network: Network, x: float, c: int, label: str = '',
    inouts: Dict = {}):
    """Cut a (periodic) network at position x_c with c>0.

    New pores are labelled and introduced at throat intersection points
    with the cutting plane. Before cutting, periodic throat connections
    in the c-direction are opened by calling open_periodic_network,
    which makes use of the inouts pore dictionary.
    """
    d = len(network.lb) # number of spatial dimensions
    # open periodic pores in direction c
    open_periodic_network(network, c, inouts)
    # cut throats incl. cut periodic ones being periodic normal to direction c
    throats = network.throats.copy() # shallow copy (iteration changes set)
    while len(throats) > 0:
        throats_next = set() # cut possible periodic pores in a 2nd round
        # cut throats ...
        for throat in throats:
            p1 = throat.pore1; p2 = throat.pore2
            # ... intersecting with cutting plane at x_c
            if (p1.pos[c-1] - x) * (x - p2.pos[c-1]) > 0:
                # treat periodic pores
                if len(throat.label) != 0:
                    direct = \
                        [int(k) for k in throat.label.split()] # periodicity
                    for k, cn in enumerate(direct):
                        if cn != 0: break # find first +/-1 periodic component
                    throats_periodic = __open_periodic_throat( \
                        network, throat, direct, k+1, inouts)
                    throats_next = throats_next | throats_periodic
                    continue
                # cut position
                f = (x - p1.pos[c-1]) / (p2.pos[c-1] - p1.pos[c-1])
                pos = [(p2.pos[j]-p1.pos[j])*f + p1.pos[j] for j in range(d)]
                pos[c-1] = x # to avoid tiny rounding errors
                # introduce and connect pore
                pore = Pore(pos, r=0.0, label=LABELS[3] + str(c) + label)
                network.add_pore(pore)
                # throat p1 - pore (reconnect throat from p2 to p1)
                throat.pore2 = pore; pore.throats.add(throat)
                # throat pore - p2 (new throat object)
                p2.throats.remove(throat)
                network.connect_pores(pore, p2, throat.r)
        # cut newly opened periodic throats
        throats = throats_next


def erase_network(network: Network, x: float, c: int, direct: bool, label: str):
    """Remove pores and connected throats from a network.
    
    Pores with positions < or > x_c for direct = True or False, respectively,
    are removed (c>0). At the same time, pores at x_c are labelled and network
    bounds updated.
    """
    for pore in network.pores.copy():
        # remove pores (and connected throats)
        if ((pore.pos[c-1] < x) and direct) or \
            ((x < pore.pos[c-1]) and not direct): network.remove_pore(pore)
        # relabel pores on cutting plane x_c
        elif pore.pos[c-1] == x: pore.label = label
    # update network bounds
    if direct: network.lb[c-1] = x
    else: network.ub[c-1] = x


def save_network_to(filename: str, network: Network):
    """Save pore network to hdf5 file."""
    import h5py
    f = h5py.File(filename, 'x')
    # write network label
    label = network.label.encode('ascii','ignore')
    f.create_dataset('network_label', (1,), dtype='S'+str(len(label)),
        data=[label])
    # bounds, Lmax
    f.create_dataset('lb', data=network.lb)
    f.create_dataset('ub', data=network.ub)
    f.create_dataset('Lmax', data=[network.Lmax])

    # write pore data
    p_grp = f.create_group('pores')
    # id
    wrk = [pore.id for pore in network.pores]
    p_grp.create_dataset('id', data=wrk)
    # radius
    wrk = [pore.r for pore in network.pores]
    p_grp.create_dataset('r', data=wrk)
    # position
    for k in range(len(network.lb)):
        wrk = [pore.pos[k] for pore in network.pores]
        p_grp.create_dataset('pos/x' + str(k), data=wrk)
    # label
    wrk = [pore.id for pore in network.pores if len(pore.label) != 0]
    if len(wrk) > 0:
        p_grp.create_dataset('label/id', data=wrk)
        wrk = [pore.label.encode('ascii','ignore') \
            for pore in network.pores if len(pore.label) != 0]
        from functools import reduce
        length = reduce(lambda a,b: max(a,b), [len(lbl) for lbl in wrk])
        p_grp.create_dataset('label/strg', (len(wrk),), dtype='S'+str(length),
            data=wrk)
        
    # write throat data
    t_grp = f.create_group('throats')
    # id
    wrk = [throat.id for throat in network.throats]
    t_grp.create_dataset('id', data=wrk)
    # radius
    wrk = [throat.r for throat in network.throats]
    t_grp.create_dataset('r', data=wrk)
    # pores
    wrk = [throat.pore1.id for throat in network.throats]
    t_grp.create_dataset('pore1', data=wrk)
    wrk = [throat.pore2.id for throat in network.throats]
    t_grp.create_dataset('pore2', data=wrk)
    # label
    wrk = [throat.id for throat in network.throats \
        if len(throat.label) != 0]
    if len(wrk) > 0:
        t_grp.create_dataset('label/id', data=wrk)
        wrk = [throat.label.encode('ascii','ignore') \
            for throat in network.throats if len(throat.label) != 0]
        from functools import reduce
        length = reduce(lambda a,b: max(a,b), [len(lbl) for lbl in wrk])
        t_grp.create_dataset('label/strg', (len(wrk),), dtype='S'+str(length),
            data=wrk)


sfritschi's avatar
sfritschi committed
1239
def load_network_from(filename: str) -> Network:
1240
1241
1242
1243
1244
1245
    """Load pore network from hdf5 file."""
    import h5py
    import numpy as np
    f = h5py.File(filename, 'r')
    # global network properties
    network = Network(label=f['network_label'][0].decode('utf-8'),
sfritschi's avatar
sfritschi committed
1246
        lb=list(f['lb']), ub=list(f['ub']), Lmax=f['Lmax'][0])
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264

    # pores
    pid = np.array(f['pores/id'])
    pores = {} # id:pore dictionary
    for id in pid:
        pore = Pore(pos=[0.0 for k in network.lb], r=0.0, id=id)
        pores[id] = pore
        network.add_pore(pore)
    # radius
    for id, r in zip(pid, np.array(f['pores/r'])): pores[id].r = r
    # position
    for k in range(len(network.lb)):
        for id, x in zip(pid, np.array(f['pores/pos/x' + str(k)])):
            pores[id].pos[k] = x
    # labels
    if 'label' in f['pores']: # are there any labels != ''?
        for id, strg in zip(f['pores/label/id'],f['pores/label/strg']):
            pores[id].label = strg.decode('utf-8')
1265
        
1266
1267
1268
    # throats
    tid = np.array(f['throats/id'])
    throats = {} # id:throat dictionary
sfritschi's avatar
sfritschi committed
1269
    for id, pore1, pore2, r in zip(tid,
1270
        [pores[id] for id in np.array(f['throats/pore1'])],
sfritschi's avatar
sfritschi committed
1271
1272
        [pores[id] for id in np.array(f['throats/pore2'])],
        np.array(f['throats/r'])):
sfritschi's avatar
sfritschi committed
1273
        throats[id] = network.connect_pores(pore1, pore2, r, id)
sfritschi's avatar
sfritschi committed
1274

1275
1276
1277
1278
1279
1280
    # labels
    if 'label' in f['throats']: # are there any labels != ''?
        for id, strg in zip(f['throats/label/id'],f['throats/label/strg']):
            throats[id].label = strg.decode('utf-8')

    return network