README.md 10.2 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/example2_network.png" alt="frackit example 2" width="800"/>
Dennis Gläser's avatar
Dennis Gläser committed
4
5
6
7
8
</p>

Example 2
=========

Dennis Gläser's avatar
Dennis Gläser committed
9
__Note that this description focuses on the C++ implementation given in the main
10
11
12
file `example2.cc`, but in `example2.py` it is illustrated 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
<b> _In this example, you will learn how to:_ </b>

* define a cylindrical domain
* compute the part of entities that is contained inside the domain
* confine the generated entity network onto the domain

Dennis Gläser's avatar
Dennis Gläser committed
19
20
In this exemplary application, a network of quadrilaterals, following the same
distributions as the network of [example 1][0], is created. However, in this
21
example we want to confine the network to a cylindrical domain (see image above).
Dennis Gläser's avatar
Dennis Gläser committed
22
23
24
25
26
27
28
29
30
31
We represent the domain via an instance of the internal cylinder class:

```cpp
// In this example we consider a cylindrical domain
Cylinder<ctype> domain(/*radius*/0.5, /*height*/1.0);
```

This creates an axis-aligned cylinder where the x- and y- axis form the basis of
the bottom face of the cylinder, while the cylinder axis is aligned with the
z-axis. Note that there are other ways to construct arbitrarily-oriented cylinders,
Dennis Gläser's avatar
Dennis Gläser committed
32
and for details on this we refer to the [class documentation][1].
Dennis Gläser's avatar
Dennis Gläser committed
33
34
35
36
37
38
39
40
41
42

In addition to constraints among entities of the different orientations, we also
want to enforce constraints with respect to the boundary of the cylindrical domain.
In this example, we simply reuse the constraints of entities of the same orientation
and write

```cpp
auto constraintsOnBoundary = constraintsOnSelf;
```

43
to create a new constraints object. We create this copy here to improve readability
44
45
of the code when evaluating the constraints. Within the loop in which the entities
are created, all constraints are enforced by writing:
Dennis Gläser's avatar
Dennis Gläser committed
46
47

```cpp
48
49
50
51
52
53
if (!constraintsOnSelf.evaluate(entitySet, quad))
{ status.increaseRejectedCounter(); continue; }

if (!constraintsOnOther.evaluate(otherEntitySet, quad))
{ status.increaseRejectedCounter(); continue; }

Dennis Gläser's avatar
Dennis Gläser committed
54
55
56
57
58
// we want to enforce constraints also w.r.t. to the cylinder boundary
if (!constraintsOnBoundary.evaluate(domain.topFace(), quad))
{ status.increaseRejectedCounter(); continue; }
if (!constraintsOnBoundary.evaluate(domain.bottomFace(), quad))
{ status.increaseRejectedCounter(); continue; }
Dennis Gläser's avatar
Dennis Gläser committed
59
if (!constraintsOnBoundary.evaluate(OCCUtilities::getShape(domain.lateralFace()), quad))
Dennis Gläser's avatar
Dennis Gläser committed
60
61
62
{ status.increaseRejectedCounter(); continue; }
```

63
64
65
66
In the first two if clauses, the constraints against other entities are checked, while
subsequently, the constraints against the boundary faces of the domain cylinder are evaluated.
As you can see, the constraints against the other entities are evaluated on the basis of the
raw entities. This can potentially fail to detect small length scales that appear after confining
Dennis Gläser's avatar
Dennis Gläser committed
67
68
the entities to the domain. A possibilty to address this and make the algorithm more robust
(but computationally more demanding) is discussed at the end of this page.
69
70
71
72

As you can see in the above code snippet, the `Cylinder` class provides functions for obtaining th
representations of the top, bottom and lateral boundaries. The top and bottom boundaries are
represented by instances of the `Disk` class, while the lateral surface is described by the class
Dennis Gläser's avatar
Dennis Gläser committed
73
74
75
76
77
78
79
80
81
82
`CylinderSurface`. Please also note that for the lateral surface, we enforce the constraints
using its representation in OpenCascade, which we obtain by calling `OCCUtilities::getShape`.
We could also write `if (!constraintsOnBoundary.evaluate(domain.lateralFace(), quad)) {...}`,
but the behaviour is different. This is because in the OpenCascade representation of cylindrical
surfaces, there is an edge connecting the upper and lower circles. This has implications on the
constraint regarding the minimum allowed length of intersection edges, as these possibly
intersect with the edge of the cylindrical surface. The figure below illustrates the intersections
of a fracture quadrilateral with both the internal and the OpenCascade representations of cylindrical
surfaces.

83
<!--- Intersection with cylinder, internal vs OCC representation --->
Dennis Gläser's avatar
Dennis Gläser committed
84
85
86
87
88
89
90
91
92
93
<p align="center">
    <img src="../../doc/img/example2_isection.png" alt="frackit example 2 - cylindrical surface intersection" width="800"/>
</p>

As you can see, when the internal representation of cylindrical surfaces is used, the result
is a single intersection edge, while this edge is split into two sub-edges in case the
OpenCascade representation is used. As a result, when evaluating the constraints against
the OpenCascade representation of the cylindrical surface, the constraint on the minimum
allowed length of intersection edges is checked for each sub-edge. Thus, this approach should
be preferred when using cylindrical domains in order to detect small length scales originating
Dennis Gläser's avatar
Dennis Gläser committed
94
95
96
from intersections with the auxiliary edge used by OpenCascade. This is important because
ultimately, when the network and the domain are fragmented and written into a CAD file, the
OpenCascade representation is used.
Dennis Gläser's avatar
Dennis Gläser committed
97
98

Moreover, we want to reject all those quadrilaterals of which only
Dennis Gläser's avatar
Dennis Gläser committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
a very small portion is inside the cylinder. This is done in the lines

