example3.py 10.3 KB
Newer Older
1
2
import sys

3
4
5
6
7
8
9
10
11
12
# 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. ##
####################################################
Dennis Gläser's avatar
Dennis Gläser committed
13
from frackit.occutilities import readShape, getSolids, getShells, getFaces
14
15
16
17
domainShape = readShape("layers.brep")

# obtain the three solids contained in the file
solids = getSolids(domainShape)
18
if len(solids) != 3: sys.exit("Expected the .brep file to contain 3 solids")
19
20
21

# 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.
22
from frackit.geometry import computeMagnitude
23
networkDomain = solids[1]
Dennis Gläser's avatar
Dennis Gläser committed
24
25
26
27
domainVolume = computeMagnitude(networkDomain)
domainShell = getShells(networkDomain)
if len(domainShell) != 1: sys.exit("Expected a single shell bounding the domain")
domainBoundaryFaces = getFaces(domainShell[0])
28
29
30
31
32
33
34
35


##############################################################################
## 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
Dennis Gläser's avatar
Dennis Gläser committed
36
37
from frackit.geometry import OCCShapeWrapper, getBoundingBox
domainBBox = getBoundingBox(networkDomain)
38
39
40
41
42
43
44
45
46
47
48
49

# 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

Dennis Gläser's avatar
Dennis Gläser committed
50
51
import math
from frackit.common import Id
52
from frackit.sampling import DiskSampler, QuadrilateralSampler, makeUniformPointSampler
Dennis Gläser's avatar
Dennis Gläser committed
53

54
# sampler for disks (orientation 1)
Dennis Gläser's avatar
Dennis Gläser committed
55
56
57
58
59
60
diskSampler = DiskSampler(pointSampler   = makeUniformPointSampler(domainBBox),
                          majAxisSampler = uniformSampler(30.0, 6.5),
                          minAxisSampler = uniformSampler(24.0, 4.5),
                          xAngleSampler  = normalSampler(math.radians(0.0), math.radians(7.5)),
                          yAngleSampler  = normalSampler(math.radians(0.0), math.radians(7.5)),
                          zAngleSampler  = normalSampler(math.radians(0.0), math.radians(7.5)))
61
62

# sampler for quadrilaterals (orientation 2)
Dennis Gläser's avatar
Dennis Gläser committed
63
64
65
66
67
quadSampler = QuadrilateralSampler(pointSampler        = makeUniformPointSampler(domainBBox),
                                   strikeAngleSampler  = normalSampler(math.radians(45.0), math.radians(5.0)),
                                   dipAngleSampler     = normalSampler(math.radians(90.0), math.radians(5.0)),
                                   strikeLengthSampler = uniformSampler(30.0, 60.0),
                                   dipLengthSampler    = uniformSampler(30.0, 60.0))
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

# 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)
Dennis Gläser's avatar
Dennis Gläser committed
85
    c.setMinIntersectingAngle(math.radians(25.0))
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    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)
Dennis Gläser's avatar
Dennis Gläser committed
101
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
102
103
104

# We can use the constraints matrix to facilitate constraint evaluation
constraintsMatrix = EntityNetworkConstraintsMatrix()
Dennis Gläser's avatar
Dennis Gläser committed
105
106
constraintsMatrix.addConstraints(constraints1,           # constraint instance
                                 [diskSetId, diskSetId]) # sets between which to use these constraints
107

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

Dennis Gläser's avatar
Dennis Gläser committed
111
112
113
constraintsMatrix.addConstraints(constraintsOnOther,       # constraint instance
                                 [[diskSetId, quadSetId],  # sets between which to use these constraints
                                  [quadSetId, diskSetId]]) # sets between which to use these constraints
114
115
116

# Moreover, we define constraints w.r.t. the domain boundary
constraintsOnDomain = makeDefaultConstraints()
Dennis Gläser's avatar
Dennis Gläser committed
117
constraintsOnDomain.setMinIntersectingAngle(math.radians(15.0))
118
119
120
121
122
123
124
125
126
127


###########################
## 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()
Dennis Gläser's avatar
Dennis Gläser committed
128
129
status.setTargetCount(diskSetId, 12) # number of desired entities of orientation 1
status.setTargetCount(quadSetId, 16) # number of desired entities of orientation 2
130
131

# store all entity sets in a dictionary
Dennis Gläser's avatar
Dennis Gläser committed
132
entitySets = {diskSetId: [], quadSetId: []}
133
134
135

# Alternate between set 1 & set 2 during sampling phase
sampleIntoSet1 = True
Dennis Gläser's avatar
Dennis Gläser committed
136
containedNetworkArea = 0.0
137
138
while not status.finished():
    id = diskSetId if sampleIntoSet1 else quadSetId
