README.md 11.1 KB
Newer Older
Dennis Gläser's avatar
Dennis Gläser committed
1
2
<!--- Example picture --->
<p align="center">
Dennis Gläser's avatar
Dennis Gläser committed
3
    <img src="../../doc/img/example3_network.png" alt="frackit example 3" width="800"/>
Dennis Gläser's avatar
Dennis Gläser committed
4
5
6
7
8
</p>

Example 3
=========

Dennis Gläser's avatar
Dennis Gläser committed
9
__As in the previous examples, this description focuses on the C++ implementation
10
11
12
given in the main file `example3.cc`, but in `example3.py` you can find how to
realize this example using the Frackit python bindings.__

Dennis Gläser's avatar
Dennis Gläser committed
13
14
15
In contrast to the [previous example][0], we now want to create a network within
a complex domain that cannot be represented by the internal geometry classes.
The domain has been created with [Gmsh][1] and has been saved in .brep file format.
Dennis Gläser's avatar
Dennis Gläser committed
16
In a first step, we need to read the file and parse it into an instance of the
Dennis Gläser's avatar
Dennis Gläser committed
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
`TopoDS_Shape` class of [OpenCascade][2]. To this end, we use the utility function
provided by Frackit in the namespace OCCUtilities:

```cpp
const auto domainShape = OCCUtilities::readShape(BREPFILE);
```

BREPFILE is a preprocessor variable that is substituted by cmake with the actual
path to the .brep file used for this test (see the file `CMakeLists.txt` in this folder).
The domain specified in the .brep file consists of three layers (see above image),
and we want to construct an entity network only in the center one. The variable
`domainShape` contains the information of all layers, and we can use the utility
function

```cpp
const auto solids = OCCUtilities::getSolids(domainShape);
```

to extract the three layers, represented by instances of the `TopoDS_Solid` class
of [OpenCascade][2]. Note that the function `getSolids()` returns an instance of
`std::vector<TopoDS_Solid>`, and thus, we obtain the middle layer with

```cpp
const auto& networkDomain = solids[1];
```

Note that this requires knowledge about the ordering of the solids in the .brep file.
As in the previous examples, we need a point sampler that samples the points to be
used as the center points of the fracture entities. In this example, we again want to
uniformly sample the points within a box, for which we use the bounding box of
`networkDomain`. This can be obtained using the utility function:

```cpp
const auto domainBBox = OCCUtilities::getBoundingBox(networkDomain);
```

With this, we create a quadrilateral sampler in a similar way as in the previous examples

```cpp
Dennis Gläser's avatar
Dennis Gläser committed
56
57
58
59
60
61
62
63
64
// 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>;

QuadrilateralSampler<3> quadSampler(makeUniformPointSampler(domainBBox),           // sampler for quadrilateral center points
                                    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
Dennis Gläser's avatar
Dennis Gläser committed
65
66
67
68
69
70
```

For the entities of the other orientation, we want to use `Disk` objects in this case.
The instantiation of the corrensponding sampler class looks like this:

```cpp
Dennis Gläser's avatar
Dennis Gläser committed
71
72
73
74
75
76
DiskSampler diskSampler(makeUniformPointSampler(domainBBox),           // sampler for disk center points
                        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
Dennis Gläser's avatar
Dennis Gläser committed
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

```

Here, the second and third constructor arguments specify the distributions used
for the major and minor axis lengths of the elliptical disks. The last three
constructor arguments specify the distributions used to determine the orientation
of the disks (for details we refer to the [class documentation][3]).

We use the `MultiGeometrySampler` class to facilitate sampling from several sampler
class with possibly varying geometry types. Arbitrary many samplers can be added to
this class, with the following syntax:

