example2.py 5.91 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
import frackit.geometry as geometry

# we want to use a cylindrical domain with radius=0.5 and height=1.0
domain = geometry.Cylinder(0.5, 1.0)

# 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
box = geometry.Box(-0.5, -0.5, 0.0, 0.5, 0.5, 1.0)
pointSampler = makeUniformPointSampler(box)

# returns a sampler from a gaussian distribution with mean and std deviation
def gaussianSampler(mean, stdDev):
    import random
    def sample():
        return random.gauss(mean, stdDev)
    return sample

# we sample quadrialeterals within the box with gaussian distributions for the parameters
from frackit.common import toRadians
from frackit.sampling import QuadrilateralSampler as QuadSampler
quadSampler1 = QuadSampler(pointSampler,
                           gaussianSampler(toRadians(45.0), toRadians(5.0)), # strike angle
                           gaussianSampler(toRadians(90.0), toRadians(5.0)), # dip angle
                           gaussianSampler(0.5, 0.1),                        # edge length
                           0.05)                                             # threshold for minimum edge length

# sampler for quadrilaterals of the secondary orientation
quadSampler2 = QuadSampler(pointSampler,
                           gaussianSampler(toRadians(0.0), toRadians(5.0)), # strike angle
                           gaussianSampler(toRadians(0.0), toRadians(5.0)), # dip angle
                           gaussianSampler(0.5, 0.1),                       # edge length
                           0.05)                                            # threshold for minimum edge length

# constructs a constraints object with default settings for this example
from frackit.entitynetwork import EntityNetworkConstraints
def makeDefaultConstraints():
    c = EntityNetworkConstraints()
    c.setMinDistance(0.05)
    c.setMinIntersectingAngle(toRadians(30.0))
    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()
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));

# 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
    quad = quadSampler1.sample() if sampleIntoSet1 else quadSampler2.sample()
    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
    if not constraintsOnBoundary.evaluate(getShape(domain.bottomFace()), quad):
        status.increaseRejectedCounter()
        continue

    # reject entities whose contained part is too small
    from frackit.magnitude import computeContainedMagnitude
    containedArea = computeContainedMagnitude(quad, domain);
    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
builder.addConfiningSubDomain(domain, Id(1));

# define all entities to be embedded in this domain
builder.addSubDomainEntities(entities1, Id(1));
builder.addSubDomainEntities(entities2, Id(1));

# let the builder construct the network and write it to gmsh file format
print("\n --- Constructing entity network from the raw entities ---\n")
network = builder.build();

print("\n --- Writing .geo file ---\n")
from frackit.io import GmshWriter
writer = GmshWriter(network);
writer.write("network", # filename of the .geo files (will add extension .geo automatically)
             0.1);      # element size to be used
print("\n --- Finished writing .geo file ---\n")