example3.cc 12.3 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 -> normal distributions for all parameters
    using Distro = std::normal_distribution<ctype>;

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
Dennis Gläser's avatar
Dennis Gläser committed
82
83
84
85
86
                            Distro(30.0, 6.5),                        // major axis length: mean value & standard deviation
                            Distro(24.0, 4.5),                        // minor axis length: mean value & standard deviation
                            Distro(toRadians(0.0), toRadians(7.5)),   // rotation around x-axis: mean value & standard deviation
                            Distro(toRadians(0.0), toRadians(7.5)),   // rotation around y-axis: mean value & standard deviation
                            Distro(toRadians(0.0), toRadians(7.5)));  // rotation around z-axis: mean value & standard deviation
87
88
89

    // sampler for quadrilaterals (orientation 2)
    QuadrilateralSampler<3> quadSampler(makeUniformPointSampler(domainBBox),                               // sampler for quadrilateral center points
Dennis Gläser's avatar
Dennis Gläser committed
90
91
92
                                        Distro(toRadians(45.0), toRadians(5.0)),  // strike angle: mean value & standard deviation
                                        Distro(toRadians(90.0), toRadians(5.0)),  // dip angle: mean value & standard deviation
                                        Distro(45.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
    Constraints constraints1;
Dennis Gläser's avatar
Dennis Gläser committed
117
    constraints1.setMinDistance(2.5);
118
    constraints1.setMinIntersectingAngle(toRadians(25.0));
Dennis Gläser's avatar
Dennis Gläser committed
119
120
    constraints1.setMinIntersectionMagnitude(5.0);
    constraints1.setMinIntersectionDistance(2.5);
121
122

    // Define constraints between entities of orientation 2
Dennis Gläser's avatar
Dennis Gläser committed
123
124
125
    // 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
126

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

Dennis Gläser's avatar
Dennis Gläser committed
132
    // We can use the constraints matrix to facilitate constraint evaluation
133
    EntityNetworkConstraintsMatrix<Constraints> constraintsMatrix;
134
135
136
137
138
    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
139
140
141
142
143
144

    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
145
146
    auto constraintsOnDomain = constraints1;
    constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0));
Dennis Gläser's avatar
Dennis Gläser committed
147

148
149
150
151
152
153
154



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

155
156
157
    // 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
158

159
160
161
    // 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
162
163
    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
164

165
    // The actual network generation loop
166
167
168
169
170
171
172
    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
173

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

178
        // Moreover, we want to avoid small fragments (< 250 m²)
179
        const auto containedArea = computeContainedMagnitude(geom, networkDomain);
Dennis Gläser's avatar
Dennis Gläser committed
180
        if (containedArea < 350.0)
181
182
183
184
185
        { 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
186
187

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

191
        // the geometry is admissible
192
193
194
195
196
197
        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
198
199
    }

200
    // print the final entity density
201
202
    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
203

204
205
206



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

212
213
    // construct and write a contained network, i.e. write out both network and domain.
    std::cout << "Building and writing contained, confined network" << std::endl;
214
    ContainedEntityNetworkBuilder<ctype> builder;
215
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

    // 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;
251
    EntityNetworkBuilder<ctype> uncontainedBuilder;
252
253
254
255
256
257
258
259
260
261
262
    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
263
264
265

    return 0;
}