example2.py 5.89 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 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 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 ``````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 `````` Dennis Gläser committed May 28, 2020 81 `````` quad = quadSampler1() if sampleIntoSet1 else quadSampler2() `````` Dennis Gläser committed May 27, 2020 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 `````` 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 104 `````` if not constraintsOnBoundary.evaluate(getShape(domain.lateralFace()), quad): `````` Dennis Gläser committed May 27, 2020 105 106 107 108 109 110 111 112 113 114 115 116 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 `````` 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")``````