test_generate_disk_network_shape.cc 9.02 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
#include <iostream>
#include <random>
#include <fstream>

#include <BRepTools.hxx>
#include <BRep_Builder.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Compound.hxx>
#include <TopoDS_Solid.hxx>

#include <frackit/geometry/box.hh>
#include <frackit/geometry/disk.hh>
#include <frackit/precision/precision.hh>
#include <frackit/magnitude/magnitude.hh>
#include <frackit/magnitude/containedmagnitude.hh>
#include <frackit/occ/breputilities.hh>

#include <frackit/sampling/pointsampling.hh>
#include <frackit/sampling/geometrysampling.hh>

#include <frackit/entitynetwork/constraints.hh>

//! create a disk network embedded in a cylinder
int main(int argc, char** argv)
{
    //! print welcome message
    std::cout << "\n\n"
              << "##########################################\n"
              << "## Creating a confined network of disks ##\n"
              << "##########################################\n"
              << "\n\n";

    // brep file name is specified in CMakeLists.txt
    std::ifstream domainShapeFile(BREPFILE);
    if (!domainShapeFile)
        throw std::runtime_error(std::string("Could not open shape file"));

    TopoDS_Shape domainShape;
    BRepTools::Read(domainShape, domainShapeFile, BRep_Builder());

    // obtain the (single) solid contained in the file & get its boundaries
    using namespace Frackit;
    const auto solids = OCCUtilities::getSolids(domainShape);

    assert(solids.size() == 1);
    const auto& domain = solids[0];
    const auto domainFaces = OCCUtilities::getFaces(domain);
    const auto domainVolume = computeMagnitude(domain);

    // create the disk samplers
    using ctype = double;
    using Disk = Disk<ctype>;
    using DiskSampler = GeometrySampler<Disk>;
    using PointSampler = GeometryPointSampler< Box<ctype> >;

    // sample points within bounding box of domain
    PointSampler pointSampler(OCCUtilities::getBoundingBox(domain));

    // sampler for disks of orientation 1
    DiskSampler diskSampler_1(std::normal_distribution<ctype>(0.35, 0.1),
                              std::normal_distribution<ctype>(0.225, 0.05),
                              std::normal_distribution<ctype>(toRadians(25.0), toRadians(5.0)),
                              std::normal_distribution<ctype>(toRadians(0.0),  toRadians(5.0)),
                              std::normal_distribution<ctype>(toRadians(45.0), toRadians(5.0)));

    // sampler for disks of orientation 1
    DiskSampler diskSampler_2(std::normal_distribution<ctype>(0.35, 0.1),
                              std::normal_distribution<ctype>(0.225, 0.05),
                              std::normal_distribution<ctype>(toRadians(-35.0), toRadians(5.0)),
                              std::normal_distribution<ctype>(toRadians(0.0),   toRadians(5.0)),
                              std::normal_distribution<ctype>(toRadians(-45.0), toRadians(5.0)));

    //! containers to store created entities
    std::vector<Disk> diskSet1;
    std::vector<Disk> diskSet2;

    //! enforce some constraints on the network
    Frackit::EntityNetworkConstraints<ctype> constraintsOnSelf;
    Frackit::EntityNetworkConstraints<ctype> constraintsOnOther;
    Frackit::EntityNetworkConstraints<ctype> constraintsOnDomain;

    // constraints among disks of the same set
    constraintsOnSelf.setMinDistance(0.1);
    constraintsOnSelf.setMinIntersectingAngle(toRadians(45.0));
    constraintsOnSelf.setMinIntersectionMagnitude(0.01);
    constraintsOnSelf.setMinIntersectionDistance(0.05);

    // constraints with the other disk set
    constraintsOnOther.setMinDistance(0.01);
    constraintsOnOther.setMinIntersectingAngle(toRadians(15.0));
    constraintsOnOther.setMinIntersectionMagnitude(0.01);
    constraintsOnOther.setMinIntersectionDistance(0.05);

    // constraints with the other disk set
    constraintsOnDomain.setMinDistance(0.01);
    constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0));
    constraintsOnDomain.setMinIntersectionMagnitude(0.01);
    constraintsOnDomain.setMinIntersectionDistance(0.05);

    //! create pseudo-random number generator to select set during creation
    std::default_random_engine randomNumber{std::random_device{}()};

    //! keep track of number of created & accepted entities
    std::size_t total = 0;
    std::size_t accepted = 0;
    std::size_t accepted_1 = 0;
    std::size_t accepted_2 = 0;

    //! print header of status output
    std::cout << "##############################################################################################################################\n"
              << "# No. of entities   |   No. rejected entities   |   Acceptance Ratio   |   Current entity density [m²/m³]   |   Progress [%] #\n"
              << "#----------------------------------------------------------------------------------------------------------------------------#\n";

    ctype currentDensity = 0.0;
    ctype currentDiskArea = 0.0;
    const std::size_t numTargetEntities_1 = 5;
    const std::size_t numTargetEntities_2 = 5;
    while (accepted_1 != numTargetEntities_1 || accepted_2 != numTargetEntities_2)
    {
        bool createSecondary = randomNumber()%2;
        if (!createSecondary && accepted_1 == numTargetEntities_1)
            createSecondary = true;
        else if (createSecondary && accepted_2 == numTargetEntities_2)
            createSecondary = false;

        auto disk = createSecondary ? diskSampler_1(pointSampler)
                                    : diskSampler_2(pointSampler);
        total++;

        // We don't want ellipses of too large aspect ratio
        if (disk.majorAxisLength() > disk.minorAxisLength()*3.0) continue;

        auto& diskSetSelf = createSecondary ? diskSet2 : diskSet1;
        const auto& diskSetOther = createSecondary ? diskSet1 : diskSet2;

        // enforce constraints w.r.t. to the own disk set
        if (!constraintsOnSelf.evaluate(diskSetSelf, disk)) continue;
        // enforce constraints w.r.t. to the other disk set
        if (!constraintsOnOther.evaluate(diskSetOther, disk)) continue;

        // // enforce constraints w.r.t. the domain boundaries
        // TODO: INTERSECTION ALGOS WITH SHAPE CLASSES
        // bool violates = false;
        // for (const auto& face : domainFaces)
        //     if (!constraintsOnDomain.evaluate(face, disk))
        //     { violates = true; break; }
        // if (violates) continue;

        // reject if intersection with domain is too small (here: 0.2m²)
        const auto containedArea = computeContainedMagnitude(disk, domain);
        if (containedArea < 0.1) continue;

        // if we get here, the disk is admissible
        diskSetSelf.push_back(std::move(disk));
        accepted++;

        if (createSecondary) accepted_2++;
        else accepted_1++;

        // compute new density (use minimum w.r.t. domain/network volume)
        currentDiskArea += containedArea;
        currentDensity = currentDiskArea/domainVolume;
        std::cout << "  "   << std::setw(18) << std::setfill(' ')
                            << std::to_string(diskSet1.size() + diskSet2.size()) + std::string(9, ' ')
                  << "|   " << std::setprecision(6) << std::setw(24) << std::setfill(' ')
                            << std::to_string(total-accepted) + std::string(12, ' ')
                  << "|   " << std::setprecision(6) << std::setw(19) << std::setfill(' ')
                            << std::to_string( double(double(accepted)/double(total)) ) + std::string(7, ' ')
                  << "|   " << std::setprecision(6) << std::setw(33) << std::setfill(' ')
                            << std::to_string(currentDensity) + std::string(17, ' ')
                  << "|   " << std::string(3, ' ') + std::to_string( double(accepted_1+accepted_2)/double(numTargetEntities_1+numTargetEntities_2) ) << std::endl;
    }

    // create a single compound and write to .brep file
    BRep_Builder b;
    TopoDS_Compound c;
    b.MakeCompound(c);
    for (const auto& disk : diskSet1) b.Add(c, OCCUtilities::getShape(disk));
    for (const auto& disk : diskSet2) b.Add(c, OCCUtilities::getShape(disk));
    b.Add(c, domain);
    BRepTools::Write(c, "disknetwork.brep");

    // fragment the network and write in separate file
    std::vector<TopoDS_Shape> allDisks;
    for (const auto& disk : diskSet1) allDisks.push_back(OCCUtilities::getShape(disk));
    for (const auto& disk : diskSet2) allDisks.push_back(OCCUtilities::getShape(disk));

    const auto fragmentedNetwork = OCCUtilities::fragment(allDisks, Frackit::Precision<ctype>::confusion());
    BRepTools::Write(fragmentedNetwork, "disknetwork_fragmented.brep");

    const auto confinedNetwork = OCCUtilities::intersect(fragmentedNetwork,
                                                         domain,
                                                         Frackit::Precision<ctype>::confusion());
    BRepTools::Write(confinedNetwork, "disknetwork_confined.brep");

    std::vector<TopoDS_Shape> allShapes({confinedNetwork, domain});
    BRepTools::Write(OCCUtilities::fragment(allShapes, Frackit::Precision<ctype>::confusion()),
                    "disknetwork_medium.brep");

    std::cout << "All tests passed" << std::endl;

    return 0;
}