example2.py 6.1 KB
 Dennis Gläser committed May 27, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ``````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 `````` Dennis Gläser committed Oct 27, 2020 19 20 21 22 23 24 25 ``````#returns a sampler from a uniform distribution between min and max def uniformSampler(min, max): import random def sample(): return random.uniform(min, max) return sample `````` Dennis Gläser committed May 27, 2020 26 27 28 29 30 31 ``````# 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 `````` Dennis Gläser committed Oct 27, 2020 32 33 `````` uniformSampler(0.4, 0.6), # strike length uniformSampler(0.4, 0.6)) # dip length `````` Dennis Gläser committed May 27, 2020 34 35 36 37 38 `````` # 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 `````` Dennis Gläser committed Oct 27, 2020 39 40 `````` uniformSampler(0.4, 0.6), # strike length uniformSampler(0.4, 0.6)) # dip length `````` Dennis Gläser committed May 27, 2020 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 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 85 86 87 `````` # 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 `````` Dennis Gläser committed May 28, 2020 88 `````` quad = quadSampler1() if sampleIntoSet1 else quadSampler2() `````` Dennis Gläser committed May 27, 2020 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 `````` 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 111 `````` if not constraintsOnBoundary.evaluate(getShape(domain.lateralFace()), quad): `````` Dennis Gläser committed May 27, 2020 112 113 114 115 `````` status.increaseRejectedCounter() continue # reject entities whose contained part is too small `````` Dennis Gläser committed Nov 05, 2020 116 `````` from frackit.geometry import computeContainedMagnitude `````` Dennis Gläser committed May 27, 2020 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 `````` 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); `````` Dennis Gläser committed Dec 10, 2020 151 ``````writer.setMeshSize(GmshWriter.GeometryTag.entity, 0.1) `````` Dennis Gläser committed Dec 10, 2020 152 ``````writer.setMeshSize(GmshWriter.GeometryTag.subDomain, 0.1) `````` Dennis Gläser committed Dec 10, 2020 153 ``````writer.write("network") # filename of the .geo files (will add extension .geo automatically) `````` Dennis Gläser committed May 27, 2020 154 ``print("\n --- Finished writing .geo file ---\n")``