example3.cc 12.4 KB
Newer Older
Dennis Gläser's avatar
Dennis Gläser committed
1
2
3
4
5
6
7
#include <iostream>
#include <stdexcept>
#include <random>

// basic geometry types
#include <frackit/geometry/box.hh>
#include <frackit/geometry/disk.hh>
8
#include <frackit/geometry/quadrilateral.hh>
Dennis Gläser's avatar
Dennis Gläser committed
9
10
11
12
13
14

// headers containing functions to compute lengths/areas/volumes
#include <frackit/magnitude/magnitude.hh>
#include <frackit/magnitude/containedmagnitude.hh>

// sampler for points and disks
15
#include <frackit/sampling/makeuniformpointsampler.hh>
16
#include <frackit/sampling/disksampler.hh>
17
18
19
#include <frackit/sampling/quadrilateralsampler.hh>
#include <frackit/sampling/multigeometrysampler.hh>
#include <frackit/sampling/status.hh>
Dennis Gläser's avatar
Dennis Gläser committed
20
21
22

// constraints to be enforced on the network (distance, angles, etc.)
#include <frackit/entitynetwork/constraints.hh>
23
24
#include <frackit/entitynetwork/constraintsmatrix.hh>
#include <frackit/entitynetwork/multigeometryentityset.hh>
Dennis Gläser's avatar
Dennis Gläser committed
25
26
27
28
29
30
31
32
33
34
35

// builder class for creating networks of entities confined in (sub-)domains
#include <frackit/entitynetwork/networkbuilder.hh>

// writes an entity network to a meshable Gmsh .geo file format
#include <frackit/io/gmshwriter.hh>

