example2.py 6.07 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 `````` `````` Dennis Gläser committed Dec 18, 2020 5 ``````domain = geometry.Cylinder(radius=0.5, height=1.0) `````` Dennis Gläser committed May 27, 2020 6 7 8 9 `````` # 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 `````` Dennis Gläser committed Dec 18, 2020 10 11 ``````box = geometry.Box(xmin=-0.5, ymin=-0.5, zmin=0.0, xmax=0.5, ymax=0.5, zmax=1.0) `````` Dennis Gläser committed May 27, 2020 12 13 14 15 ``````pointSampler = makeUniformPointSampler(box) # returns a sampler from a gaussian distribution with mean and std deviation def gaussianSampler(mean, stdDev): `````` Dennis Gläser committed Dec 18, 2020 16 `````` def sample(): return random.gauss(mean, stdDev) `````` Dennis Gläser committed May 27, 2020 17 18 `````` return sample `````` Dennis Gläser committed Oct 27, 2020 19 20 ``````#returns a sampler from a uniform distribution between min and max def uniformSampler(min, max): `````` Dennis Gläser committed Dec 18, 2020 21 `````` def sample(): return random.uniform(min, max) `````` Dennis Gläser committed Oct 27, 2020 22 23 `````` return sample `````` Dennis Gläser committed May 27, 2020 24 25 ``````# we sample quadrialeterals within the box with gaussian distributions for the parameters from frackit.sampling import QuadrilateralSampler as QuadSampler `````` Dennis Gläser committed Dec 18, 2020 26 27 28 29 30 ``````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)) `````` Dennis Gläser committed May 27, 2020 31 32 `````` # sampler for quadrilaterals of the secondary orientation `````` Dennis Gläser committed Dec 18, 2020 33 34 35 36 37 ``````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)) `````` Dennis Gläser committed May 27, 2020 38 39 40 41 42 43 `````` # constructs a constraints object with default settings for this example from frackit.entitynetwork import EntityNetworkConstraints def makeDefaultConstraints(): c = EntityNetworkConstraints() c.setMinDistance(0.05) `````` Dennis Gläser committed Dec 18, 2020 44 `````` c.setMinIntersectingAngle(math.radians(30.0)) `````` Dennis Gläser committed May 27, 2020 45 46 47 48 49 50 51 52 53 54 55 56 57 58 `````` 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() `````` Dennis Gläser committed Dec 18, 2020 59 ``````constraintsOnOther.setMinIntersectingAngle(math.radians(40.0)) `````` Dennis Gläser committed May 27, 2020 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 `````` # 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 `````` Dennis Gläser committed May 28, 2020 85 `````` quad = quadSampler1() if sampleIntoSet1 else quadSampler2() `````` Dennis Gläser committed May 27, 2020 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 `````` 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 `````` Dennis Gläser committed Oct 07, 2020 108 `````` if not constraintsOnBoundary.evaluate(getShape(domain.lateralFace()), quad): `````` Dennis Gläser committed May 27, 2020 109 110 111 112 `````` status.increaseRejectedCounter() continue # reject entities whose contained part is too small `````` Dennis Gläser committed Nov 05, 2020 113 `````` from frackit.geometry import computeContainedMagnitude `````` Dennis Gläser committed Dec 18, 2020 114 `````` containedArea = computeContainedMagnitude(quad, domain) `````` Dennis Gläser committed May 27, 2020 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 `````` 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 `````` Dennis Gläser committed Dec 18, 2020 135 ``````builder.addConfiningSubDomain(domain, Id(1)) `````` Dennis Gläser committed May 27, 2020 136 137 `````` # define all entities to be embedded in this domain `````` Dennis Gläser committed Dec 18, 2020 138 139 ``````builder.addSubDomainEntities(entities1, Id(1)) builder.addSubDomainEntities(entities2, Id(1)) `````` Dennis Gläser committed May 27, 2020 140 141 142 `````` # 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 143 ``````network = builder.build() `````` Dennis Gläser committed May 27, 2020 144 145 146 `````` print("\n --- Writing .geo file ---\n") from frackit.io import GmshWriter `````` Dennis Gläser committed Dec 18, 2020 147 ``````writer = GmshWriter(network) `````` Dennis Gläser committed Dec 10, 2020 148 ``````writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1) `````` Dennis Gläser committed Dec 10, 2020 149 ``````writer.setMeshSize(GmshWriter.GeometryTag.subDomain, 0.1) `````` Dennis Gläser committed Dec 10, 2020 150 ``````writer.write("network") # filename of the .geo files (will add extension .geo automatically) `````` Dennis Gläser committed May 27, 2020 151 ``print("\n --- Finished writing .geo file ---\n")``