example2.py 6.07 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

Dennis Gläser's avatar
Dennis Gläser committed
5
domain = geometry.Cylinder(radius=0.5, height=1.0)
6
7
8
9

# we sample points uniformly within the shifted unit cube such
# that the origin is in the center of the bottom boundary of the domain
from frackit.sampling import makeUniformPointSampler
Dennis Gläser's avatar
Dennis Gläser committed
10
11
box = geometry.Box(xmin=-0.5, ymin=-0.5, zmin=0.0,
                   xmax=0.5,  ymax=0.5,  zmax=1.0)
12
13
14
15
pointSampler = makeUniformPointSampler(box)

# 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
16
    def sample(): return random.gauss(mean, stdDev)
17
18
    return sample

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

24
25
# we sample quadrialeterals within the box with gaussian distributions for the parameters
from frackit.sampling import QuadrilateralSampler as QuadSampler
Dennis Gläser's avatar
Dennis Gläser committed
26
27
28
29
30
quadSampler1 = QuadSampler(pointSampler        = pointSampler,
                           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.6),
                           dipLengthSampler    = uniformSampler(0.4, 0.6))
31
32

# sampler for quadrilaterals of the secondary orientation
Dennis Gläser's avatar
Dennis Gläser committed
33
34
35
36
37
quadSampler2 = QuadSampler(pointSampler        = pointSampler,
                           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.6),
                           dipLengthSampler    = uniformSampler(0.4, 0.6))
38
39
40
41
42
43

# constructs a constraints object with default settings for this example
from frackit.entitynetwork import EntityNetworkConstraints
def makeDefaultConstraints():
    c = EntityNetworkConstraints()
    c.setMinDistance(0.05)
Dennis Gläser's avatar
Dennis Gläser committed
44
    c.setMinIntersectingAngle(math.radians(30.0))
45
46
47
48
49
50
51
52
53
54
55
56
57
58
    c.setMinIntersectionMagnitude(0.05)
    c.setMinIntersectionDistance(0.05)
    return c

# 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.
constraintsOnSelf = makeDefaultConstraints()

# with respect to entities of the other set, we want to have larger intersection angles
constraintsOnOther = makeDefaultConstraints()
Dennis Gläser's avatar
Dennis Gläser committed
59
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
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

# we use the default constraints w.r.t. to the domain boundary
constraintsOnBoundary = makeDefaultConstraints()

# 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
85
    quad = quadSampler1() if sampleIntoSet1 else quadSampler2()
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    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

    # enforce constraints w.r.t. to the domain boundary
    if not constraintsOnBoundary.evaluate(domain.topFace(), quad):
        status.increaseRejectedCounter()
        continue
    if not constraintsOnBoundary.evaluate(domain.bottomFace(), quad):
        status.increaseRejectedCounter()
        continue

    from frackit.occutilities import getShape
108
    if not constraintsOnBoundary.evaluate(getShape(domain.lateralFace()), quad):
109
110
111
112
        status.increaseRejectedCounter()
        continue

    # reject entities whose contained part is too small
113
    from frackit.geometry import computeContainedMagnitude
Dennis Gläser's avatar
Dennis Gläser committed
114
    containedArea = computeContainedMagnitude(quad, domain)
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
    if (containedArea < 0.01):
        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 ContainedEntityNetworkBuilder
builder = ContainedEntityNetworkBuilder()

# our domain confines all entities and receives the id 1
Dennis Gläser's avatar
Dennis Gläser committed
135
builder.addConfiningSubDomain(domain, Id(1))
136
137

# define all entities to be embedded in this domain
Dennis Gläser's avatar
Dennis Gläser committed
138
139
builder.addSubDomainEntities(entities1, Id(1))
builder.addSubDomainEntities(entities2, Id(1))
140
141
142

# 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
143
network = builder.build()
144
145
146

print("\n --- Writing .geo file ---\n")
from frackit.io import GmshWriter
Dennis Gläser's avatar
Dennis Gläser committed
147
writer = GmshWriter(network)
148
writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
149
writer.setMeshSize(GmshWriter.GeometryTag.subDomain, 0.1)
150
writer.write("network") # filename of the .geo files (will add extension .geo automatically)
151
print("\n --- Finished writing .geo file ---\n")