example1.py 4.89 KB
Newer Older
1
import frackit.geometry as geometry
Dennis Gläser's avatar
Dennis Gläser committed
2
3
import random
import math
4
5

# we use the unit cube as domain
Dennis Gläser's avatar
Dennis Gläser committed
6
7
box = geometry.Box(xmin=0.0, ymin=0.0, zmin=0.0,
                   xmax=1.0, ymax=1.0, zmax=1.0)
8
9
10
11
12
13

# we sample points uniformly within the domain
from frackit.sampling import makeUniformPointSampler

# returns a sampler from a gaussian distribution with mean and std deviation
def gaussianSampler(mean, stdDev):
Dennis Gläser's avatar
Dennis Gläser committed
14
    def sample(): return random.gauss(mean, stdDev)
15
16
    return sample

Dennis Gläser's avatar
Dennis Gläser committed
17
18
# returns a sampler from a uniform distribution between min and max
def uniformSampler(min, max):
Dennis Gläser's avatar
Dennis Gläser committed
19
    def sample(): return random.uniform(min, max)
Dennis Gläser's avatar
Dennis Gläser committed
20
21
    return sample

Dennis Gläser's avatar
Dennis Gläser committed
22
# we sample quadrilaterals within the box with gaussian distributions for the parameters
23
from frackit.sampling import QuadrilateralSampler as QuadSampler
Dennis Gläser's avatar
Dennis Gläser committed
24
25
26
27
28
quadSampler1 = QuadSampler(pointSampler        = makeUniformPointSampler(box),
                           strikeAngleSampler  = gaussianSampler(math.radians(45.0), math.radians(5.0)),
                           dipAngleSampler     = gaussianSampler(math.radians(90.0), math.radians(5.0)),
                           strikeLengthSampler = uniformSampler(0.4, 0.8),
                           dipLengthSampler    = uniformSampler(0.4, 0.8))
29
30

# sampler for quadrilaterals of the secondary orientation
Dennis Gläser's avatar
Dennis Gläser committed
31
32
33
34
35
quadSampler2 = QuadSampler(pointSampler        = makeUniformPointSampler(box),
                           strikeAngleSampler  = gaussianSampler(math.radians(0.0), math.radians(5.0)),
                           dipAngleSampler     = gaussianSampler(math.radians(0.0), math.radians(5.0)),
                           strikeLengthSampler = uniformSampler(0.4, 0.8),
                           dipLengthSampler    = uniformSampler(0.4, 0.8))
36
37
38
39
40
41
42
43
44
45

# We want to enforce some constraints on the set of quadrilaterals.
# In particular, for entities of the same set we want a minimum spacing
# distance of 5cm, and the quadrilaterals must not intersect in angles
# less than 30°. Moreover, if they intersect, we don't want intersection
# edges whose length is smaller than 5cm, and, the intersection should not
# be too close to the boundary of one of two intersecting quadrilaterals. Here: 5cm.
from frackit.entitynetwork import EntityNetworkConstraints
constraintsOnSelf = EntityNetworkConstraints()
constraintsOnSelf.setMinDistance(0.05)
Dennis Gläser's avatar
Dennis Gläser committed
46
constraintsOnSelf.setMinIntersectingAngle(math.radians(30.0))
47
48
49
50
51
52
constraintsOnSelf.setMinIntersectionMagnitude(0.05)
constraintsOnSelf.setMinIntersectionDistance(0.05)

# with respect to entities of the other set, we want to have larger intersection angles
constraintsOnOther = EntityNetworkConstraints()
constraintsOnOther.setMinDistance(0.05)
Dennis Gläser's avatar
Dennis Gläser committed
53
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
constraintsOnOther.setMinIntersectionMagnitude(0.05)
constraintsOnOther.setMinIntersectionDistance(0.05)

# we use unique identifiers for the entities of the two orientations
from frackit.common import Id
idSet1 = Id(1)
idSet2 = Id(2)

# use the status class to define when to stop sampling
from frackit.sampling import SamplingStatus
status = SamplingStatus()
status.setTargetCount(idSet1, 15) # we want 15 entities in each set
status.setTargetCount(idSet2, 15) # we want 15 entities in each set

# start sampling into set 1 and keep alternating
sampleIntoSet1 = True

# lists to store the sampled entities
entities1 = []
entities2 = []

print("\n --- Start entity sampling ---\n")
while not status.finished():
    # sample a quadrilateral, alternating between sampler 1 and sampler 2
78
    quad = quadSampler1() if sampleIntoSet1 else quadSampler2()
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
    entitySet = entities1 if sampleIntoSet1 else entities2
    otherEntitySet = entities2 if sampleIntoSet1 else entities1

    # skip the rest if constraints are violated
    if not constraintsOnSelf.evaluate(entitySet, quad):
        status.increaseRejectedCounter()
        continue

    # skip the rest if constraints are violated
    if not constraintsOnOther.evaluate(otherEntitySet, quad):
        status.increaseRejectedCounter()
        continue

    # the entity is admissible
    entitySet.append(quad)
    id = idSet1 if sampleIntoSet1 else idSet2
    status.increaseCounter(id)
    status.print()

    # sample into other set the next time
    sampleIntoSet1 = not sampleIntoSet1

print("\n --- Finished entity sampling ---\n")

# We can now create an entity network from the two sets
from frackit.entitynetwork import EntityNetworkBuilder
builder = EntityNetworkBuilder()
Dennis Gläser's avatar
Dennis Gläser committed
106
107
builder.addEntities(entities1)
builder.addEntities(entities2)
108
109
110

# let the builder construct the network and write it to gmsh file format
print("\n --- Constructing entity network from the raw entities ---\n")
Dennis Gläser's avatar
Dennis Gläser committed
111
network = builder.build()
112
113
114

print("\n --- Writing .geo file ---\n")
from frackit.io import GmshWriter
Dennis Gläser's avatar
Dennis Gläser committed
115
writer = GmshWriter(network)
116
117
118
writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
writer.write("network") # filename of the .geo files (will add extension .geo automatically)

119
print("\n --- Finished writing .geo file ---\n")