int main(int argc, char** argv)
{
    //! print welcome message
    std::cout << "\n\n"
36
37
38
              << "##############################################################################\n"
              << "## Example: Creation of entities (disks/quadrilaterals) in a layered medium ##\n"
              << "##############################################################################\n"
Dennis Gläser's avatar
Dennis Gläser committed
39
40
              << "\n\n";

41
42
43
44
45
46
    // Define some types used here
    using namespace Frackit;
    using ctype = double;
    using Disk = Disk<ctype>;
    using Quad = Quadrilateral<ctype, 3>;

Dennis Gläser's avatar
Dennis Gläser committed
47
48
49
50
    /////////////////////////////////////////////////////
    // 1. Read in the domain geometry from .brep file. //
    //    The file name is defined in CMakeLists.txt   //
    /////////////////////////////////////////////////////
51
    const auto domainShape = OCCUtilities::readShape(BREPFILE);
Dennis Gläser's avatar
Dennis Gläser committed
52
53
54
55

    // obtain the three solids contained in the file
    const auto solids = OCCUtilities::getSolids(domainShape);
    if (solids.size() != 3)
56
        throw std::runtime_error("Expected the .brep file to contain 3 solids");
Dennis Gläser's avatar
Dennis Gläser committed
57

58
59
    // The sub-domain we want to create a network in is the center one.
    // Compute its volume and get the boundary faces for constraint evaluation.
Dennis Gläser's avatar
Dennis Gläser committed
60
61
    const auto& networkDomain = solids[1];
    const auto domainVolume = computeMagnitude(networkDomain);
62
63
64
    const auto shells = OCCUtilities::getShells(networkDomain);
    if (shells.size() != 1) throw std::runtime_error("Expected a single shell bounding the domain");
    const auto domainBoundaryFaces = OCCUtilities::getFaces(shells[0]);
Dennis Gläser's avatar
Dennis Gläser committed
65

66
67
68



69
70
71
72
    //////////////////////////////////////////////////////////////////////////////
    // 2. Make sampler classes to randomly sample points (entity centers)       //
    //    and disk/quadrilaterals as entities (using the sampled center points) //
    //////////////////////////////////////////////////////////////////////////////
Dennis Gläser's avatar
Dennis Gläser committed
73

Dennis Gläser's avatar
Dennis Gläser committed
74
75
76
    // we use the default sampler types, thus, default distributions (see traits classes)
    using NormalDistro = std::normal_distribution<ctype>;
    using UniformDistro = std::uniform_real_distribution<ctype>;
Dennis Gläser's avatar
Dennis Gläser committed
77

78
    // Bounding box of the domain in which we want to place the entities
Dennis Gläser's avatar
Dennis Gläser committed
79
80
    const auto domainBBox = OCCUtilities::getBoundingBox(networkDomain);

81
82
    // sampler for disks (orientation 1)
    DiskSampler diskSampler(makeUniformPointSampler(domainBBox),                               // sampler for disk center points
Dennis Gläser's avatar
Dennis Gläser committed
83
84
85
86
87
                            NormalDistro(30.0, 6.5),                        // major axis length: mean value & standard deviation
                            NormalDistro(24.0, 4.5),                        // minor axis length: mean value & standard deviation
                            NormalDistro(toRadians(0.0), toRadians(7.5)),   // rotation around x-axis: mean value & standard deviation
                            NormalDistro(toRadians(0.0), toRadians(7.5)),   // rotation around y-axis: mean value & standard deviation
                            NormalDistro(toRadians(0.0), toRadians(7.5)));  // rotation around z-axis: mean value & standard deviation
88
89
90

    // sampler for quadrilaterals (orientation 2)
    QuadrilateralSampler<3> quadSampler(makeUniformPointSampler(domainBBox),                               // sampler for quadrilateral center points
Dennis Gläser's avatar
Dennis Gläser committed
91
92
93
94
                                        NormalDistro(toRadians(45.0), toRadians(5.0)),  // strike angle: mean value & standard deviation
                                        NormalDistro(toRadians(90.0), toRadians(5.0)),  // dip angle: mean value & standard deviation
                                        UniformDistro(30.0, 60.0),                      // strike length
                                        UniformDistro(30.0, 60.0));                     // dip length
95
96
97
98
99
100
101
102
103
104
105
106

    // Define ids for the two entity sets
    const Id diskSetId(1); // we give the set of orientation one, consisting of disks, the id 1
    const Id quadSetId(2); // we give the set of orientation two, consisting of quadrilaterals, the id 2

    // Sampler that samples both geometry types (disks & quads)
    // In this one can define an arbitrary number of samplers, each
    // of which is associated with an entity set with a unique id
    MultiGeometrySampler<Disk, Quad> multiSampler;
    multiSampler.addGeometrySampler(diskSampler, diskSetId);
    multiSampler.addGeometrySampler(quadSampler, quadSetId);

107
108
109
110
111
112
113
114



    ///////////////////////////////////////////////////////////////////////
    // 3. Define constraints that should be fulfilled among the entities //
    //    of different orientations.                                     //
    ///////////////////////////////////////////////////////////////////////

115
    // Define constraints between entities of orientation 1
116
    using Constraints = EntityNetworkConstraints<ctype>;
117
    Constraints constraints1;
Dennis Gläser's avatar
Dennis Gläser committed
118
    constraints1.setMinDistance(2.5);
119
    constraints1.setMinIntersectingAngle(toRadians(25.0));
Dennis Gläser's avatar
Dennis Gläser committed
120
121
    constraints1.setMinIntersectionMagnitude(5.0);
    constraints1.setMinIntersectionDistance(2.5);
122
123

    // Define constraints between entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
124
125
126
    // we want to enforce larger spacing between those entities
    auto constraints2 = constraints1;
    constraints2.setMinDistance(5.0);
Dennis Gläser's avatar
Dennis Gläser committed
127

128
    // Define constraints between entities of different sets
Dennis Gläser's avatar
Dennis Gläser committed
129
    auto constraintsOnOther = constraints1;
130
131
    constraintsOnOther.setMinDistance(2.5);
    constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
Dennis Gläser's avatar
Dennis Gläser committed
132

Dennis Gläser's avatar
Dennis Gläser committed
133
    // We can use the constraints matrix to facilitate constraint evaluation
134
    EntityNetworkConstraintsMatrix<Constraints> constraintsMatrix;
135
136
137
138
139
    constraintsMatrix.addConstraints(constraints1,                  // constraint instance
                                     IdPair(diskSetId, diskSetId)); // set between which to use these constraints

    constraintsMatrix.addConstraints(constraints2,                  // constraint instance
                                     IdPair(quadSetId, quadSetId)); // set between which to use these constraints
140
141
142
143
144
145

    constraintsMatrix.addConstraints(constraintsOnOther,              // constraint instance
                                     {IdPair(diskSetId, quadSetId),   // sets between which to use these constraints
                                      IdPair(quadSetId, diskSetId)}); // sets between which to use these constraints

    // Moreover, we define constraints w.r.t. the domain boundary
Dennis Gläser's avatar
Dennis Gläser committed
146
147
    auto constraintsOnDomain = constraints1;
    constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0));
Dennis Gläser's avatar
Dennis Gläser committed
148

149
150
151
152
153
154
155



    ///////////////////////////
    // 4. Network generation //
    ///////////////////////////

156
157
158
    // Helper class to store an arbitrary number
    // of entity sets of different geometry types.
    MultiGeometryEntitySet<Disk, Quad> entitySets;
Dennis Gläser's avatar
Dennis Gläser committed
159

160
161
162
    // Helper class for terminal output of the creation
    // progress and definition of stop criterion etc
    SamplingStatus status;
Dennis Gläser's avatar
Dennis Gläser committed
163
164
    status.setTargetCount(diskSetId, 12); // we want 11 entities of orientation 1
    status.setTargetCount(quadSetId, 16); // we want 13 entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
165

166
    // The actual network generation loop
