example3.py 10.7 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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
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
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
# print welcome message
print("\n\n"
      "##############################################################################\n"
      "## Example: Creation of entities (disks/quadrilaterals) in a layered medium ##\n"
      "##############################################################################\n"
      "\n\n")

####################################################
# 1. Read in the domain geometry from .brep file. ##
#    The file name is defined in CMakeLists.txt   ##
####################################################
from frackit.occutilities import readShape, getSolids, getShells, getFaces, getBoundingBox
domainShape = readShape("layers.brep")

# obtain the three solids contained in the file
solids = getSolids(domainShape)
if len(solids) != 3: raise RuntimeError("Expected the .brep file to contain 3 solids")

# The sub-domain we want to create a network in is the center one.
# Compute its volume and get the boundary faces for constraint evaluation.
from frackit.magnitude import computeMagnitude
networkDomain = solids[1]
domainVolume = computeMagnitude(networkDomain);
shells = getShells(networkDomain)
if len(shells) != 1: raise RuntimeError("Expected a single shell bounding the domain")

domainBoundaryFaces = getFaces(shells[0])
if len(domainBoundaryFaces) != 6: raise RuntimeError("Expected 6 faces to bound the domain")


##############################################################################
## 2. Make sampler classes to randomly sample points (entity centers)       ##
##    and disk/quadrilaterals as entities (using the sampled center points) ##
##############################################################################

# Bounding box of the domain in which we want to place the entities
# In order for this to work we parse the solid back into a wrapper around a BRep shape
from frackit.geometry import Box
from frackit.occutilities import OCCShapeWrapper
domainBBox = getBoundingBox( OCCShapeWrapper(networkDomain) )

# creates a uniform sampler within an interval
import random
def uniformSampler(a, b):
    def sample(): return random.uniform(a, b)
    return sample

# creates a sampler from a normal distribution with the given mean and std deviation
def normalSampler(mean, stdDev):
    def sample(): return random.gauss(mean, stdDev)
    return sample

from frackit.common import toRadians, Id
from frackit.sampling import DiskSampler, QuadrilateralSampler, makeUniformPointSampler
# sampler for disks (orientation 1)
diskSampler = DiskSampler(makeUniformPointSampler(domainBBox),           # sampler for disk center points
                          uniformSampler(30.0, 6.5),                     # major axis length: mean value & standard deviation
                          uniformSampler(24.0, 4.5),                     # minor axis length: mean value & standard deviation
                          normalSampler(toRadians(0.0), toRadians(7.5)), # rotation around x-axis: mean value & standard deviation
                          normalSampler(toRadians(0.0), toRadians(7.5)), # rotation around y-axis: mean value & standard deviation
                          normalSampler(toRadians(0.0), toRadians(7.5))) # rotation around z-axis: mean value & standard deviation

# sampler for quadrilaterals (orientation 2)
quadSampler = QuadrilateralSampler(makeUniformPointSampler(domainBBox),             # sampler for quadrilateral center points
                                   normalSampler(toRadians(45.0), toRadians(5.0)),  # strike angle: mean value & standard deviation
                                   normalSampler(toRadians(90.0), toRadians(5.0)),  # dip angle: mean value & standard deviation
                                   normalSampler(45.0, 5.0),                        # edge length: mean value & standard deviation
                                   5.0)                                             # threshold for minimum edge length

# Define ids for the two entity sets
diskSetId = Id(1) # we give the set of orientation one, consisting of disks, the id 1
quadSetId = Id(2) # we give the set of orientation two, consisting of quadrilaterals, the id 2


#######################################################################
## 3. Define constraints that should be fulfilled among the entities ##
##    of different orientations.                                     ##
#######################################################################

from frackit.entitynetwork import EntityNetworkConstraints, EntityNetworkConstraintsMatrix

# constructs default constraints for this test to be modified after
def makeDefaultConstraints():
    c = EntityNetworkConstraints()
    c.setMinDistance(2.5)
    c.setMinIntersectingAngle(toRadians(25.0))
    c.setMinIntersectionMagnitude(5.0)
    c.setMinIntersectionDistance(2.5)
    return c

# Define constraints between entities of orientation 1
constraints1 = makeDefaultConstraints()

# Define constraints between entities of orientation 2
# we want to enforce larger spacing between those entities
constraints2 = makeDefaultConstraints()
constraints2.setMinDistance(5.0)

# Define constraints between entities of different sets
constraintsOnOther = makeDefaultConstraints()
constraintsOnOther.setMinDistance(2.5)
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0))

# We can use the constraints matrix to facilitate constraint evaluation
constraintsMatrix = EntityNetworkConstraintsMatrix()
constraintsMatrix.addConstraints(constraints1,                       # constraint instance
                                 [diskSetId.get(), diskSetId.get()]) # sets between which to use these constraints