139
    geom = diskSampler() if sampleIntoSet1 else quadSampler()
140
141
142
143

    # If the set this geometry belongs to is finished, skip the rest
    if status.finished(id):
        sampleIntoSet1 = not sampleIntoSet1
Dennis Gläser's avatar
Dennis Gläser committed
144
        status.increaseRejectedCounter("set finished")
145
146
        continue

Dennis Gläser's avatar
Dennis Gläser committed
147
    # Moreover, we want to avoid small fragments
148
    from frackit.geometry import computeContainedMagnitude
Dennis Gläser's avatar
Dennis Gläser committed
149
150
    containedArea = computeContainedMagnitude(geom, networkDomain)
    if containedArea < 250.0:
Dennis Gläser's avatar
Dennis Gläser committed
151
        status.increaseRejectedCounter("minimum contained area violation")
152
153
154
        continue

    # enforce constraints w.r.t. to the other entities
155
156
157
    eval = constraintsMatrix.evaluate(entitySets, geom, id)
    if eval.violationDetected():
        status.increaseRejectedCounter(eval.violationLabel())
158
159
160
        continue

    # enforce constraints w.r.t. the domain boundaries
161
162
163
    eval = constraintsOnDomain.evaluate(domainBoundaryFaces, geom)
    if eval.violationDetected():
        status.increaseRejectedCounter(eval.violationLabel() + " (boundary)")
164
165
166
        continue

    # the geometry is admissible
Dennis Gläser's avatar
Dennis Gläser committed
167
    entitySets[id].append(geom)
168
169
170
171
    status.increaseCounter(id)
    status.print()

    # keep track of entity area
Dennis Gläser's avatar
Dennis Gläser committed
172
    containedNetworkArea += containedArea
173
174
175
176
177

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

# print the final entity density
Dennis Gläser's avatar
Dennis Gläser committed
178
density = containedNetworkArea/domainVolume
179
180
print("\nEntity density of the contained network: {:f} m²/m³\n".format(density))

Dennis Gläser's avatar
Dennis Gläser committed
181
# print info on rejection events
Dennis Gläser's avatar
Dennis Gläser committed
182
status.printRejectionData()
Dennis Gläser's avatar
Dennis Gläser committed
183
print("")
184
185
186
187
188
189
190

##########################################################################
## 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
191
print("Building and writing contained, confined network")
192
from frackit.entitynetwork import EntityNetworkBuilder, ContainedEntityNetworkBuilder
Dennis Gläser's avatar
Dennis Gläser committed
193
builder = ContainedEntityNetworkBuilder()
194
195

# add sub-domains
Dennis Gläser's avatar
Dennis Gläser committed
196
197
198
builder.addConfiningSubDomain(solids[0],     Id(1))
builder.addConfiningSubDomain(networkDomain, Id(2))
builder.addConfiningSubDomain(solids[2],     Id(3))
199
200
201
202
203
204

# 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
Dennis Gläser's avatar
Dennis Gläser committed
205
gmshWriter = GmshWriter(builder.build())
206
207
208
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
gmshWriter.write("contained_confined") # body of the filename to be used (will add .geo)
209
210
211

# 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
212
print("Building and writing contained, unconfined network")
Dennis Gläser's avatar
Dennis Gläser committed
213
214
215
216
builder.clear()
builder.addSubDomain(solids[0],     Id(1))
builder.addSubDomain(networkDomain, Id(2))
builder.addSubDomain(solids[2],     Id(3))
217
218
for setId in entitySets: builder.addSubDomainEntities(entitySets[setId], Id(2))

Dennis Gläser's avatar
Dennis Gläser committed
219
gmshWriter = GmshWriter(builder.build())
220
221
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
Dennis Gläser's avatar
Dennis Gläser committed
222
gmshWriter.write("contained_unconfined")
223
224
225
226
227

# 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()
Dennis Gläser's avatar
Dennis Gläser committed
228
uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2))
229
230
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))

Dennis Gläser's avatar
Dennis Gläser committed
231
gmshWriter = GmshWriter(uncontainedBuilder.build())
232
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
Dennis Gläser's avatar
Dennis Gläser committed
233
gmshWriter.write("uncontained_confined")
234
235

# ... or not confining it
Dennis Gläser's avatar
Dennis Gläser committed
236
print("Building and writing uncontained, unconfined network")
Dennis Gläser's avatar
Dennis Gläser committed
237
uncontainedBuilder.clear()
238
239
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))

Dennis Gläser's avatar
Dennis Gläser committed
240
gmshWriter = GmshWriter(uncontainedBuilder.build())
241
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
Dennis Gläser's avatar
Dennis Gläser committed
242
gmshWriter.write("uncontained_unconfined")