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

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

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

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

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

// 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"
39
40
41
              << "##############################################################################\n"
              << "## Example: Creation of entities (disks/quadrilaterals) in a layered medium ##\n"
              << "##############################################################################\n"
Dennis Gläser's avatar
Dennis Gläser committed
42
43
              << "\n\n";

44
45
46
47
48
49
    // 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
50
51
52
53
    /////////////////////////////////////////////////////
    // 1. Read in the domain geometry from .brep file. //
    //    The file name is defined in CMakeLists.txt   //
    /////////////////////////////////////////////////////
54
    const auto domainShape = OCCUtilities::readShape(BREPFILE);
Dennis Gläser's avatar
Dennis Gläser committed
55
56
57
58

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

61
62
    // 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
63
64
    const auto& networkDomain = solids[1];
    const auto domainVolume = computeMagnitude(networkDomain);
65
66
67
    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
68

69
70
71



72
73
74
75
    //////////////////////////////////////////////////////////////////////////////
    // 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
76

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

80
81
    // sampler for disks (orientation 1)
    DiskSampler diskSampler(makeUniformPointSampler(domainBBox),                               // sampler for disk center points
82
83
                            std::normal_distribution<ctype>(30.0, 6.5),                        // major axis length: mean value & standard deviation
                            std::normal_distribution<ctype>(24.0, 4.5),                        // minor axis length: mean value & standard deviation
84
85
86
87
88
89
90
91
                            std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)),   // rotation around x-axis: mean value & standard deviation
                            std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)),   // rotation around y-axis: mean value & standard deviation
                            std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)));  // rotation around z-axis: mean value & standard deviation

    // sampler for quadrilaterals (orientation 2)
    QuadrilateralSampler<3> quadSampler(makeUniformPointSampler(domainBBox),                               // sampler for quadrilateral center points
                                        std::normal_distribution<ctype>(toRadians(0.0), toRadians(5.0)),   // strike angle: mean value & standard deviation
                                        std::normal_distribution<ctype>(toRadians(45.0), toRadians(5.0)),  // dip angle: mean value & standard deviation
92
                                        std::normal_distribution<ctype>(35.0, 5.0),                        // edge length: mean value & standard deviation
93
94
95
96
97
98
99
100
101
102
103
104
105
                                        5.0);                                                              // threshold for minimum edge length

    // 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);

106
107
108
109
110
111
112
113



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

114
    // Define constraints between entities of orientation 1
115
    using Constraints = EntityNetworkConstraints<ctype>;
116
117
118
119
120
121
122
123
124
125
126
127
    Constraints constraints1;
    constraints1.setMinDistance(1.5);
    constraints1.setMinIntersectingAngle(toRadians(25.0));
    constraints1.setMinIntersectionMagnitude(2.5);
    constraints1.setMinIntersectionDistance(2.0);

    // Define constraints between entities of orientation 2
    Constraints constraints2;
    constraints2.setMinDistance(10.0);
    constraints2.setMinIntersectingAngle(toRadians(25.0));
    constraints2.setMinIntersectionMagnitude(2.5);
    constraints2.setMinIntersectionDistance(2.0);
Dennis Gläser's avatar
Dennis Gläser committed
128

129
130
131
132
    // Define constraints between entities of different sets
    Constraints constraintsOnOther;
    constraintsOnOther.setMinDistance(2.5);
    constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
Dennis Gläser's avatar
Dennis Gläser committed
133
134
135
    constraintsOnOther.setMinIntersectionMagnitude(2.5);
    constraintsOnOther.setMinIntersectionDistance(2.0);

136
137
    // We can use a constraint matrix to facilitate constraint evaluation
    EntityNetworkConstraintsMatrix<Constraints> constraintsMatrix;
138
139
140
141
142
    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
143
144
145
146
147
148
149

    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
    EntityNetworkConstraints constraintsOnDomain;
Dennis Gläser's avatar
Dennis Gläser committed
150
151
152
153
154
    constraintsOnDomain.setMinDistance(1.5);
    constraintsOnDomain.setMinIntersectingAngle(toRadians(5.0));
    constraintsOnDomain.setMinIntersectionMagnitude(20.0);
    constraintsOnDomain.setMinIntersectionDistance(2.0);

155
156
157
158
159
160
161



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

162
163
164
    // 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
165

166
167
168
    // Helper class for terminal output of the creation
    // progress and definition of stop criterion etc
    SamplingStatus status;
169
170
    status.setTargetCount(diskSetId, 10); // we want 10 entities of orientation 1
    status.setTargetCount(quadSetId, 10); // we want 10 entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
171

172
    // The actual network generation loop
173
174
175
176
177
178
179
    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
180

181
182
183
184
        // If the set this geometry belongs to is finished, skip the rest
        if (status.finished(id))
        { status.increaseRejectedCounter(); continue; }

185
        // Moreover, we want to avoid small fragments (< 250 m²)
186
        const auto containedArea = computeContainedMagnitude(geom, networkDomain);
187
        if (containedArea < 250.0)
188
189
190
191
192
        { 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
193
194

        // enforce constraints w.r.t. the domain boundaries
195
196
197
198
199
200
201
202
203
204
        if (!constraintsOnDomain.evaluate(domainBoundaryFaces, geom))
        { status.increaseRejectedCounter(); continue; }

        // the disk is admissible
        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
205
206
    }

207
    // print the final entity density
208
209
    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
210

211
212
213



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

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
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
    // construct and write a contained network, i.e. write out both network and domain.
    std::cout << "Building and writing contained, confined network" << std::endl;
    ContainedEntityNetworkBuilder builder;

    // 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;
    EntityNetworkBuilder uncontainedBuilder;
    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
270
271
272

    return 0;
}