constraintsMatrix.addConstraints(constraints2,                       # constraint instance
                                 [quadSetId.get(), quadSetId.get()]) # sets between which to use these constraints

constraintsMatrix.addConstraints(constraintsOnOther,                   # constraint instance
                                 [[diskSetId.get(), quadSetId.get()],  # sets between which to use these constraints
                                  [quadSetId.get(), diskSetId.get()]]) # sets between which to use these constraints

# Moreover, we define constraints w.r.t. the domain boundary
constraintsOnDomain = makeDefaultConstraints()
constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0))


###########################
## 4. Network generation ##
###########################

# Helper class for terminal output of the creation
# progress and definition of stop criterion etc
from frackit.sampling import SamplingStatus
status = SamplingStatus()
status.setTargetCount(diskSetId, 12) # we want 11 entities of orientation 1
status.setTargetCount(quadSetId, 16) # we want 13 entities of orientation 2

# store all entity sets in a dictionary
entitySets = {diskSetId.get(): [], quadSetId.get(): []}

# Alternate between set 1 & set 2 during sampling phase
sampleIntoSet1 = True
containedNetworkArea = 0.0;
while not status.finished():
    id = diskSetId if sampleIntoSet1 else quadSetId
    geom = diskSampler.sample() if sampleIntoSet1 else quadSampler.sample()

    # If the set this geometry belongs to is finished, skip the rest
    if status.finished(id):
        sampleIntoSet1 = not sampleIntoSet1
        status.increaseRejectedCounter()
        continue

    # Moreover, we want to avoid small fragments (< 250 m²)
    from frackit.magnitude import computeContainedMagnitude
    containedArea = computeContainedMagnitude(geom, networkDomain);
    if containedArea < 350.0:
        status.increaseRejectedCounter()
        continue

    # enforce constraints w.r.t. to the other entities
    if not constraintsMatrix.evaluate(entitySets, geom, id):
        status.increaseRejectedCounter()
        continue

    # enforce constraints w.r.t. the domain boundaries
    if not constraintsOnDomain.evaluate(domainBoundaryFaces, geom):
        status.increaseRejectedCounter()
        continue

    # the geometry is admissible
    entitySets[id.get()].append(geom)
    status.increaseCounter(id)
    status.print()

    # keep track of entity area
    containedNetworkArea += containedArea;

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

# print the final entity density
density = containedNetworkArea/domainVolume;
print("\nEntity density of the contained network: {:f} m²/m³\n".format(density))


##########################################################################
## 5. The entities of the network have been created. We can now         ##
##    construct different types of networks (contained, confined, etc.) ##
##########################################################################

# construct and write a contained network, i.e. write out both network and domain.
print("Building and writing contained, confined network\n")
from frackit.entitynetwork import EntityNetworkBuilder, ContainedEntityNetworkBuilder
builder = ContainedEntityNetworkBuilder();

# add sub-domains
builder.addConfiningSubDomain(solids[0],     Id(1));
builder.addConfiningSubDomain(networkDomain, Id(2));
builder.addConfiningSubDomain(solids[2],     Id(3));

# The entites that we created all belong to subdomain 2
for setId in entitySets: builder.addSubDomainEntities(entitySets[setId], Id(2))

# now we can build and write out the network in Gmsh file format
from frackit.io import GmshWriter
gmshWriter = GmshWriter(builder.build());
gmshWriter.write("contained_confined", # body of the filename to be used (will add .geo)
                 2.5,                  # mesh size to be used on entities
                 5.0);                 # mesh size to be used on domain boundaries

# we can also not confine the network to its sub-domain,
# simply by adding the sub-domains as non-confining
print("Building and writing contained, unconfined network\n")
builder.clear();
builder.addSubDomain(solids[0],     Id(1));
builder.addSubDomain(networkDomain, Id(2));
builder.addSubDomain(solids[2],     Id(3));
for setId in entitySets: builder.addSubDomainEntities(entitySets[setId], Id(2))

gmshWriter = GmshWriter(builder.build());
gmshWriter.write("contained_unconfined", 2.5, 5.0);

# We could also only write out the network, without the domain
# For example, confining the network to the sub-domain...
print("Building and writing uncontained, confined network")
uncontainedBuilder = EntityNetworkBuilder()
uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2));
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))

gmshWriter = GmshWriter(uncontainedBuilder.build());
gmshWriter.write("uncontained_confined", 2.5);

# ... or not confining it
print("Building and writing uncontained, unconfined network\n")
uncontainedBuilder.clear();
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))

gmshWriter = GmshWriter(uncontainedBuilder.build());
gmshWriter.write("uncontained_unconfined", 2.5);