README.md 12.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
16
17
18
19
20
21
22
<b> _In this example, you will learn how to:_ </b>

* read in a complex domain geometry from an CAD file
* get the bounding box of an arbitrary geometry
* use the `MultiGeometrySampler` and `MultiGeometryEntitySet` classes to conveniently compose a network out of multiple entity geometries
* use the `ConstraintsMatrix` class to conveniently define several constraints to be fulfilled among different entity sets
* define multiple confining/non-confining sub-domains to place the entities in

### Read in a domain geometry

Dennis Gläser's avatar
Dennis Gläser committed
23
24
25
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
26
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
27
28
29
30
31
32
33
34
35
36
`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),
Dennis Gläser's avatar
Dennis Gläser committed
37
and here we want to construct an entity network only in the center one. The variable
Dennis Gläser's avatar
Dennis Gläser committed
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
`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.
Dennis Gläser's avatar
Dennis Gläser committed
54
55
56

### Define entity samplers

Dennis Gläser's avatar
Dennis Gläser committed
57
58
59
60
61
62
63
64
65
66
67
68
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
69
70
71
72
73
74
75
76
77
// 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
78
79
80
81
82
83
```

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
84
85
86
87
88
89
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
90
91
92
93
94
95
96
97
98

```

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
Dennis Gläser's avatar
Dennis Gläser committed
99
classes with possibly varying geometry types. Arbitrary many samplers can be added to
Dennis Gläser's avatar
Dennis Gläser committed
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
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);
```

Dennis Gläser's avatar
Dennis Gläser committed
118
119
Each sampler is associated with a unique identifier and an error is thrown if one
tries to add a sampler with an identifier that is already taken. Sampling
Dennis Gläser's avatar
Dennis Gläser committed
120
121
122
123
124
125
126
127
128
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);
```

Dennis Gläser's avatar
Dennis Gläser committed
129
130
131
132
133
134
135
136
137
stores in `id` the identifier of the sampler from which `geom` was sampled
(in this case this is either `diskSetId` or `quadSetId`). Note that, per default,
a `SequentialSamplingStrategy` is used, which means that the provided ids are used
successively in the calls to the `()` operator. However, users can implement their
own strategies and pass them to the constructor of `MultiGeometrySampler`.
In the above code snippet, the variable `geom` holds an instance of an abstract
geometry class as the return type of the `()` operator of the `MultiGeometrySampler`
must be uniquely 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
138
139
140

In this context, another useful class is the `MultiGeometryEntitySet`. It can
store arbitrarily many sets of entities of different types. These
Dennis Gläser's avatar
Dennis Gläser committed
141
142
types have to be provided as template arguments to the class. In this example,
we use `Quadrilateral`s and `Disk`s, and therefore, we construct the class as follows:
Dennis Gläser's avatar
Dennis Gläser committed
143
144
145
146
147

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

Dennis Gläser's avatar
Dennis Gläser committed
148
Entities are again added by passing a unique identifier, that is, the code
Dennis Gläser's avatar
Dennis Gläser committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

```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
166
This allows the generation and storage of multiple geometry types and multiple orientations
Dennis Gläser's avatar
Dennis Gläser committed
167
168
in a compact way.

Dennis Gläser's avatar
Dennis Gläser committed
169
170
### Constraints definitions

Dennis Gläser's avatar
Dennis Gläser committed
171
In the previous examples, the constraints for a new entity candidate were evaluated
Dennis Gläser's avatar
Dennis Gläser committed
172
manually against the previously admitted entities obtained from the different samplers.
Dennis Gläser's avatar
Dennis Gläser committed
173
This becomes increasingly cumbersome the more orientations are to be considered, especially
Dennis Gläser's avatar
Dennis Gläser committed
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
204
205
206
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
Dennis Gläser's avatar
Dennis Gläser committed
207
208
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
Dennis Gläser's avatar
Dennis Gläser committed
209
210
211
212
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`.

Dennis Gläser's avatar
Dennis Gläser committed
213
214
### Network construction

Dennis Gläser's avatar
Dennis Gläser committed
215
216
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
217
file format. In a first step, we create an instance of the network builder
Dennis Gläser's avatar
Dennis Gläser committed
218
219
220
class for contained networks and define the sub-domains:

```cpp
221
ContainedEntityNetworkBuilder<ctype> builder;
Dennis Gläser's avatar
Dennis Gläser committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237

// 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
238
be added to the builder (see [example 2][0]). If `entitySets` holds entity sets
Dennis Gläser's avatar
Dennis Gläser committed
239
240
241
242
243
244
245
246
247
248
249
250
251
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
252
253
| [: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
254

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