example1.py 4.89 KB
 Dennis Gläser committed May 27, 2020 1 ``````import frackit.geometry as geometry `````` Dennis Gläser committed Dec 18, 2020 2 3 ``````import random import math `````` Dennis Gläser committed May 27, 2020 4 5 `````` # we use the unit cube as domain `````` Dennis Gläser committed Dec 18, 2020 6 7 ``````box = geometry.Box(xmin=0.0, ymin=0.0, zmin=0.0, xmax=1.0, ymax=1.0, zmax=1.0) `````` Dennis Gläser committed May 27, 2020 8 9 10 11 12 13 `````` # 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): `````` Dennis Gläser committed Dec 18, 2020 14 `````` def sample(): return random.gauss(mean, stdDev) `````` Dennis Gläser committed May 27, 2020 15 16 `````` return sample `````` Dennis Gläser committed Oct 27, 2020 17 18 ``````# returns a sampler from a uniform distribution between min and max def uniformSampler(min, max): `````` Dennis Gläser committed Dec 18, 2020 19 `````` def sample(): return random.uniform(min, max) `````` Dennis Gläser committed Oct 27, 2020 20 21 `````` return sample `````` Dennis Gläser committed Dec 18, 2020 22 ``````# we sample quadrilaterals within the box with gaussian distributions for the parameters `````` Dennis Gläser committed May 27, 2020 23 ``````from frackit.sampling import QuadrilateralSampler as QuadSampler `````` Dennis Gläser committed Dec 18, 2020 24 25 26 27 28 ``````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)) `````` Dennis Gläser committed May 27, 2020 29 30 `````` # sampler for quadrilaterals of the secondary orientation `````` Dennis Gläser committed Dec 18, 2020 31 32 33 34 35 ``````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)) `````` Dennis Gläser committed May 27, 2020 36 37 38 39 40 41 42 43 44 45 `````` # 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) `````` Dennis Gläser committed Dec 18, 2020 46 ``````constraintsOnSelf.setMinIntersectingAngle(math.radians(30.0)) `````` Dennis Gläser committed May 27, 2020 47 48 49 50 51 52 ``````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) `````` Dennis Gläser committed Dec 18, 2020 53 ``````constraintsOnOther.setMinIntersectingAngle(math.radians(40.0)) `````` Dennis Gläser committed May 27, 2020 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 ``````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 `````` Dennis Gläser committed May 28, 2020 78 `````` quad = quadSampler1() if sampleIntoSet1 else quadSampler2() `````` Dennis Gläser committed May 27, 2020 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 `````` 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() `````` Dennis Gläser committed Dec 18, 2020 106 107 ``````builder.addEntities(entities1) builder.addEntities(entities2) `````` Dennis Gläser committed May 27, 2020 108 109 110 `````` # let the builder construct the network and write it to gmsh file format print("\n --- Constructing entity network from the raw entities ---\n") `````` Dennis Gläser committed Dec 18, 2020 111 ``````network = builder.build() `````` Dennis Gläser committed May 27, 2020 112 113 114 `````` print("\n --- Writing .geo file ---\n") from frackit.io import GmshWriter `````` Dennis Gläser committed Dec 18, 2020 115 ``````writer = GmshWriter(network) `````` Dennis Gläser committed Dec 10, 2020 116 117 118 ``````writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1) writer.write("network") # filename of the .geo files (will add extension .geo automatically) `````` Dennis Gläser committed May 27, 2020 119 ``print("\n --- Finished writing .geo file ---\n")``