example3.py 10.9 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

    # 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
146
        status.increaseRejectedCounter("set finished")
147
148
149
        continue

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

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

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

    # the geometry is admissible
Dennis Gläser's avatar
Dennis Gläser committed
169
    entitySets[id].append(geom)
170
171
172
173
174
175
176
177
178
179
180
181
182
    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))

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

##########################################################################
## 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
193
print("Building and writing contained, confined network")
194
195
196
197
198
199
200
201
202
203
204
205
206
207
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());
208
209
210
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)
211
212
213

# 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
214
print("Building and writing contained, unconfined network")
215
216
217
218
219
220
221
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());
222
223
224
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
gmshWriter.write("contained_unconfined");
225
226
227
228
229
230
231
232
233

# 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());
234
235
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.write("uncontained_confined");
236
237

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

gmshWriter = GmshWriter(uncontainedBuilder.build());
243
244
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.write("uncontained_unconfined");