README.md 13 KB
Newer Older
 Dennis Gläser committed Nov 06, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ``````

Example 4 ========= __As in the previous examples, this description focuses on the C++ implementation given in the main file `example4.cc`, but in `example4.py` you can find how to realize this example using the Frackit python bindings.__ _In this example you will learn how to:_ * create a domain geometry by intersection/cut operations on basic geometry types * generate random polygons using the `PolygonSampler` class `````` Dennis Gläser committed Dec 18, 2020 17 ``````* design a custom, incremental network generation algorithm `````` Dennis Gläser committed Nov 06, 2020 18 19 20 21 22 23 24 25 26 `````` Let us consider some idealized, hexahedral geological layer with an embedded tunnel structure. The excavation of tunnels can have an effect on the stress state in the surrounding rock and can lead to the generation of fractures, in particular in the vicinity of the tunnels. This example aims at an illustration of how such situations can modeled using `Frackit`. ### Step 1: domain generation `````` Dennis Gläser committed Dec 18, 2020 27 ``````We want to focus on a spherical region around the center of a hexahedral `````` Dennis Gläser committed Nov 06, 2020 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 ``````layer. To this end, we construct instances of the `Box` and `Sphere` classes to represent these two geometries: ```cpp // the layer has a thickness of 40m, and initially we consider 100m lateral extent const Box layer(-50, -50, -20.0, 50, 50, 20.0); const Sphere domainSphere(Point(0.0, 0.0, 0.0), /*radius*/25.0); // a tunnel with a radius of 5m crosses the domain const ctype tunnelRadius = 5.0; const Direction tunnelDirection(Vector(1.0, 0.0, 0.0)); const Circle tunnelBase(Point(-25.0, 0.0, 0.0), tunnelDirection, tunnelRadius); const Cylinder tunnel(tunnelBase, /*height*/50.0); ``` To construct the domain in which we want to place the entities of our network, we use the `intersect` function in the `OCCUtilities` namespace to compute the geometry of the volume that is the common part of `domainSphere` and `layer`. The OpenCascade algorithms work on OpenCascade shapes, and therefore, we convert the internal representations into shapes when passing them to the algorithm: ```cpp const auto layerSphereCut = OCCUtilities::intersect(OCCUtilities::getShape(layer), OCCUtilities::getShape(domainSphere), /*eps*/1e-6); auto solids = OCCUtilities::getSolids(layerSphereCut); if (solids.size() != 1) throw std::runtime_error("Intersection operation with domainSphere should yield a single solid"); ``` The result of `intersect` is a generic shape object. However, here we know that the result is a single three-dimensional volume, and thus, we obtain the solid shapes of the intersection result, of which we know there should only be a single one. The function `getSolids` returns a vector of all solids contained in a shape, and here we check if the result contains a single one, as is expected. Subsequently, we cut the tunnel out of the intersection result (with the `cut` function in the `OCCUtilities` namespace) in order to obtain the actual model domain: ```cpp const auto tunnelCut = OCCUtilities::cut(solids[0], OCCUtilities::getShape(tunnel), /*eps*/1e-6); solids = OCCUtilities::getSolids(tunnelCut); if (solids.size() != 1) throw std::runtime_error("Cut operation with tunnel should yield a single solid"); // get the domain shape const auto& domain = solids[0]; ``` `````` Dennis Gläser committed Dec 18, 2020 74 ``````### Step 2: entity sampler definitions `````` Dennis Gläser committed Nov 06, 2020 75 76 77 78 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 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 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 `````` Instead of randomly sample entities within the domain volume, in this example we want to use a different approach. Here, the goal is to generate a network in the vicinity of the tunnel. To this end, we use a three-stage algorithm to construct the network, where in each stage we generate entities of a different orientation. We use polygonal entities, and at this point we define a convenience lambda that, provided a point sampler, returns an instance of `PolygonSampler` of the given orientation: ```cpp // we use the default sampler types, thus, default distributions (see traits classes) // This means, normal distributions for angles and uniform distributions for the size using PolygonSampler = Frackit::PolygonSampler; // lambda function to construct polygon sampler with given point sampler & orientation auto getSampler = [] (const auto& pointSampler, unsigned int orientation) -> PolygonSampler { using NormalDistro = std::normal_distribution; auto sizeDistro = std::uniform_real_distribution(3.0, 6.0); auto numCornerDistro = std::uniform_int_distribution(4, 9); switch (orientation) { case 1: return PolygonSampler(pointSampler, // sampler for polygon centers NormalDistro(toRadians(15.0), toRadians(10.0)), // strike angle: mean value & standard deviation NormalDistro(toRadians(90.0), toRadians(10.0)), // dip angle: mean value & standard deviation /*strikeLength/*/sizeDistro, /*dipLength*/sizeDistro, numCornerDistro); case 2: return PolygonSampler(pointSampler, // sampler for polygon centers NormalDistro(toRadians(0.0), toRadians(10.0)), // strike angle: mean value & standard deviation NormalDistro(toRadians(0.0), toRadians(10.0)), // dip angle: mean value & standard deviation /*strikeLength/*/sizeDistro, /*dipLength*/sizeDistro, numCornerDistro); case 3: return PolygonSampler(pointSampler, // sampler for polygon centers NormalDistro(toRadians(-75.0), toRadians(10.0)), // strike angle: mean value & standard deviation NormalDistro(toRadians(90.0), toRadians(10.0)), // dip angle: mean value & standard deviation /*strikeLength/*/sizeDistro, /*dipLength*/sizeDistro, numCornerDistro); default: throw std::runtime_error("Unsupported orientation index"); } }; ``` As can be seen, an instance of the `PolygonSampler` class is constructed with statistical distributions on the strike and dip angles (polygon orientations), the size of the polygon in strike and dip direction (polygon sizes) and the number of corners. We use the default distribution types here, which means normal distributions for the orientation angles and uniform distributions for the sizes and number of corners. As for the other samplers, one can choose different distributions by implementing and passing a traits class as template argument to `PolygonSampler`. For more details we refer to the corrensponding [class documentation](https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/frackit/sampling/polygonsampler.hh) ### Step 3: constraints definitions The first three examples already explained how to define and use the `Constraints` class of `Frackit`. As in those examples, here we also use several instances of that class to define constraints to be fulfilled between entities of different orientations, or between entities and the domain boundary. For more details on this, please look at the previous examples. Here, we want to focus on additional custom constraints that we would like to be fulfilled. That is, we want to ensure that the intersections between different entities do not produce very small length scales on the polygon boundaries as a result of being very close to a polygon corner. To this end, we define the following lambda function to be used during network generation: ```cpp // When polygons intersect, the intersection might be very close to // one of the corner points. We want to avoid small length scales caused // by this, and this lambda provides the check for it. auto producesSmallLengthScale = [] (const auto& geometry, const auto& is) -> bool { // lambda to check if the intersection is too close (< 1mm) to a vertex of the geometry auto isTooClose = [&is] (const auto& v) { return computeDistance(v, is) < 1e-3; }; const auto shape = OCCUtilities::getShape(geometry); const auto vertices = OCCUtilities::getVertices(shape); return std::any_of(vertices.begin(), vertices.end(), isTooClose); }; ``` The lambda `producesSmallLengthScale` checks if the intersection `is` is less than 1 mm away from any of the vertices of the given geometry. Note that in this lambda we again operate on shapes rather than the internal representation of polygons. The reason for this is that the above implementation is more generic and works for any given geometry type. In this example, it is called both with polygons as well as OpenCascade shapes. ### Step 4: entity sampling As mentioned before, here we want to construct the network in three stages: * Stage 1: generate entities of orientation 1 close to the tunnel, each entity must intersect the tunnel * Stage 2: for each entity `e1` of stage 1, generate an entity of orientation 2 that intersects `e1` * Stage 3: for each entity `e2` of stage 2, generate an entity of orientation 3 that intersects `e2` This ensures that all entities of the network are connected to at least one other, and that the overall network remains in the vicinity of the tunnel. __Stage 1__: In order to efficiently sample entities that intersect the tunnel, we search for new candidates in a small region around it. Recall that our lambda function `getSampler` expects a point sampler to be passed, which is used to sample the spatial location of a new polygonal entity candidate. In this first stage, we use a point sampler that uniformly samples points within a hollow cylinder around the tunnel. This is done in following two lines of code. ```cpp const HollowCylinder sampleCylinder1(tunnel.bottomFace().center(), tunnel.direction(), tunnel.radius(), // inner radius tunnel.radius() + 2.0, // outer radius tunnel.height()); auto polygonSampler1 = getSampler(makeUniformPointSampler(sampleCylinder1), /*orientationIdx*/1); ``` Subsequently, `polygonSampler1` is used to get entity candidates until the desired number of entities has been accepted, checking each entity against a number of constraints. In particular, entities that do not intersect the tunnel are rejected immediately. __Stage 2__: For each entity `e` of stage 1, we want to generate an entity of orientation 2 that intersects `e`. Therefore, we loop over all entities of stage 1 and sample new candidates in their vicinity. This is achieved by computing the bounding box of `e`, `````` Dennis Gläser committed Dec 18, 2020 199 200 ``````which is then used to define a region in which new candidates are sampled. This is done within the `getVicinityPointSampler` lambda: `````` Dennis Gläser committed Nov 06, 2020 201 202 `````` ```cpp `````` Dennis Gläser committed Dec 18, 2020 203 204 ``````// this lambda function generates a point sampler within the vicinity of the given geometry auto getVicinityPointSampler = [&] (const auto& g, unsigned int orientationIdx) `````` Dennis Gläser committed Nov 06, 2020 205 206 ``````{ using std::sqrt; `````` Dennis Gläser committed Dec 18, 2020 207 208 `````` const auto rawBBox = getBoundingBox(g); const auto charLength = sqrt(computeMagnitude(g)); `````` Dennis Gläser committed Nov 06, 2020 209 210 211 `````` const Box sampleBox(rawBBox.xMin()-charLength, rawBBox.yMin()-charLength, rawBBox.zMin()-charLength, rawBBox.xMax()+charLength, rawBBox.yMax()+charLength, rawBBox.zMax()+charLength); `````` Dennis Gläser committed Dec 18, 2020 212 213 `````` return getSampler(makeUniformPointSampler(sampleBox), orientationIdx); }; `````` Dennis Gläser committed Nov 06, 2020 214 215 ````````` `````` Dennis Gläser committed Dec 18, 2020 216 217 218 219 ``````Reusing the `getSampler` lambda, this returns an instance of the `PolygonSampler` class, which generates polygons of the given orientation within the box `sampleBox` that describes the vicinity of the given geometry `g`. Such samplers are then successively created for each entity of stage 1 in order to find one intersecting entity of orientation 2. `````` Dennis Gläser committed Nov 06, 2020 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 `````` __Stage 3__: this stage follows the same logic as stage 2, but it samples entities of orientation 3 and looks for entity candidates in the vicinity of the entities of stage 2: Note that due to the restrictive constraints used in this example, the entity generation step can take up to a few minutes. ### Step 5: network generation As in the previous example, we use the `ContainedEntityNetworkBuilder` class to construct the network and embed it in our domain: ```cpp // construct and write a contained network, i.e. write out both network and domain. std::cout << "Building and writing network" << std::endl; ContainedEntityNetworkBuilder builder; `````` Dennis Gläser committed Dec 18, 2020 238 239 240 241 ``````// add sub-domain builder.addConfiningSubDomain(domain, Id(1)); // embed all entities in that domain `````` Dennis Gläser committed Nov 06, 2020 242 243 244 245 246 247 248 249 250 251 252 253 ``````builder.addSubDomainEntities(entitiesSet1, Id(1)); builder.addSubDomainEntities(entitiesSet2, Id(1)); builder.addSubDomainEntities(entitiesSet3, Id(1)); ``` Here, we give the domain the `Id` "1" and define all entities to be embedded within it. The following lines of code then perform the actual construction of the network in the call `builer.build()`, and write out both the network and the domain to Gmsh file format: ```cpp // now we can build and write out the network in Gmsh file format `````` Dennis Gläser committed Dec 10, 2020 254 255 256 ``````gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 0.5); gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0); gmshWriter.write("network"); `````` Dennis Gläser committed Nov 06, 2020 257 ```````