README.md 9.85 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
In this exemplary application, a network of quadrilaterals, following the same
distributions as the network of [example 1][0], is created. However, in this
15
example we want to confine the network to a cylindrical domain (see image above).
Dennis Gläser's avatar
Dennis Gläser committed
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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,
and for details on this we refer to the [class documentation][2].

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

37
to create a new constraints object. We create this copy here to improve readability
38
39
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
40
41

```cpp
42
43
44
45
46
47
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
48
49
50
51
52
// 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
53
if (!constraintsOnBoundary.evaluate(OCCUtilities::getShape(domain.lateralFace()), quad))
Dennis Gläser's avatar
Dennis Gläser committed
54
55
56
{ status.increaseRejectedCounter(); continue; }
```

57
58
59
60
61
62
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
the entities to the domain. A possibilty to address this and make the algorithm more robust (but
less efficient) is discussed at the end of this page.

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
67
68
69
70
71
72
73
74
75
76
`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.

77
<!--- Intersection with cylinder, internal vs OCC representation --->
Dennis Gläser's avatar
Dennis Gläser committed
78
79
80
81
82
83
84
85
86
87
<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
88
89
90
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
91
92

Moreover, we want to reject all those quadrilaterals of which only
Dennis Gläser's avatar
Dennis Gläser committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106
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
107
into an instance of the class `ContainedEntityNetwork`, using the corresponding builder class:
Dennis Gläser's avatar
Dennis Gläser committed
108
109

```cpp
110
ContainedEntityNetworkBuilder<ctype> builder;
Dennis Gläser's avatar
Dennis Gläser committed
111
112
113
114
115
116
117

// 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
118
on sub-domains, and which entities are embedded in which sub-domain. Each sub-domain
Dennis Gläser's avatar
Dennis Gläser committed
119
receives a unique identifier by means of the `Id` class. By calling
120
121
122
`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
123
one could call `builder.addSubDomain(domain, Id(1))`, in which case embedded networks
124
will not be confined (see [example 3][2]). Entities are associated with the sub-
Dennis Gläser's avatar
Dennis Gläser committed
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
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.

148
149
# Enforcing constraints on confined entities

Dennis Gläser's avatar
Dennis Gläser committed
150
As mentioned before, enforcing the constraints between the raw entities, small length
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
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,
thus, using the final geometry of an entity after its confinement. This could be achieved
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;
173
const auto containedShape = intersect(getShape(quad), getShape(domain), /*epsilon*/1e-6);
174
175
176
const auto containedFaces = getFaces(containedShape);

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

180
if (!constraintsOnOther.evaluate(otherEntitySet, containedFaces))
181
182
183
184
{ status.increaseRejectedCounter(); continue; }
```

Note that the variable `containedFaces` is of type `std::vector<TopoDS_Face>`,
Dennis Gläser's avatar
Dennis Gläser committed
185
where `TopoDS_Face` is the OpenCascade type used to represent general two-dimensional
186
faces. In the above code snippets we have assumed `entitySet` and `otherEntitySet`
187
188
189
to also be of type `std::vector<TopoDS_Face>`, and to be composed of the previously
accepted faces.

190
191
192
193
194
195
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.

196
197
[go to example 3][2]

Dennis Gläser's avatar
Dennis Gläser committed
198
199
200
[0]: https://git.iws.uni-stuttgart.de/DennisGlaeser/frackit/tree/master/appl/example1
[1]: https://git.iws.uni-stuttgart.de/DennisGlaeser/frackit/tree/master/geometry/cylinder.hh
[2]: https://git.iws.uni-stuttgart.de/DennisGlaeser/frackit/tree/master/appl/example3