Commit 6e894a3f authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/simplify-examples' into 'master'

Feature/cleanup-examples

See merge request tools/frackit!206
parents 95f6ab29 3ed481ef
Pipeline #2846 passed with stages
in 9 minutes and 49 seconds
......@@ -36,11 +36,10 @@ int main(int argc, char** argv)
// These points will be used as the quadrilateral centers.
auto pointSampler = makeUniformPointSampler(domain);
// Sampler class for quadrilaterals. Per default, this uses
// uniform distributions for all parameters defining the quadrilaterals.
// Quadrilateral samplers require distributions for strike angle, dip angle,
// edge length (see class description for more details). Moreover, we define
// a minimum edge length.
// Sampler class for quadrilaterals. Per default, this uses uniform distributions
// for all parameters defining the quadrilaterals. Quadrilateral samplers require
// distributions for strike angle, dip angle, edge length (see class description
// for more details).
using QuadSampler = QuadrilateralSampler<worldDimension>;
using NormalDistro = std::normal_distribution<ctype>;
using UniformDistro = std::uniform_real_distribution<ctype>;
......
import frackit.geometry as geometry
import random
import math
# we use the unit cube as domain
box = geometry.Box(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
box = geometry.Box(xmin=0.0, ymin=0.0, zmin=0.0,
xmax=1.0, ymax=1.0, zmax=1.0)
# we sample points uniformly within the domain
from frackit.sampling import makeUniformPointSampler
# 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)
def sample(): return random.gauss(mean, stdDev)
return sample
# returns a sampler from a uniform distribution between min and max
def uniformSampler(min, max):
import random
def sample():
return random.uniform(min, max)
def sample(): return random.uniform(min, max)
return sample
# we sample quadrialeterals within the box with gaussian distributions for the parameters
from frackit.common import toRadians
# we sample quadrilaterals within the box with gaussian distributions for the parameters
from frackit.sampling import QuadrilateralSampler as QuadSampler
quadSampler1 = QuadSampler(makeUniformPointSampler(box),
gaussianSampler(toRadians(45.0), toRadians(5.0)), # strike angle
gaussianSampler(toRadians(90.0), toRadians(5.0)), # dip angle
uniformSampler(0.4, 0.8), # strike length
uniformSampler(0.4, 0.8)) # dip length
quadSampler1 = QuadSampler(pointSampler = makeUniformPointSampler(box),
strikeAngleSampler = gaussianSampler(math.radians(45.0), math.radians(5.0)),
dipAngleSampler = gaussianSampler(math.radians(90.0), math.radians(5.0)),
strikeLengthSampler = uniformSampler(0.4, 0.8),
dipLengthSampler = uniformSampler(0.4, 0.8))
# sampler for quadrilaterals of the secondary orientation
quadSampler2 = QuadSampler(makeUniformPointSampler(box),
gaussianSampler(toRadians(0.0), toRadians(5.0)), # strike angle
gaussianSampler(toRadians(0.0), toRadians(5.0)), # dip angle
uniformSampler(0.4, 0.8), # strike length
uniformSampler(0.4, 0.8)) # dip length
quadSampler2 = QuadSampler(pointSampler = makeUniformPointSampler(box),
strikeAngleSampler = gaussianSampler(math.radians(0.0), math.radians(5.0)),
dipAngleSampler = gaussianSampler(math.radians(0.0), math.radians(5.0)),
strikeLengthSampler = uniformSampler(0.4, 0.8),
dipLengthSampler = uniformSampler(0.4, 0.8))
# We want to enforce some constraints on the set of quadrilaterals.
# In particular, for entities of the same set we want a minimum spacing
......@@ -45,14 +43,14 @@ quadSampler2 = QuadSampler(makeUniformPointSampler(box),
from frackit.entitynetwork import EntityNetworkConstraints
constraintsOnSelf = EntityNetworkConstraints()
constraintsOnSelf.setMinDistance(0.05)
constraintsOnSelf.setMinIntersectingAngle(toRadians(30.0))
constraintsOnSelf.setMinIntersectingAngle(math.radians(30.0))
constraintsOnSelf.setMinIntersectionMagnitude(0.05)
constraintsOnSelf.setMinIntersectionDistance(0.05)
# with respect to entities of the other set, we want to have larger intersection angles
constraintsOnOther = EntityNetworkConstraints()
constraintsOnOther.setMinDistance(0.05)
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
constraintsOnOther.setMinIntersectionMagnitude(0.05)
constraintsOnOther.setMinIntersectionDistance(0.05)
......@@ -105,16 +103,16 @@ print("\n --- Finished entity sampling ---\n")
# We can now create an entity network from the two sets
from frackit.entitynetwork import EntityNetworkBuilder
builder = EntityNetworkBuilder()
builder.addEntities(entities1);
builder.addEntities(entities2);
builder.addEntities(entities1)
builder.addEntities(entities2)
# 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();
network = builder.build()
print("\n --- Writing .geo file ---\n")
from frackit.io import GmshWriter
writer = GmshWriter(network);
writer = GmshWriter(network)
writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
writer.write("network") # filename of the .geo files (will add extension .geo automatically)
......
......@@ -67,7 +67,7 @@ raw entities. This can potentially fail to detect small length scales that appea
the entities to the domain. A possibilty to address this and make the algorithm more robust
(but computationally more demanding) is discussed at the end of this page.
As you can see in the above code snippet, the `Cylinder` class provides functions for obtaining th
As you can see in the above code snippet, the `Cylinder` class provides functions for obtaining the
representations of the top, bottom and lateral boundaries. The top and bottom boundaries are
represented by instances of the `Disk` class, while the lateral surface is described by the class
`CylinderSurface`. Please also note that for the lateral surface, we enforce the constraints
......
import frackit.geometry as geometry
import random
import math
# we want to use a cylindrical domain with radius=0.5 and height=1.0
domain = geometry.Cylinder(0.5, 1.0)
domain = geometry.Cylinder(radius=0.5, height=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)
box = geometry.Box(xmin=-0.5, ymin=-0.5, zmin=0.0,
xmax=0.5, ymax=0.5, zmax=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)
def sample(): return random.gauss(mean, stdDev)
return sample
#returns a sampler from a uniform distribution between min and max
def uniformSampler(min, max):
import random
def sample():
return random.uniform(min, max)
def sample(): return random.uniform(min, max)
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
uniformSampler(0.4, 0.6), # strike length
uniformSampler(0.4, 0.6)) # dip length
quadSampler1 = QuadSampler(pointSampler = pointSampler,
strikeAngleSampler = gaussianSampler(math.radians(45.0), math.radians(5.0)),
dipAngleSampler = gaussianSampler(math.radians(90.0), math.radians(5.0)),
strikeLengthSampler = uniformSampler(0.4, 0.6),
dipLengthSampler = uniformSampler(0.4, 0.6))
# 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
uniformSampler(0.4, 0.6), # strike length
uniformSampler(0.4, 0.6)) # dip length
quadSampler2 = QuadSampler(pointSampler = pointSampler,
strikeAngleSampler = gaussianSampler(math.radians(0.0), math.radians(5.0)),
dipAngleSampler = gaussianSampler(math.radians(0.0), math.radians(5.0)),
strikeLengthSampler = uniformSampler(0.4, 0.6),
dipLengthSampler = uniformSampler(0.4, 0.6))
# 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.setMinIntersectingAngle(math.radians(30.0))
c.setMinIntersectionMagnitude(0.05)
c.setMinIntersectionDistance(0.05)
return c
......@@ -59,7 +56,7 @@ constraintsOnSelf = makeDefaultConstraints()
# with respect to entities of the other set, we want to have larger intersection angles
constraintsOnOther = makeDefaultConstraints()
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
# we use the default constraints w.r.t. to the domain boundary
constraintsOnBoundary = makeDefaultConstraints()
......@@ -114,7 +111,7 @@ while not status.finished():
# reject entities whose contained part is too small
from frackit.geometry import computeContainedMagnitude
containedArea = computeContainedMagnitude(quad, domain);
containedArea = computeContainedMagnitude(quad, domain)
if (containedArea < 0.01):
status.increaseRejectedCounter()
continue
......@@ -135,19 +132,19 @@ from frackit.entitynetwork import ContainedEntityNetworkBuilder
builder = ContainedEntityNetworkBuilder()
# our domain confines all entities and receives the id 1
builder.addConfiningSubDomain(domain, 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));
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();
network = builder.build()
print("\n --- Writing .geo file ---\n")
from frackit.io import GmshWriter
writer = GmshWriter(network);
writer = GmshWriter(network)
writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1)
writer.setMeshSize(GmshWriter.GeometryTag.subDomain, 0.1)
writer.write("network") # filename of the .geo files (will add extension .geo automatically)
......
......@@ -17,6 +17,8 @@ realize this example using the Frackit python bindings.__
* use the `MultiGeometrySampler` and `MultiGeometryEntitySet` classes to conveniently compose a network out of multiple entity geometries
* use the `ConstraintsMatrix` class to conveniently define several constraints to be fulfilled among different entity sets
* define multiple confining/non-confining sub-domains to place the entities in
* register the reasons for entity rejections within the `Status` class and print a summary
* Build and write networks that are confined/unconfined to their embedding domains
### Read in a domain geometry
......
......@@ -59,9 +59,9 @@ int main(int argc, char** argv)
// Compute its volume and get the boundary faces for constraint evaluation.
const auto& networkDomain = solids[1];
const auto domainVolume = computeMagnitude(networkDomain);
const auto shells = OCCUtilities::getShells(networkDomain);
if (shells.size() != 1) throw std::runtime_error("Expected a single shell bounding the domain");
const auto domainBoundaryFaces = OCCUtilities::getFaces(shells[0]);
const auto domainShell = OCCUtilities::getShells(networkDomain);
if (domainShell.size() != 1) throw std::runtime_error("Expected a single shell bounding the domain");
const auto domainBoundaryFaces = OCCUtilities::getFaces(domainShell[0]);
......@@ -157,14 +157,13 @@ int main(int argc, char** argv)
// of entity sets of different geometry types.
MultiGeometryEntitySet<Disk, Quad> entitySets;
// Helper class for terminal output of the creation
// progress and definition of stop criterion etc
// Check if a value has been passed via the command line,
// otherwise use the defaults. This is used in the test suite
// to speed up testing time.
// Helper class for terminal output of the creation progress
// and definition of stop criterion etc. Check if a value has
// been passed via the command line, otherwise use the defaults.
// This is used in the test suite to speed up testing time.
SamplingStatus status;
status.setTargetCount(diskSetId, (argc > 1 ? std::stoi(argv[1]) : 12)); // we want 11 entities of orientation 1
status.setTargetCount(quadSetId, (argc > 1 ? std::stoi(argv[1]) : 16)); // we want 13 entities of orientation 2
status.setTargetCount(diskSetId, (argc > 1 ? std::stoi(argv[1]) : 12)); // desired number of entities of orientation 1
status.setTargetCount(quadSetId, (argc > 1 ? std::stoi(argv[1]) : 16)); // desired number of entities of orientation 2
// The actual network generation loop
ctype containedNetworkArea = 0.0;
......@@ -179,9 +178,9 @@ int main(int argc, char** argv)
if (status.finished(id))
{ status.increaseRejectedCounter("set finished"); continue; }
// Moreover, we want to avoid small fragments (< 250 m²)
// Moreover, we want to avoid small fragments
const auto containedArea = computeContainedMagnitude(geom, networkDomain);
if (containedArea < 350.0)
if (containedArea < 250.0)
{ status.increaseRejectedCounter("minimum contained area violation"); continue; }
// enforce constraints w.r.t. to the other entities
......
......@@ -9,9 +9,8 @@ print("\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
from frackit.occutilities import readShape, getSolids, getShells, getFaces
domainShape = readShape("layers.brep")
# obtain the three solids contained in the file
......@@ -22,12 +21,10 @@ if len(solids) != 3: sys.exit("Expected the .brep file to contain 3 solids")
# Compute its volume and get the boundary faces for constraint evaluation.
from frackit.geometry import computeMagnitude
networkDomain = solids[1]
domainVolume = computeMagnitude(networkDomain);
shells = getShells(networkDomain)
if len(shells) != 1: sys.exit("Expected a single shell bounding the domain")
domainBoundaryFaces = getFaces(shells[0])
if len(domainBoundaryFaces) != 6: sys.exit("Expected 6 faces to bound the domain")
domainVolume = computeMagnitude(networkDomain)
domainShell = getShells(networkDomain)
if len(domainShell) != 1: sys.exit("Expected a single shell bounding the domain")
domainBoundaryFaces = getFaces(domainShell[0])
##############################################################################
......@@ -36,9 +33,8 @@ if len(domainBoundaryFaces) != 6: sys.exit("Expected 6 faces to bound the domain
##############################################################################
# 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, OCCShapeWrapper
domainBBox = getBoundingBox( OCCShapeWrapper(networkDomain) )
from frackit.geometry import OCCShapeWrapper, getBoundingBox
domainBBox = getBoundingBox(networkDomain)
# creates a uniform sampler within an interval
import random
......@@ -51,22 +47,24 @@ def normalSampler(mean, stdDev):
def sample(): return random.gauss(mean, stdDev)
return sample
from frackit.common import toRadians, Id
import math
from frackit.common import 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
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)))
# 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
uniformSampler(30.0, 60.0), # strike length
uniformSampler(30.0, 60.0)) # dip length
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))
# Define ids for the two entity sets
diskSetId = Id(1) # we give the set of orientation one, consisting of disks, the id 1
......@@ -84,7 +82,7 @@ from frackit.entitynetwork import EntityNetworkConstraints, EntityNetworkConstra
def makeDefaultConstraints():
c = EntityNetworkConstraints()
c.setMinDistance(2.5)
c.setMinIntersectingAngle(toRadians(25.0))
c.setMinIntersectingAngle(math.radians(25.0))
c.setMinIntersectionMagnitude(5.0)
c.setMinIntersectionDistance(2.5)
return c
......@@ -100,7 +98,7 @@ constraints2.setMinDistance(5.0)
# Define constraints between entities of different sets
constraintsOnOther = makeDefaultConstraints()
constraintsOnOther.setMinDistance(2.5)
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0))
constraintsOnOther.setMinIntersectingAngle(math.radians(40.0))
# We can use the constraints matrix to facilitate constraint evaluation
constraintsMatrix = EntityNetworkConstraintsMatrix()
......@@ -116,7 +114,7 @@ constraintsMatrix.addConstraints(constraintsOnOther, # constraint instance
# Moreover, we define constraints w.r.t. the domain boundary
constraintsOnDomain = makeDefaultConstraints()
constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0))
constraintsOnDomain.setMinIntersectingAngle(math.radians(15.0))
###########################
......@@ -127,15 +125,15 @@ constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0))
# 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
status.setTargetCount(diskSetId, 12) # number of desired entities of orientation 1
status.setTargetCount(quadSetId, 16) # number of desired entities of orientation 2
# store all entity sets in a dictionary
entitySets = {diskSetId: [], quadSetId: []}
# Alternate between set 1 & set 2 during sampling phase
sampleIntoSet1 = True
containedNetworkArea = 0.0;
containedNetworkArea = 0.0
while not status.finished():
id = diskSetId if sampleIntoSet1 else quadSetId
geom = diskSampler() if sampleIntoSet1 else quadSampler()
......@@ -146,10 +144,10 @@ while not status.finished():
status.increaseRejectedCounter("set finished")
continue
# Moreover, we want to avoid small fragments (< 250 m²)
# Moreover, we want to avoid small fragments
from frackit.geometry import computeContainedMagnitude
containedArea = computeContainedMagnitude(geom, networkDomain);
if containedArea < 350.0:
containedArea = computeContainedMagnitude(geom, networkDomain)
if containedArea < 250.0:
status.increaseRejectedCounter("minimum contained area violation")
continue
......@@ -171,17 +169,17 @@ while not status.finished():
status.print()
# keep track of entity area
containedNetworkArea += containedArea;
containedNetworkArea += containedArea
# sample from the other set next time
sampleIntoSet1 = not sampleIntoSet1
# print the final entity density
density = containedNetworkArea/domainVolume;
density = containedNetworkArea/domainVolume
print("\nEntity density of the contained network: {:f} m²/m³\n".format(density))
# print info on rejection events
status.printRejectionData();
status.printRejectionData()
print("")
##########################################################################
......@@ -192,19 +190,19 @@ print("")
# construct and write a contained network, i.e. write out both network and domain.
print("Building and writing contained, confined network")
from frackit.entitynetwork import EntityNetworkBuilder, ContainedEntityNetworkBuilder
builder = ContainedEntityNetworkBuilder();
builder = ContainedEntityNetworkBuilder()
# add sub-domains
builder.addConfiningSubDomain(solids[0], Id(1));
builder.addConfiningSubDomain(networkDomain, Id(2));
builder.addConfiningSubDomain(solids[2], Id(3));
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 = GmshWriter(builder.build())
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)
......@@ -212,33 +210,33 @@ gmshWriter.write("contained_confined") # body of the filename to be used (will a
# 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")
builder.clear();
builder.addSubDomain(solids[0], Id(1));
builder.addSubDomain(networkDomain, Id(2));
builder.addSubDomain(solids[2], Id(3));
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 = GmshWriter(builder.build())
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.setMeshSize(GmshWriter.GeometryTag.subDomain, 5.0)
gmshWriter.write("contained_unconfined");
gmshWriter.write("contained_unconfined")
# 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));
uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2))
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))
gmshWriter = GmshWriter(uncontainedBuilder.build());
gmshWriter = GmshWriter(uncontainedBuilder.build())
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.write("uncontained_confined");
gmshWriter.write("uncontained_confined")
# ... or not confining it
print("Building and writing uncontained, unconfined network")
uncontainedBuilder.clear();
uncontainedBuilder.clear()
for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))
gmshWriter = GmshWriter(uncontainedBuilder.build());
gmshWriter = GmshWriter(uncontainedBuilder.build())
gmshWriter.setMeshSize(GmshWriter.GeometryTag.entity, 2.5)
gmshWriter.write("uncontained_unconfined");
gmshWriter.write("uncontained_unconfined")
......@@ -14,7 +14,7 @@ realize this example using the Frackit python bindings.__
* create a domain geometry by intersection/cut operations on basic geometry types
* generate random polygons using the `PolygonSampler` class
* design your own incremental network generation algorithm
* design a custom, incremental network generation algorithm
Let us consider some idealized, hexahedral geological layer with an embedded
tunnel structure. The excavation of tunnels can have an effect on the stress
......@@ -24,7 +24,7 @@ of how such situations can modeled using `Frackit`.
### Step 1: domain generation
We want to focus on a spherical region around the center of the hexahedral
We want to focus on a spherical region around the center of a hexahedral
layer. To this end, we construct instances of the `Box` and `Sphere` classes
to represent these two geometries:
......@@ -71,7 +71,7 @@ if (solids.size() != 1) throw std::runtime_error("Cut operation with tunnel shou
const auto& domain = solids[0];
```
## Step 2: entity sampler definitions
### Step 2: entity sampler definitions
Instead of randomly sample entities within the domain volume, in this example we
want to use a different approach. Here, the goal is to generate a network in the
......@@ -196,45 +196,32 @@ rejected immediately.
__Stage 2__: For each entity `e` of stage 1, we want to generate an entity of orientation 2
that intersects `e`. Therefore, we loop over all entities of stage 1 and sample new
candidates in their vicinity. This is achieved by computing the bounding box of `e`,
which is then used to define a region in which new candidates are sampled:
which is then used to define a region in which new candidates are sampled. This
is done within the `getVicinityPointSampler` lambda:
```cpp
for (const auto& e : entitiesSet1)
// this lambda function generates a point sampler within the vicinity of the given geometry
auto getVicinityPointSampler = [&] (const auto& g, unsigned int orientationIdx)
{
using std::sqrt;
// sample an entity of orientation 2 in the vicinity of e
const auto rawBBox = getBoundingBox(e);
const auto charLength = sqrt(computeMagnitude(e)); // characteristic length of e
const auto rawBBox = getBoundingBox(g);
const auto charLength = sqrt(computeMagnitude(g));
const Box sampleBox(rawBBox.xMin()-charLength, rawBBox.yMin()-charLength, rawBBox.zMin()-charLength,
rawBBox.xMax()+charLength, rawBBox.yMax()+charLength, rawBBox.zMax()+charLength);
auto sampler = getSampler(makeUniformPointSampler(sampleBox), /*orientationIdx*/2);
return getSampler(makeUniformPointSampler(sampleBox), orientationIdx);
};
```
The variable `sampler` now holds an instance of the `PolygonSampler` class, which generates
polygons of orientation 2 within the box `sampleBox` that describes the vicinity of `e`.
We then sample candidates until we have found an entity that fulfills all desired constraints.
After that, it is continued with the next entity of `entitiesSet1` until all have been visited.
Reusing the `getSampler` lambda, this returns an instance of the `PolygonSampler` class, which generates
polygons of the given orientation within the box `sampleBox` that describes the vicinity of the given geometry `g`.
Such samplers are then successively created for each entity of stage 1 in order to
find one intersecting entity of orientation 2.
__Stage 3__: this stage follows the same logic as stage 2, but it samples entities of
orientation 3 and looks for entity candidates in the vicinity of the entities of stage 2:
```cpp
for (const auto& e : entitiesSet2)
{
using std::sqrt;
// sample an entity of orientation 2 in the vicinity of e
const auto rawBBox = getBoundingBox(e);
const auto charLength = sqrt(computeMagnitude(e));
const Box sampleBox(rawBBox.xMin()-charLength, rawBBox.yMin()-charLength, rawBBox.zMin()-charLength,
rawBBox.xMax()+charLength, rawBBox.yMax()+charLength, rawBBox.zMax()+charLength);
auto sampler = getSampler(makeUniformPointSampler(sampleBox), /*orientationIdx*/3);
```
Note that due to the restrictive constraints used in this example, the entity
generation step can take up to a few minutes.