example3.py 10.5 KB
Newer Older
1
2
import sys

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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)
19
if len(solids) != 3: sys.exit("Expected the .brep file to contain 3 solids")
20
21
22

# 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.
23
from frackit.geometry import computeMagnitude
24
25
26
networkDomain = solids[1]
domainVolume = computeMagnitude(networkDomain);
shells = getShells(networkDomain)
27
if len(shells) != 1: sys.exit("Expected a single shell bounding the domain")
28
29

domainBoundaryFaces = getFaces(shells[0])
30
if len(domainBoundaryFaces) != 6: sys.exit("Expected 6 faces to bound the domain")
31
32
33
34
35
36
37
38
39


##############################################################################
## 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
40
from frackit.geometry import Box, OCCShapeWrapper
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
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
Dennis Gläser's avatar
Dennis Gläser committed
68
69
                                   uniformSampler(30.0, 60.0),                      # strike length
                                   uniformSampler(30.0, 60.0))                      # dip length
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

# 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()
Dennis Gläser's avatar
Dennis Gläser committed
107
108
constraintsMatrix.addConstraints(constraints1,           # constraint instance
                                 [diskSetId, diskSetId]) # sets between which to use these constraints
109

Dennis Gläser's avatar
Dennis Gläser committed
110
111
constraintsMatrix.addConstraints(constraints2,           # constraint instance
                                 [quadSetId, quadSetId]) # sets between which to use these constraints
112

Dennis Gläser's avatar
Dennis Gläser committed
113
114
115
constraintsMatrix.addConstraints(constraintsOnOther,       # constraint instance
                                 [[diskSetId, quadSetId],  # sets between which to use these constraints
                                  [quadSetId, diskSetId]]) # sets between which to use these constraints
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

# 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
Dennis Gläser's avatar
Dennis Gläser committed
134
entitySets = {diskSetId: [], quadSetId: []}
135
136
137
138
139
140

# Alternate between set 1 & set 2 during sampling phase
sampleIntoSet1 = True
containedNetworkArea = 0.0;
while not status.finished():
    id = diskSetId if sampleIntoSet1 else quadSetId
141
    geom = diskSampler() if sampleIntoSet1 else quadSampler()
142
143
144
145
146
147
148
149

    # 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²)
150
    from frackit.geometry import computeContainedMagnitude
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    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
Dennis Gläser's avatar
Dennis Gläser committed
167
    entitySets[id].append(geom)
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    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.
Dennis Gläser's avatar
Dennis Gläser committed
188
print("Building and writing contained, confined network")
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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
Dennis Gläser's avatar
Dennis Gläser committed
209
print("Building and writing contained, unconfined network")
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
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
Dennis Gläser's avatar
Dennis Gläser committed
230
print("Building and writing uncontained, unconfined network")
231
232
233
234
235
uncontainedBuilder.clear();
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))

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