example3.cc 13 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);
Dennis Gläser's avatar
Dennis Gläser committed
62
63
64
    const auto domainShell = OCCUtilities::getShells(networkDomain);
    if (domainShell.size() != 1) throw std::runtime_error("Expected a single shell bounding the domain");
    const auto domainBoundaryFaces = OCCUtilities::getFaces(domainShell[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

Dennis Gläser's avatar
Dennis Gläser committed
160
161
162
163
    // Helper class for terminal output of the creation progress
    // and definition of stop criterion etc. 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.
164
    SamplingStatus status;
Dennis Gläser's avatar
Dennis Gläser committed
165
166
    status.setTargetCount(diskSetId, (argc > 1 ? std::stoi(argv[1]) : 12)); // desired number of entities of orientation 1
    status.setTargetCount(quadSetId, (argc > 1 ? std::stoi(argv[1]) : 16)); // desired number of entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
167

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

177
178
        // 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
179
        { status.increaseRejectedCounter("set finished"); continue; }
180

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

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

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

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

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

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


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

217
218
    // construct and write a contained network, i.e. write out both network and domain.
    std::cout << "Building and writing contained, confined network" << std::endl;
219
    ContainedEntityNetworkBuilder<ctype> builder;
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236

    // 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());
237
238
239
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
    gmshWriter.write("contained_confined");
240
241
242
243
244
245
246
247
248
249
250

    // 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());
251
252
253
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::subDomain, 5.0);
    gmshWriter.write("contained_unconfined");
254
255
256
257

    // 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;
258
    EntityNetworkBuilder<ctype> uncontainedBuilder;
259
260
261
    uncontainedBuilder.addConfiningSubDomain(networkDomain, Id(2));
    entitySets.exportEntitySets(uncontainedBuilder, Id(2));
    gmshWriter = GmshWriter(uncontainedBuilder.build());
262
263
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
    gmshWriter.write("uncontained_confined");
264
265
266
267
268
269

    // ... or not confining it
    std::cout << "Building and writing uncontained, unconfined network" << std::endl;
    uncontainedBuilder.clear();
    entitySets.exportEntitySets(uncontainedBuilder);
    gmshWriter = GmshWriter(uncontainedBuilder.build());
270
271
    gmshWriter.setMeshSize(GmshWriter::GeometryTag::entity, 2.5);
    gmshWriter.write("uncontained_unconfined");
Dennis Gläser's avatar
Dennis Gläser committed
272
273
274

    return 0;
}