```cpp
// 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

using Disk = Disk<ctype>;
using Quad = Quadrilateral<ctype, 3>;

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

Each sampler is associated with a unique identifier and an error is thrown if it
is tried to add a sampler with an identifier that is already taken. Sampling
occurs again by using the `()` operator, but, in this case it receives an instance
of `Id`, in which it will store the identifier of the sampler from which the
geometry has been generated. That is, the code

```cpp
Id id;
auto geom = multiSampler(id);
```

stores in `id` the identifier of the sampler from which `geom` was sampled.
The variable `geom` holds an instance of an abstract geometry class as the
return type of the `()` operator of the `MultiGeometrySampler` must be uniquely
Dennis Gläser's avatar
Dennis Gläser committed
119
120
defined. However, instances of the abstract geometry class can be cast back into
the actual geometry (in this case `Disk` or `Quad`).
Dennis Gläser's avatar
Dennis Gläser committed
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

In this context, another useful class is the `MultiGeometryEntitySet`. It can
store arbitrarily many sets of entities of different types. These
types have to be provided as template arguments to the class:

```cpp
MultiGeometryEntitySet<Disk, Quad> entitySets;
```

Entities are added passing a unique identifier, that is, the code

```cpp
Disk disk = diskSampler();
entitySets.addEntity(disk, diskSetId);
```

adds a disk to an entity set with the identifier `diskSetId`. A respective overload
for this function is also implemented for the abstract geometry class mentioned above.
Thus, the geometries sampled from the `MultiGeometrySampler` can directly be inserted
into the instance of `MultiGeometryEntitySet`:

```cpp
Id id;
auto geom = multiSampler(id);
entitySets.addEntity(geom, id);
```

Dennis Gläser's avatar
Dennis Gläser committed
148
This allows the generation and storage of multiple geometry types and multiple orientations
Dennis Gläser's avatar
Dennis Gläser committed
149
150
in a compact way.

Dennis Gläser's avatar
Dennis Gläser committed
151
In the previous examples, the constraints for a new entity candidate were evaluated
Dennis Gläser's avatar
Dennis Gläser committed
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
manually against the previously admitted entities obtained from the different samplers.
This becomes increasingly cumbersome the more directions are to be considered, especially
if individual constraints are chosen between entities obtained from different samplers.
In such cases, the `EntityNetworkConstraintsMatrix` class can be used to facilitate
the constraints evaluation. It can be used in conjunction with the `MultiGeometryEntitySet`
class and allows for the definition of an arbitrary number of constraints,
which are to be evaluated among entities of specific orientations. In this example, we
create three different instances of the `EntityNetworkConstraints` class, and add them
to the matrix with:

```cpp
EntityNetworkConstraintsMatrix<Constraints> constraintsMatrix;
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

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
```

The `IdPair` objects are used to tell the constraints matrix between entities of
which sets, identified by ids, the individual constraints must hold. The evaluation
of all constraints can then be written in one line:

```cpp
// enforce constraints w.r.t. to the other entities
if (!constraintsMatrix.evaluate(entitySets, geom, id))
{ status.increaseRejectedCounter(); continue; }
```

In this function call, `id` holds the identifier of the set for which the candidate
`geom` was sampled. Internally, the `EntityNetworkConstraintsMatrix` instance will
check `geom` against the other entity sets using all constraints that were defined for this identifier. For instance, in this example a disk-shaped geometry sampled from
the sampler with id `diskSetId`, will be checked against all entities of the set with
id `diskSetId` using `constraints1` and against all entities of the set with id `quadSetId`
using `constraintsOnOther`.

After successful generation of the desired number of entities, we again want to
construct an entity network from the raw entities and write it out in [Gmsh][1]
Dennis Gläser's avatar
Dennis Gläser committed
194
file format. In a first step, we create an instance of the network builder
Dennis Gläser's avatar
Dennis Gläser committed
195
196
197
class for contained networks and define the sub-domains:

```cpp
198
ContainedEntityNetworkBuilder<ctype> builder;
Dennis Gläser's avatar
Dennis Gläser committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

// add sub-domains
builder.addConfiningSubDomain(solids[0],     Id(1));
builder.addConfiningSubDomain(networkDomain, Id(2));
builder.addConfiningSubDomain(solids[2],     Id(3));
```

We have given the identifier with index 2 to the (sub-)domain for which we have created
the entity network. Using the convenience function available in the `MultiGeometryEntitySet`
class, we can add all entities to the builder by writing

```cpp
entitySets.exportEntitySets(builder, Id(2));
```

where the second argument states the (sub-)domain for which the entities should
Dennis Gläser's avatar
Dennis Gläser committed
215
be added to the builder (see [example 2][0]). If `entitySets` holds entity sets
Dennis Gläser's avatar
Dennis Gläser committed
216
217
218
219
220
221
222
223
224
225
226
227
228
associated with multiple (sub-)domains, or, if only some of the entity sets should
be added to the builder, one can use:

```cpp
entitySets.exportEntitySets({diskSetId, quadSetId}, builder, Id(2));
```

With this call only the entities of those sets that are given in the list of ids
inside the brackets (`{}`) are exported. An entity network can then be built and
written out in the same manner as shown in the other examples. At the end of the
main file `example3.cc`, it is shown how different types of networks
(contained-confined, contained-unconfined, etc.) can be created from the raw entities.

Dennis Gläser's avatar
Dennis Gläser committed
229
230
| [:arrow_right: Go to example 4](https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/appl/example4) |
|---:|
Dennis Gläser's avatar
Dennis Gläser committed
231

Dennis Gläser's avatar
Dennis Gläser committed
232
[0]: https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/appl/example2/README.md
Dennis Gläser's avatar
Dennis Gläser committed
233
234
[1]: http://gmsh.info/
[2]: https://www.opencascade.com/content/download-center
Dennis Gläser's avatar
Dennis Gläser committed
235
[3]: https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/frackit/sampling/disksampler.hh