167
168
169
170
171
172
173
    ctype containedNetworkArea = 0.0;
    while (!status.finished())
    {
        // randomly sample a geometry and obtain the id
        // of the entity set for which it was sampled.
        Id id;
        auto geom = multiSampler(id);
Dennis Gläser's avatar
Dennis Gläser committed
174

175
176
177
178
        // If the set this geometry belongs to is finished, skip the rest
        if (status.finished(id))
        { status.increaseRejectedCounter(); continue; }

179
        // Moreover, we want to avoid small fragments (< 250 m²)
180
        const auto containedArea = computeContainedMagnitude(geom, networkDomain);
Dennis Gläser's avatar
Dennis Gläser committed
181
        if (containedArea < 350.0)
182
183
184
185
186
        { status.increaseRejectedCounter(); continue; }

        // enforce constraints w.r.t. to the other entities
        if (!constraintsMatrix.evaluate(entitySets, geom, id))
        { status.increaseRejectedCounter(); continue; }
Dennis Gläser's avatar
Dennis Gläser committed
187
188

        // enforce constraints w.r.t. the domain boundaries
189
190
191
        if (!constraintsOnDomain.evaluate(domainBoundaryFaces, geom))
        { status.increaseRejectedCounter(); continue; }

192
        // the geometry is admissible
193
194
195
196
197
198
        entitySets.addEntity(geom, id);
        status.increaseCounter(id);
        status.print();

        // keep track of entity area
        containedNetworkArea += containedArea;
Dennis Gläser's avatar
Dennis Gläser committed
199
200
    }

201
    // print the final entity density
202
203
    const auto density = containedNetworkArea/domainVolume;
    std::cout << "\nEntity density of the contained network: " << density << " m²/m³" << std::endl;
Dennis Gläser's avatar
Dennis Gläser committed
204

205
206
207



Dennis Gläser's avatar
Dennis Gläser committed
208
    //////////////////////////////////////////////////////////////////////////
209
    // 5. The entities of the network have been created. We can now         //
Dennis Gläser's avatar
Dennis Gläser committed
210
211
212
    //    construct different types of networks (contained, confined, etc.) //
    //////////////////////////////////////////////////////////////////////////

213
214
    // construct and write a contained network, i.e. write out both network and domain.
    std::cout << "Building and writing contained, confined network" << std::endl;
215
    ContainedEntityNetworkBuilder<ctype> builder;
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251

    // add sub-domains
    builder.addConfiningSubDomain(solids[0],     Id(1));
    builder.addConfiningSubDomain(networkDomain, Id(2));
    builder.addConfiningSubDomain(solids[2],     Id(3));

    // The entites that we created all belong to subdomain 2
    // Use the convenience function to add all entities to the network builder.
    // The last argument provides the id of the sub-domain to which the entities belong.
    entitySets.exportEntitySets(builder, Id(2));

    // Note that one can also call:
    // entitySets.exportEntitySets({diskSetId, quadSetId}, containedConfinedBuilder, Id(2));
    // This overload can be used to define specific entity sets that should be added to the builder

    // now we can build and write out the network in Gmsh file format
    GmshWriter gmshWriter(builder.build());
    gmshWriter.write("contained_confined", // body of the filename to be used (will add .geo)
                     2.5,                  // mesh size to be used on entities
                     5.0);                 // mesh size to be used on domain boundaries

    // we can also not confine the network to its sub-domain,
    // simply by adding the sub-domains as non-confining
    std::cout << "Building and writing contained, unconfined network" << std::endl;
    builder.clear();
    builder.addSubDomain(solids[0],     Id(1));
    builder.addSubDomain(networkDomain, Id(2));
    builder.addSubDomain(solids[2],     Id(3));
    entitySets.exportEntitySets(builder, Id(2));

    gmshWriter = GmshWriter(builder.build());
    gmshWriter.write("contained_unconfined", 2.5, 5.0);

    // We could also only write out the network, without the domain
    // For example, confining the network to the sub-domain...
    std::cout << "Building and writing uncontained, confined network" << std::endl;
252
    EntityNetworkBuilder<ctype> uncontainedBuilder;
253
254
255
256
257
258
259
260
261
262
263
    uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2));
    entitySets.exportEntitySets(uncontainedBuilder, Id(2));
    gmshWriter = GmshWriter(uncontainedBuilder.build());
    gmshWriter.write("uncontained_confined", 2.5);

    // ... or not confining it
    std::cout << "Building and writing uncontained, unconfined network" << std::endl;
    uncontainedBuilder.clear();
    entitySets.exportEntitySets(uncontainedBuilder);
    gmshWriter = GmshWriter(uncontainedBuilder.build());
    gmshWriter.write("uncontained_unconfined", 2.5);
Dennis Gläser's avatar
Dennis Gläser committed
264
265
266

    return 0;
}