Skip to content
GitLab
Commit 2b31244a authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[ex3] modify example code

parent 1bb24325
Hide whitespace changes
Inline Side-by-side
appl/example3/README.md
View file @ 2b31244a
......@@ -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
......
appl/example3/example3.cc
View file @ 2b31244a
......@@ -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
......
appl/example3/example3.py
View file @ 2b31244a
......@@ -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")
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment