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

Feature/python examples

See merge request tools/frackit!97

appl/example1/example1.py

+import frackit.geometry as geometry
+
+# we use the unit cube as domain
+box = geometry.Box(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
+
+# we sample points uniformly within the domain
+from frackit.sampling import makeUniformPointSampler
+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)
+    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
+                           gaussianSampler(0.5, 0.1),                        # edge length
+                           0.05)                                             # threshold for minimum edge length
+
+# 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
+                           gaussianSampler(0.5, 0.1),                       # edge length
+                           0.05)                                            # threshold for minimum edge length
+
+# We want to enforce some constraints on the set of quadrilaterals.
+# In particular, for entities of the same set we want a minimum spacing
+# distance of 5cm, and the quadrilaterals must not intersect in angles
+# less than 30°. Moreover, if they intersect, we don't want intersection
+# edges whose length is smaller than 5cm, and, the intersection should not
+# be too close to the boundary of one of two intersecting quadrilaterals. Here: 5cm.
+from frackit.entitynetwork import EntityNetworkConstraints
+constraintsOnSelf = EntityNetworkConstraints()
+constraintsOnSelf.setMinDistance(0.05)
+constraintsOnSelf.setMinIntersectingAngle(toRadians(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.setMinIntersectionMagnitude(0.05)
+constraintsOnOther.setMinIntersectionDistance(0.05)
+
+# we use unique identifiers for the entities of the two orientations
+from frackit.common import Id
+idSet1 = Id(1)
+idSet2 = Id(2)
+
+# use the status class to define when to stop sampling
+from frackit.sampling import SamplingStatus
+status = SamplingStatus()
+status.setTargetCount(idSet1, 15) # we want 15 entities in each set
+status.setTargetCount(idSet2, 15) # we want 15 entities in each set
+
+# start sampling into set 1 and keep alternating
+sampleIntoSet1 = True
+
+# lists to store the sampled entities
+entities1 = []
+entities2 = []
+
+print("\n --- Start entity sampling ---\n")
+while not status.finished():
+    # sample a quadrilateral, alternating between sampler 1 and sampler 2
+    quad = quadSampler1.sample() if sampleIntoSet1 else quadSampler2.sample()
+    entitySet = entities1 if sampleIntoSet1 else entities2
+    otherEntitySet = entities2 if sampleIntoSet1 else entities1
+
+    # skip the rest if constraints are violated
+    if not constraintsOnSelf.evaluate(entitySet, quad):
+        status.increaseRejectedCounter()
+        continue
+
+    # skip the rest if constraints are violated
+    if not constraintsOnOther.evaluate(otherEntitySet, quad):
+        status.increaseRejectedCounter()
+        continue
+
+    # the entity is admissible
+    entitySet.append(quad)
+    id = idSet1 if sampleIntoSet1 else idSet2
+    status.increaseCounter(id)
+    status.print()
+
+    # sample into other set the next time
+    sampleIntoSet1 = not sampleIntoSet1
+
+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);
+
+# 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();
+
+print("\n --- Writing .geo file ---\n")
+from frackit.io import GmshWriter
+writer = GmshWriter(network);
+writer.write("network", # filename of the .geo files (will add extension .geo automatically)
+             0.1);      # element size to be used
+print("\n --- Finished writing .geo file ---\n")

appl/example2/example2.py

+import frackit.geometry as geometry
+
+# we want to use a cylindrical domain with radius=0.5 and height=1.0
+domain = geometry.Cylinder(0.5, 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)
+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)
+    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
+                           gaussianSampler(0.5, 0.1),                        # edge length
+                           0.05)                                             # threshold for minimum edge length
+
+# 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
+                           gaussianSampler(0.5, 0.1),                       # edge length
+                           0.05)                                            # threshold for minimum edge length
+
+# 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.setMinIntersectionMagnitude(0.05)
+    c.setMinIntersectionDistance(0.05)
+    return c
+
+# We want to enforce some constraints on the set of quadrilaterals.
+# In particular, for entities of the same set we want a minimum spacing
+# distance of 5cm, and the quadrilaterals must not intersect in angles
+# less than 30°. Moreover, if they intersect, we don't want intersection
+# edges whose length is smaller than 5cm, and, the intersection should not
+# be too close to the boundary of one of two intersecting quadrilaterals. Here: 5cm.
+constraintsOnSelf = makeDefaultConstraints()
+
+# with respect to entities of the other set, we want to have larger intersection angles
+constraintsOnOther = makeDefaultConstraints()
+constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
+
+# we use the default constraints w.r.t. to the domain boundary
+constraintsOnBoundary = makeDefaultConstraints()
+
+# we use unique identifiers for the entities of the two orientations
+from frackit.common import Id
+idSet1 = Id(1)
+idSet2 = Id(2)
+
+# use the status class to define when to stop sampling
+from frackit.sampling import SamplingStatus
+status = SamplingStatus()
+status.setTargetCount(idSet1, 15) # we want 15 entities in each set
+status.setTargetCount(idSet2, 15) # we want 15 entities in each set
+
+# start sampling into set 1 and keep alternating
+sampleIntoSet1 = True
+
+# lists to store the sampled entities
+entities1 = []
+entities2 = []
+
+print("\n --- Start entity sampling ---\n")
+while not status.finished():
+    # sample a quadrilateral, alternating between sampler 1 and sampler 2
+    quad = quadSampler1.sample() if sampleIntoSet1 else quadSampler2.sample()
+    entitySet = entities1 if sampleIntoSet1 else entities2
+    otherEntitySet = entities2 if sampleIntoSet1 else entities1
+
+    # skip the rest if constraints are violated
+    if not constraintsOnSelf.evaluate(entitySet, quad):
+        status.increaseRejectedCounter()
+        continue
+
+    # skip the rest if constraints are violated
+    if not constraintsOnOther.evaluate(otherEntitySet, quad):
+        status.increaseRejectedCounter()
+        continue
+
+    # enforce constraints w.r.t. to the domain boundary
+    if not constraintsOnBoundary.evaluate(domain.topFace(), quad):
+        status.increaseRejectedCounter()
+        continue
+    if not constraintsOnBoundary.evaluate(domain.bottomFace(), quad):
+        status.increaseRejectedCounter()
+        continue
+
+    from frackit.occutilities import getShape
+    if not constraintsOnBoundary.evaluate(getShape(domain.bottomFace()), quad):
+        status.increaseRejectedCounter()
+        continue
+
+    # reject entities whose contained part is too small
+    from frackit.magnitude import computeContainedMagnitude
+    containedArea = computeContainedMagnitude(quad, domain);
+    if (containedArea < 0.01):
+        status.increaseRejectedCounter()
+        continue
+
+    # the entity is admissible
+    entitySet.append(quad)
+    id = idSet1 if sampleIntoSet1 else idSet2
+    status.increaseCounter(id)
+    status.print()
+
+    # sample into other set the next time
+    sampleIntoSet1 = not sampleIntoSet1
+
+print("\n --- Finished entity sampling ---\n")
+
+# We can now create an entity network from the two sets
+from frackit.entitynetwork import ContainedEntityNetworkBuilder
+builder = ContainedEntityNetworkBuilder()
+
+# our domain confines all entities and receives the 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));
+
+# 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();
+
+print("\n --- Writing .geo file ---\n")
+from frackit.io import GmshWriter
+writer = GmshWriter(network);
+writer.write("network", # filename of the .geo files (will add extension .geo automatically)
+             0.1);      # element size to be used
+print("\n --- Finished writing .geo file ---\n")

