example3.cc 12.8 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
    // Helper class for terminal output of the creation
    // progress and definition of stop criterion etc
162
163
164
    // Check if a value has been passed via the command line,
    // otherwise use the defaults. This is used in the test suite
    // to speed up testing time.
165
    SamplingStatus status;
166
167
    status.setTargetCount(diskSetId, (argc > 1 ? std::stoi(argv[1]) : 12)); // we want 11 entities of orientation 1
    status.setTargetCount(quadSetId, (argc > 1 ? std::stoi(argv[1]) : 16)); // we want 13 entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
168

169
    // The actual network generation loop
170
171
172
173
174
175
176
    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
177

178
179
        // If the set this geometry belongs to is finished, skip the rest
        if (status.finished(id))
Dennis Gläser's avatar
Dennis Gläser committed
180
        { status.increaseRejectedCounter("set finished"); continue; }
181

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

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

        // enforce constraints w.r.t. the domain boundaries
Dennis Gläser's avatar
Dennis Gläser committed
192
193
        if (const auto res = constraintsOnDomain.evaluate(domainBoundaryFaces, geom); !res)
        { status.increaseRejectedCounter(res.violationLabel() + " (boundary)"); continue; }
194

195
        // the geometry is admissible
196
197
198
199
200
201
        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
202
203
    }

204
    // print the final entity density
205
206
    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
207

Dennis Gläser's avatar
Dennis Gläser committed
208
209
210
    // print info on rejection events
    status.printRejectionData();
    std::cout << std::endl;
211
212


Dennis Gläser's avatar
Dennis Gläser committed
213
    //////////////////////////////////////////////////////////////////////////
214
    // 5. The entities of the network have been created. We can now         //
Dennis Gläser's avatar
Dennis Gläser committed
215
216
217
    //    construct different types of networks (contained, confined, etc.) //
    //////////////////////////////////////////////////////////////////////////

218
219
    // construct and write a contained network, i.e. write out both network and domain.
    std::cout << "Building and writing contained, confined network" << std::endl;
220
    ContainedEntityNetworkBuilder<ctype> builder;
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
252
253
254
255
256

    // 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;
257
    EntityNetworkBuilder<ctype> uncontainedBuilder;
258
259
260
261
262
263
264
265
266
267
268
    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
269
270
271

    return 0;
}