```cpp
const auto containedArea = computeContainedMagnitude(quad, domain);

// reject if this is too small (< 0.01 m^2)
if (containedArea < 0.01)
{ status.increaseRejectedCounter(); continue; }
```

where `computeContainedMagnitude()` is a free function that returns the length/area/volume
of the part of a geometry that is contained inside of another geometry.

Finally, after the desired number of entities has been created, we cast the entities
Dennis Gläser's avatar
Dennis Gläser committed
113
into an instance of the class `ContainedEntityNetwork`, using the corresponding builder class:
Dennis Gläser's avatar
Dennis Gläser committed
114
115

```cpp
116
ContainedEntityNetworkBuilder<ctype> builder;
Dennis Gläser's avatar
Dennis Gläser committed
117
118
119
120
121
122
123

// define the domain (single sub-domain) and give it a unique id
builder.addConfiningSubDomain(domain, Id(1));
```

The class `ContainedEntityNetworkBuilder` returns an instance of a `ContainedEntityNetwork`
when the `build()` function is called (see below). This network implementation contains information
124
on sub-domains, and which entities are embedded in which sub-domain. Each sub-domain
Dennis Gläser's avatar
Dennis Gläser committed
125
receives a unique identifier by means of the `Id` class. By calling
126
127
128
`builder.addConfiningSubDomain(domain, Id(1))`, we define the sub-domain to be
confining, i.e. entities that are added to this sub-domain will be confined to it
by cutting away all parts that lie outside the sub-domain. In contrast to that,
Dennis Gläser's avatar
Dennis Gläser committed
129
one could call `builder.addSubDomain(domain, Id(1))`, in which case embedded networks
130
will not be confined (see [example 3][2]). Entities are associated with the sub-
Dennis Gläser's avatar
Dennis Gläser committed
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
domain they are embedded in, and are added to the builder class by writing

```cpp
// define entities to be embedded in this domain
builder.addSubDomainEntities(entitySet1, Id(1));
builder.addSubDomainEntities(entitySet2, Id(1));

const auto network = builder.build();
```

The variable `network` now holds an instance of the `ContainedEntityNetwork`,
which we can also pass to the `GmshWriter`:

```cpp
GmshWriter writer(network);
writer.write("network", // filename of the .geo files (will add extension .geo automatically)
             0.1,       // element size to be used on the quadrilaterals
             0.2);      // element size to be used in the domain away from the quadrilaterals
```

As you can see, one can specify different mesh sizes to be used on the fracture
entities and in the rest of the domain.

154
155
# Enforcing constraints on confined entities

Dennis Gläser's avatar
Dennis Gläser committed
156
As mentioned before, enforcing the constraints between the raw entities, small length
157
158
159
160
161
162
163
164
165
166
167
168
scales appearing after confining the entities to the domain are not detected. An
illustration of such a situation is depicted in the following figure:

<!--- Intersection of two quadrilateraly inside a cylinder, unconfined vs confined --->
<p align="center">
    <img src="../../doc/img/example2_constraints.png" alt="frackit example 2 - unconfined vs confined intersections" width="800"/>
</p>

In such a case, the constraint on the minimum allowed length of intersection segments
might be fulfilled when checked on the raw entities, but might not be fulfilled on the
final geometry after confinement of all entities to the domain.
In order to overcome this, one can also evaluate the constraints on the confined entities,
Dennis Gläser's avatar
Dennis Gläser committed
169
and thus, using the final geometry of an entity after its confinement. This could be achieved
170
171
172
173
174
175
176
177
178
by adjusting the part of the code of this example where the constraints are evaluated in the following way:

```cpp
auto quad = sampleIntoSet1 ? quadSampler1() : quadSampler2();
auto& entitySet = sampleIntoSet1 ? entitySet1 : entitySet2;
const auto& otherEntitySet = sampleIntoSet1 ? entitySet2 : entitySet1;

// get the face(s) that describe(s) the part of the entity that is contained in the domain
using namespace OCCUtilities;
179
const auto containedShape = intersect(getShape(quad), getShape(domain), /*epsilon*/1e-6);
180
181
182
const auto containedFaces = getFaces(containedShape);

// check constraints for all faces making up the contained part
183
if (!constraintsOnSelf.evaluate(entitySet, containedFaces))
184
185
{ status.increaseRejectedCounter(); continue; }

186
if (!constraintsOnOther.evaluate(otherEntitySet, containedFaces))
187
188
189
190
{ status.increaseRejectedCounter(); continue; }
```

Note that the variable `containedFaces` is of type `std::vector<TopoDS_Face>`,
Dennis Gläser's avatar
Dennis Gläser committed
191
where `TopoDS_Face` is the OpenCascade type used to represent general two-dimensional
192
faces. In the above code snippets we have assumed `entitySet` and `otherEntitySet`
193
to also be of type `std::vector<TopoDS_Face>`, and to be composed of the previously
Dennis Gläser's avatar
Dennis Gläser committed
194
accepted faces (also confined to the domain).
195

196
197
198
199
200
201
An effective way to remove very small fragments or edges that were not "detected"
by the constraints, is to adjust the epsilon used for the fragmentation.
You can set a custom epsilon in the builder class by calling its `setEpsilon`
function. Vertices, for instance, that are closer than the defined epsilon are
snapped into a single vertex.

Dennis Gläser's avatar
Dennis Gläser committed
202
203
| [:arrow_right: Go to example 3](https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/appl/example3) |
|---:|
204

Dennis Gläser's avatar
Dennis Gläser committed
205
206
207
208

[0]: https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/appl/example1/README.md
[1]: https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/frackit/geometry/cylinder.hh
[2]: https://git.iws.uni-stuttgart.de/tools/frackit/tree/master/appl/example3/README.md