appl/example3/example3.cc

@@ -188,7 +188,7 @@ int main(int argc, char** argv)
         if (!constraintsOnDomain.evaluate(domainBoundaryFaces, geom))
         { status.increaseRejectedCounter(); continue; }
 
-        // the disk is admissible
+        // the geometry is admissible
         entitySets.addEntity(geom, id);
         status.increaseCounter(id);
         status.print();

appl/example3/example3.py

+# 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)
+if len(solids) != 3: raise RuntimeError("Expected the .brep file to contain 3 solids")
+
+# 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.
+from frackit.magnitude import computeMagnitude
+networkDomain = solids[1]
+domainVolume = computeMagnitude(networkDomain);
+shells = getShells(networkDomain)
+if len(shells) != 1: raise RuntimeError("Expected a single shell bounding the domain")
+
+domainBoundaryFaces = getFaces(shells[0])
+if len(domainBoundaryFaces) != 6: raise RuntimeError("Expected 6 faces to bound the domain")
+
+
+##############################################################################
+## 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
+from frackit.geometry import Box
+from frackit.occutilities import OCCShapeWrapper
+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
+                                   normalSampler(45.0, 5.0),                        # edge length: mean value & standard deviation
+                                   5.0)                                             # threshold for minimum edge length
+
+# 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()
+constraintsMatrix.addConstraints(constraints1,                       # constraint instance
+                                 [diskSetId.get(), diskSetId.get()]) # sets between which to use these constraints
+
+constraintsMatrix.addConstraints(constraints2,                       # constraint instance
+                                 [quadSetId.get(), quadSetId.get()]) # sets between which to use these constraints
+
+constraintsMatrix.addConstraints(constraintsOnOther,                   # constraint instance
+                                 [[diskSetId.get(), quadSetId.get()],  # sets between which to use these constraints
+                                  [quadSetId.get(), diskSetId.get()]]) # sets between which to use these constraints
+
+# 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
+entitySets = {diskSetId.get(): [], quadSetId.get(): []}
+
+# Alternate between set 1 & set 2 during sampling phase
+sampleIntoSet1 = True
+containedNetworkArea = 0.0;
+while not status.finished():
+    id = diskSetId if sampleIntoSet1 else quadSetId
+    geom = diskSampler.sample() if sampleIntoSet1 else quadSampler.sample()
+
+    # 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²)
+    from frackit.magnitude import computeContainedMagnitude
+    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
+    entitySets[id.get()].append(geom)
+    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.
+print("Building and writing contained, confined network\n")
+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
+print("Building and writing contained, unconfined network\n")
+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
+print("Building and writing uncontained, unconfined network\n")
+uncontainedBuilder.clear();
+for setId in entitySets: uncontainedBuilder.addSubDomainEntities(entitySets[setId], Id(2))
+
+gmshWriter = GmshWriter(uncontainedBuilder.build());
+gmshWriter.write("uncontained_unconfined", 2.5);