example1.cc 6.42 KB
Newer Older
Dennis Gläser's avatar
Dennis Gläser committed
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
#include <random>
#include <iostream>

#include <frackit/common/id.hh>
#include <frackit/io/gmshwriter.hh>

#include <frackit/sampling/makeuniformpointsampler.hh>
#include <frackit/sampling/quadrilateralsampler.hh>
#include <frackit/sampling/status.hh>

#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/box.hh>

#include <frackit/entitynetwork/constraints.hh>
#include <frackit/entitynetwork/networkbuilder.hh>

//! create a network of 3d quadrilaterals
int main()
{
    using namespace Frackit;

    // We consider 3d space here
    static constexpr int worldDimension = 3;

    // Define the type used for coordinates
    using ctype = double;

    // Internal geometry type for 3d quadrilaterals
    using Quad = Quadrilateral<ctype, worldDimension>;

    // Define a domain (here: unit cube) in which the quadrilaterals should be created.
    // Boxes are created by providing xmin, ymin, zmin and xmax, ymax and zmax in constructor.
    Box<ctype> domain(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);

    // We now create a sampler instance that uniformly samples points within this box.
    // These points will be used as the quadrilateral centers.
    auto pointSampler = makeUniformPointSampler(domain);

    // Sampler class for quadrilaterals. Per default, this uses
    // uniform distributions for all parameters defining the quadrilaterals.
    // Quadrilateral samplers require distributions for strike angle, dip angle,
    // edge length (see class description for more details). Moreover, we define
    // a minimum edge length.
    using QuadSampler = QuadrilateralSampler<worldDimension>;
    using Distro = std::normal_distribution<ctype>;
    QuadSampler quadSampler1(pointSampler,
                             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(0.5, 0.1),                         // edge length: mean value & standard deviation
                             0.05);                                    // threshold for minimum edge length

    // We use this sampler to create quadrilaterals based on the distributions.
    // However, we want to create another set of quadrilaterals, whose entities
    // are approximately orthogonal to those defined with the above sampler.
    QuadSampler quadSampler2(makeUniformPointSampler(domain),         // use a new point sampler instance!
                             Distro(toRadians(0.0), toRadians(5.0)),  // strike angle: mean value & standard deviation
                             Distro(toRadians(0.0), toRadians(5.0)),  // dip angle: mean value & standard deviation
                             Distro(0.5, 0.1),                        // edge length: mean value & standard deviation
                             0.05);                                   // threshold for minimum edge length

    // We want to enforce some constraints on the set of quadrilaterals.
    // In particular, for entities of the same set we want a minimum spacing
    // distance of 5cm, and the quadrilaterals must not intersect in angles
    // less than 30°. Moreover, if they intersect, we don't want intersection
    // edges whose length is smaller than 5cm, and, the intersection should not
    // be too close to the boundary of one of two intersecting quadrilaterals. Here: 5cm.
    EntityNetworkConstraints constraintsOnSelf;
    constraintsOnSelf.setMinDistance(0.05);
    constraintsOnSelf.setMinIntersectingAngle(toRadians(30.0));
    constraintsOnSelf.setMinIntersectionMagnitude(0.05);
    constraintsOnSelf.setMinIntersectionDistance(0.05);

    // with respect to entities of the other set, we want to have larger intersection angles
    EntityNetworkConstraints constraintsOnOther;
    constraintsOnOther.setMinDistance(0.05);
    constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
    constraintsOnOther.setMinIntersectionMagnitude(0.05);
    constraintsOnOther.setMinIntersectionDistance(0.05);

    // container to store created entities
    std::vector<Quad> entitySet1;
    std::vector<Quad> entitySet2;

    // we give unique identifiers to both entity sets
    const Id idSet1(1);
    const Id idSet2(2);

    // use the status class to define when to stop sampling
    SamplingStatus status;
    status.setTargetCount(idSet1, 10); // we want 10 entities in set 1
    status.setTargetCount(idSet2, 10); // we want 10 entities in set 2

    // start sampling into set 1 and keep alternating
    bool sampleIntoSet1 = true;

    std::cout << "\n --- Start entity sampling ---\n" << std::endl;
    while (!status.finished())
    {
        // sample a quadrilateral, alternating between sampler 1 and sampler 2
        auto quad = sampleIntoSet1 ? quadSampler1() : quadSampler2();
        auto& entitySet = sampleIntoSet1 ? entitySet1 : entitySet2;
        const auto& otherEntitySet = sampleIntoSet1 ? entitySet2 : entitySet1;

        // sample again if constraints w.r.t. other
        // entities of this set are not fulfilled
        if (!constraintsOnSelf.evaluate(entitySet, quad))
        { status.increaseRejectedCounter(); continue; }

        // sample again if constraints w.r.t. other
        // entities of the other set are not fulfilled
        if (!constraintsOnOther.evaluate(otherEntitySet, quad))
        { status.increaseRejectedCounter(); continue; }

        // the quadrilateral is admissible
        entitySet.push_back(quad);

        // tell the status we have a new entity in this set
        const auto& setId = sampleIntoSet1 ? idSet1 : idSet2;
        status.increaseCounter(setId);
        status.print();

        // sample into the other set the next time
        sampleIntoSet1 = !sampleIntoSet1;
    }
    std::cout << "\n --- Finished entity sampling ---\n" << std::endl;

    // We can now create an entity network from the two sets
    EntityNetworkBuilder builder;
    builder.addEntities(entitySet1);
    builder.addEntities(entitySet2);

    const auto network = builder.build();

    // This can be written out in Gmsh (.geo) format to be
    // meshed by a two-dimensional surface mesh
    std::cout << "\n --- Writing .geo file ---\n" << std::endl;
    GmshWriter writer(network);
    writer.write("network", // filename of the .geo files (will add extension .geo automatically)
                 0.1);      // element size to be used
    std::cout << "\n --- Finished writing .geo file ---\n" << std::endl;

    return 0;
}