Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
tools
frackit
Commits
3ed481ef
Commit
3ed481ef
authored
Dec 18, 2020
by
Dennis Gläser
Browse files
[ex4] modify example code
parent
7a9a6ae1
Pipeline
#2845
canceled with stages
in 6 seconds
Changes
3
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
appl/example4/README.md
View file @
3ed481ef
...
...
@@ -14,7 +14,7 @@ realize this example using the Frackit python bindings.__
*
create a domain geometry by intersection/cut operations on basic geometry types
*
generate random polygons using the
`PolygonSampler`
class
*
design
your own
incremental network generation algorithm
*
design
a custom,
incremental network generation algorithm
Let us consider some idealized, hexahedral geological layer with an embedded
tunnel structure. The excavation of tunnels can have an effect on the stress
...
...
@@ -24,7 +24,7 @@ of how such situations can modeled using `Frackit`.
### Step 1: domain generation
We want to focus on a spherical region around the center of
the
hexahedral
We want to focus on a spherical region around the center of
a
hexahedral
layer. To this end, we construct instances of the
`Box`
and
`Sphere`
classes
to represent these two geometries:
...
...
@@ -71,7 +71,7 @@ if (solids.size() != 1) throw std::runtime_error("Cut operation with tunnel shou
const
auto
&
domain
=
solids
[
0
];
```
## Step 2: entity sampler definitions
##
#
Step 2: entity sampler definitions
Instead of randomly sample entities within the domain volume, in this example we
want to use a different approach. Here, the goal is to generate a network in the
...
...
@@ -196,45 +196,32 @@ rejected immediately.
__Stage 2__
: For each entity
`e`
of stage 1, we want to generate an entity of orientation 2
that intersects
`e`
. Therefore, we loop over all entities of stage 1 and sample new
candidates in their vicinity. This is achieved by computing the bounding box of
`e`
,
which is then used to define a region in which new candidates are sampled:
which is then used to define a region in which new candidates are sampled. This
is done within the
`getVicinityPointSampler`
lambda:
```
cpp
for
(
const
auto
&
e
:
entitiesSet1
)
// this lambda function generates a point sampler within the vicinity of the given geometry
auto
getVicinityPointSampler
=
[
&
]
(
const
auto
&
g
,
unsigned
int
orientationIdx
)
{
using
std
::
sqrt
;
// sample an entity of orientation 2 in the vicinity of e
const
auto
rawBBox
=
getBoundingBox
(
e
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
e
));
// characteristic length of e
const
auto
rawBBox
=
getBoundingBox
(
g
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
g
));
const
Box
sampleBox
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
auto
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
/*orientationIdx*/
2
);
return
getSampler
(
makeUniformPointSampler
(
sampleBox
),
orientationIdx
);
};
```
The variable
`sampler`
now hold
s an instance of the
`PolygonSampler`
class, which generates
polygons of orientation
2
within the box
`sampleBox`
that describes the vicinity of
`e
`
.
We then sample candidates until we have found an entity that fulfills all desired constraints.
After that, it is continued with the next
entity of
`
ent
i
ti
esSet1`
until all have been visited
.
Reusing the
`getSampler`
lambda, this return
s an instance of the
`PolygonSampler`
class, which generates
polygons of
the given
orientation within the box
`sampleBox`
that describes the vicinity of
the given geometry
`g
`
.
Such samplers are then successively created for each entity of stage 1 in order to
find one intersecting
entity of
ori
ent
a
ti
on 2
.
__Stage 3__
: this stage follows the same logic as stage 2, but it samples entities of
orientation 3 and looks for entity candidates in the vicinity of the entities of stage 2:
```
cpp
for
(
const
auto
&
e
:
entitiesSet2
)
{
using
std
::
sqrt
;
// sample an entity of orientation 2 in the vicinity of e
const
auto
rawBBox
=
getBoundingBox
(
e
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
e
));
const
Box
sampleBox
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
auto
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
/*orientationIdx*/
3
);
```
Note that due to the restrictive constraints used in this example, the entity
generation step can take up to a few minutes.
...
...
@@ -248,8 +235,10 @@ to construct the network and embed it in our domain:
std
::
cout
<<
"Building and writing network"
<<
std
::
endl
;
ContainedEntityNetworkBuilder
<
ctype
>
builder
;
// add sub-domains
builder
.
addConfiningSubDomain
(
domain
,
Id
(
1
));
// add sub-domain
builder
.
addConfiningSubDomain
(
domain
,
Id
(
1
));
// embed all entities in that domain
builder
.
addSubDomainEntities
(
entitiesSet1
,
Id
(
1
));
builder
.
addSubDomainEntities
(
entitiesSet2
,
Id
(
1
));
builder
.
addSubDomainEntities
(
entitiesSet3
,
Id
(
1
));
...
...
appl/example4/example4.cc
View file @
3ed481ef
...
...
@@ -89,23 +89,25 @@ int main(int argc, char** argv)
const
Circle
tunnelBase
(
Point
(
-
25.0
,
0.0
,
0.0
),
tunnelDirection
,
tunnelRadius
);
const
Cylinder
tunnel
(
tunnelBase
,
/*height*/
50.0
);
// the domain
i
s the layer
constrained to the domainSphere and cut by the tunnel
//
restrict
the domain s
phere to
the layer
(this should yield a single solid shape)
const
auto
layerSphereCut
=
OCCUtilities
::
intersect
(
OCCUtilities
::
getShape
(
layer
),
OCCUtilities
::
getShape
(
domainSphere
),
/*eps*/
1e-6
);
auto
solids
=
OCCUtilities
::
getSolids
(
layerSphereCut
);
if
(
solids
.
size
()
!=
1
)
throw
std
::
runtime_error
(
"Intersection operation with domainSphere should yield a single solid"
);
// now cut the tunnel out of the domain (should also yield a single solid)
const
auto
tunnelCut
=
OCCUtilities
::
cut
(
solids
[
0
],
OCCUtilities
::
getShape
(
tunnel
),
/*eps*/
1e-6
);
solids
=
OCCUtilities
::
getSolids
(
tunnelCut
);
if
(
solids
.
size
()
!=
1
)
throw
std
::
runtime_error
(
"Cut operation with tunnel should yield a single solid"
);
//
get the
domain shape
//
we now have the resulting
domain shape
const
auto
&
domain
=
solids
[
0
];
// get the OCC reprensentation of the sphere surface
const
auto
sphereFaces
=
OCCUtilities
::
getFaces
(
OCCUtilities
::
getShape
(
domainSphere
));
if
(
sphereFaces
.
size
()
!=
1
)
throw
std
::
runtime_error
(
"We expect the sphere to consist of a single face"
);
// get the OCC representation of the sphere surface (for constraints evaluation)
const
auto
sphereShape
=
OCCUtilities
::
getShape
(
domainSphere
);
const
auto
&
sphereFaces
=
OCCUtilities
::
getFaces
(
sphereShape
);
if
(
sphereFaces
.
size
()
!=
1
)
throw
std
::
runtime_error
(
"Sphere should contain a single face"
);
const
auto
&
sphereSurfaceShape
=
sphereFaces
[
0
];
// get the OCC representation of tunnel surface
...
...
@@ -188,7 +190,7 @@ int main(int argc, char** argv)
auto
producesSmallLengthScale
=
[]
(
const
auto
&
geometry
,
const
auto
&
is
)
->
bool
{
// lambda to check if the intersection is too close (< 1mm) to a vertex of the geometry
auto
isTooClose
=
[
&
is
]
(
const
auto
&
v
)
{
return
computeDistance
(
v
,
is
)
<
1e-
3
;
};
auto
isTooClose
=
[
&
is
]
(
const
auto
&
v
)
{
return
computeDistance
(
v
,
is
)
<
1e-
2
;
};
const
auto
shape
=
OCCUtilities
::
getShape
(
geometry
);
const
auto
vertices
=
OCCUtilities
::
getVertices
(
shape
);
...
...
@@ -237,7 +239,8 @@ int main(int argc, char** argv)
{
status
.
increaseRejectedCounter
();
continue
;
}
// compute the part contained in the domain
const
auto
containedShape
=
OCCUtilities
::
intersect
(
OCCUtilities
::
getShape
(
candidate
),
domain
,
/*eps*/
1e-6
);
const
auto
candidateShape
=
OCCUtilities
::
getShape
(
candidate
);
const
auto
containedShape
=
OCCUtilities
::
intersect
(
candidateShape
,
domain
,
/*eps*/
1e-6
);
const
auto
containedFaces
=
OCCUtilities
::
getFaces
(
containedShape
);
// skip everything that does not lead to a single contained face
...
...
@@ -276,22 +279,25 @@ int main(int argc, char** argv)
<<
maxTriesForEntity
<<
" unsuccessful tries. --"
<<
std
::
endl
;
};
status
.
reset
();
status
.
setTargetCount
(
idSet2
,
entitiesSet1
.
size
());
for
(
const
auto
&
e
:
entitiesSet1
)
// this lambda function generates a point sampler within the vicinity of the given geometry
auto
getVicinityPointSampler
=
[
&
]
(
const
auto
&
g
,
unsigned
int
orientationIdx
)
{
using
std
::
sqrt
;
// sample an entity of orientation 2 in the vicinity of e
const
auto
rawBBox
=
getBoundingBox
(
e
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
e
));
const
auto
rawBBox
=
getBoundingBox
(
g
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
g
));
const
Box
sampleBox
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
return
getSampler
(
makeUniformPointSampler
(
sampleBox
),
orientationIdx
);
};
status
.
reset
();
status
.
setTargetCount
(
idSet2
,
entitiesSet1
.
size
());
for
(
const
auto
&
e
:
entitiesSet1
)
{
bool
accepted
=
false
;
std
::
size_t
tryCount
=
0
;
auto
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
/*orientationIdx*/
2
);
auto
sampler
=
getVicinityPointSampler
(
e
,
/*orientationIdx*/
2
);
while
(
!
accepted
&&
tryCount
<=
maxTriesForEntity
)
{
tryCount
++
;
...
...
@@ -309,7 +315,8 @@ int main(int argc, char** argv)
{
status
.
increaseRejectedCounter
();
continue
;
}
// compute the part contained in the domain
const
auto
containedShape
=
OCCUtilities
::
intersect
(
OCCUtilities
::
getShape
(
candidate
),
domain
,
/*eps*/
1e-6
);
const
auto
candidateShape
=
OCCUtilities
::
getShape
(
candidate
);
const
auto
containedShape
=
OCCUtilities
::
intersect
(
candidateShape
,
domain
,
/*eps*/
1e-6
);
const
auto
containedFaces
=
OCCUtilities
::
getFaces
(
containedShape
);
// skip everything that does not lead to a single contained face
...
...
@@ -350,18 +357,9 @@ int main(int argc, char** argv)
status
.
setTargetCount
(
idSet3
,
entitiesSet2
.
size
());
for
(
const
auto
&
e
:
entitiesSet2
)
{
using
std
::
sqrt
;
// sample an entity of orientation 2 in the vicinity of e
const
auto
rawBBox
=
getBoundingBox
(
e
);
const
auto
charLength
=
sqrt
(
computeMagnitude
(
e
));
const
Box
sampleBox
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
auto
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
/*orientationIdx*/
3
);
bool
accepted
=
false
;
std
::
size_t
tryCount
=
0
;
auto
sampler
=
getVicinityPointSampler
(
e
,
/*orientationIdx*/
3
);
while
(
!
accepted
&&
tryCount
<=
maxTriesForEntity
)
{
tryCount
++
;
...
...
@@ -379,7 +377,8 @@ int main(int argc, char** argv)
{
status
.
increaseRejectedCounter
();
continue
;
}
// compute the part contained in the domain
const
auto
containedShape
=
OCCUtilities
::
intersect
(
OCCUtilities
::
getShape
(
candidate
),
domain
,
/*eps*/
1e-6
);
const
auto
candidateShape
=
OCCUtilities
::
getShape
(
candidate
);
const
auto
containedShape
=
OCCUtilities
::
intersect
(
candidateShape
,
domain
,
/*eps*/
1e-6
);
const
auto
containedFaces
=
OCCUtilities
::
getFaces
(
containedShape
);
// skip everything that does not lead to a single contained face
...
...
@@ -431,8 +430,10 @@ int main(int argc, char** argv)
std
::
cout
<<
"Building and writing network"
<<
std
::
endl
;
ContainedEntityNetworkBuilder
<
ctype
>
builder
;
// add sub-domains
builder
.
addConfiningSubDomain
(
domain
,
Id
(
1
));
// add sub-domain
builder
.
addConfiningSubDomain
(
domain
,
Id
(
1
));
// embed all entities in that domain
builder
.
addSubDomainEntities
(
entitiesSet1
,
Id
(
1
));
builder
.
addSubDomainEntities
(
entitiesSet2
,
Id
(
1
));
builder
.
addSubDomainEntities
(
entitiesSet3
,
Id
(
1
));
...
...
appl/example4/example4.py
View file @
3ed481ef
import
sys
import
math
from
time
import
process_time
startTime
=
process_time
()
...
...
@@ -20,21 +21,25 @@ print("\n\n"
from
frackit.geometry
import
Box
,
Sphere
,
Cylinder
,
Direction
,
Vector
,
Circle
,
Point
# the layer has a thickness of 40m, and initially we consider 100m lateral extent
layer
=
Box
(
-
50.0
,
-
50.0
,
-
20.0
,
50.0
,
50.0
,
20.0
)
domainSphere
=
Sphere
(
Point
(
0.0
,
0.0
,
0.0
),
25.0
)
layer
=
Box
(
xmin
=-
50.0
,
ymin
=-
50.0
,
zmin
=-
20.0
,
xmax
=
50.0
,
ymax
=
50.0
,
zmax
=
20.0
)
domainSphere
=
Sphere
(
center
=
Point
(
0.0
,
0.0
,
0.0
),
radius
=
25.0
)
# a tunnel with a radius of 5m crosses the domain
tunnelRadius
=
5.0
tunnelDirection
=
Direction
(
Vector
(
1.0
,
0.0
,
0.0
))
tunnelBase
=
Circle
(
Point
(
-
25.0
,
0.0
,
0.0
),
tunnelDirection
,
tunnelRadius
)
tunnel
=
Cylinder
(
tunnelBase
,
50.0
)
tunnelBase
=
Circle
(
center
=
Point
(
-
25.0
,
0.0
,
0.0
),
normal
=
tunnelDirection
,
radius
=
tunnelRadius
)
tunnel
=
Cylinder
(
bottomCircle
=
tunnelBase
,
height
=
50.0
)
# the domain
is the layer constrained to the domainSphere and cut by the tunnel
#
restrict
the domain
sphere to the layer (should yield a single solid)
from
frackit.occutilities
import
getShape
,
getSolids
,
getFaces
,
cut
from
frackit.occutilities
import
intersect
as
intersectShapes
solids
=
getSolids
(
intersectShapes
(
getShape
(
layer
),
getShape
(
domainSphere
),
1e-6
)
)
if
len
(
solids
)
!=
1
:
sys
.
exit
(
"Intersection operation with domainSphere should yield a single solid"
)
# cut the tunnel out of the domain
solids
=
getSolids
(
cut
(
solids
[
0
],
getShape
(
tunnel
),
1e-6
)
)
if
len
(
solids
)
!=
1
:
sys
.
exit
(
"Cut operation with tunnel should yield a single solid"
)
...
...
@@ -66,28 +71,30 @@ def choiceSampler(a, b):
return
sample
# function to create polygon samplers for the different orientations
from
frackit.common
import
toRadians
from
frackit.sampling
import
PolygonSampler
sizeDistro
=
uniformSampler
(
3.0
,
6.0
)
numCornersDistro
=
choiceSampler
(
4
,
9
)
def
getSampler
(
pointSampler
,
orientation
):
if
orientation
==
1
:
return
PolygonSampler
(
pointSampler
,
# sampler for polygon centers
normalSampler
(
toRadians
(
15.0
),
toRadians
(
10.0
)),
# strike angle: mean value & standard deviation
normalSampler
(
toRadians
(
90.0
),
toRadians
(
10.0
)),
# dip angle: mean value & standard deviation
sizeDistro
,
sizeDistro
,
numCornersDistro
)
# strike & dip length, number of corners
return
PolygonSampler
(
pointSampler
=
pointSampler
,
strikeAngleSampler
=
normalSampler
(
math
.
radians
(
15.0
),
math
.
radians
(
10.0
)),
dipAngleSampler
=
normalSampler
(
math
.
radians
(
90.0
),
math
.
radians
(
10.0
)),
strikeLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
dipLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
numCornersSampler
=
choiceSampler
(
4
,
9
))
elif
orientation
==
2
:
return
PolygonSampler
(
pointSampler
,
# sampler for polygon centers
normalSampler
(
toRadians
(
0.0
),
toRadians
(
10.0
)),
# strike angle: mean value & standard deviation
normalSampler
(
toRadians
(
0.0
),
toRadians
(
10.0
)),
# dip angle: mean value & standard deviation
sizeDistro
,
sizeDistro
,
numCornersDistro
)
# strike & dip length, number of corners
return
PolygonSampler
(
pointSampler
=
pointSampler
,
strikeAngleSampler
=
normalSampler
(
math
.
radians
(
0.0
),
math
.
radians
(
10.0
)),
dipAngleSampler
=
normalSampler
(
math
.
radians
(
0.0
),
math
.
radians
(
10.0
)),
strikeLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
dipLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
numCornersSampler
=
choiceSampler
(
4
,
9
))
elif
orientation
==
3
:
return
PolygonSampler
(
pointSampler
,
# sampler for polygon centers
normalSampler
(
toRadians
(
-
75.0
),
toRadians
(
10.0
)),
# strike angle: mean value & standard deviation
normalSampler
(
toRadians
(
90.0
),
toRadians
(
10.0
)),
# dip angle: mean value & standard deviation
sizeDistro
,
sizeDistro
,
numCornersDistro
)
# strike & dip length, number of corners
return
PolygonSampler
(
pointSampler
=
pointSampler
,
strikeAngleSampler
=
normalSampler
(
math
.
radians
(
-
75.0
),
math
.
radians
(
10.0
)),
dipAngleSampler
=
normalSampler
(
math
.
radians
(
90.0
),
math
.
radians
(
10.0
)),
strikeLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
dipLengthSampler
=
uniformSampler
(
3.0
,
6.0
),
numCornersSampler
=
choiceSampler
(
4
,
9
))
else
:
import
sys
sys
.
exit
(
"Unsupported orientation index"
)
...
...
@@ -104,7 +111,7 @@ from frackit.entitynetwork import EntityNetworkConstraints, EntityNetworkConstra
def
makeDefaultConstraints
():
c
=
EntityNetworkConstraints
()
c
.
setMinDistance
(
1.0
)
c
.
setMinIntersectingAngle
(
toR
adians
(
25.0
))
c
.
setMinIntersectingAngle
(
math
.
r
adians
(
25.0
))
c
.
setMinIntersectionMagnitude
(
0.1
)
c
.
setMinIntersectionDistance
(
0.1
)
return
c
...
...
@@ -124,21 +131,20 @@ constraints3.setMinDistance(0.15)
# Define constraints between entities of different sets
constraintsOnOther
=
makeDefaultConstraints
()
constraintsOnOther
.
setMinDistance
(
0.25
)
constraintsOnOther
.
setMinIntersectingAngle
(
toR
adians
(
30.0
))
constraintsOnOther
.
setMinIntersectingAngle
(
math
.
r
adians
(
30.0
))
# Moreover, we define constraints w.r.t. the domain boundary
constraintsOnBoundary
=
makeDefaultConstraints
()
constraintsOnBoundary
.
setMinIntersectingAngle
(
toR
adians
(
10.0
))
constraintsOnBoundary
.
setMinIntersectingAngle
(
math
.
r
adians
(
10.0
))
# When polygons intersect, the intersection might be very close to
# one of the corner points. We want to avoid small length scales caused
# by this, and this
lambda
provides the check for it.
# by this, and this
function
provides the check for it.
def
producesSmallLengthScale
(
geometry
,
isection
):
from
frackit.geometry
import
computeDistance
from
frackit.occutilities
import
getVertices
verts
=
getVertices
(
getShape
(
geometry
))
if
any
((
computeDistance
(
v
,
isection
)
<
1e-
3
for
v
in
verts
)):
return
True
if
any
((
computeDistance
(
v
,
isection
)
<
1e-
2
for
v
in
verts
)):
return
True
return
False
###########################
...
...
@@ -158,8 +164,11 @@ from frackit.geometry import HollowCylinder
numTargetEntities
=
75
status
=
SamplingStatus
()
status
.
setTargetCount
(
setId1
,
numTargetEntities
)
sampleCylinder1
=
HollowCylinder
(
tunnel
.
bottomFace
().
center
(),
tunnel
.
direction
(),
tunnel
.
radius
(),
tunnel
.
radius
()
+
2.0
,
tunnel
.
height
())
sampleCylinder1
=
HollowCylinder
(
bottomCenter
=
tunnel
.
bottomFace
().
center
(),
axis
=
tunnel
.
direction
(),
innerRadius
=
tunnel
.
radius
(),
outerRadius
=
tunnel
.
radius
()
+
2.0
,
height
=
tunnel
.
height
())
from
frackit.geometry
import
computeMagnitude
,
getBoundingBox
from
frackit.geometry
import
intersect
...
...
@@ -167,7 +176,7 @@ from frackit.occutilities import getFaces
from
frackit.sampling
import
makeUniformPointSampler
print
(
'## Step 1: generation of entities of orientation 1 around the tunnel
\n\n\n
'
)
area
=
0.0
;
area
=
0.0
polygonSampler1
=
getSampler
(
makeUniformPointSampler
(
sampleCylinder1
),
1
)
while
not
status
.
finished
():
candidate
=
polygonSampler1
()
...
...
@@ -182,15 +191,16 @@ while not status.finished():
status
.
increaseRejectedCounter
();
continue
;
# compute the part contained in the domain
containedShape
=
intersectShapes
(
getShape
(
candidate
),
domain
,
1e-6
)
candidateShape
=
getShape
(
candidate
)
containedShape
=
intersectShapes
(
candidateShape
,
domain
,
1e-6
)
containedFaces
=
getFaces
(
containedShape
)
# skip everything that does not lead to a single contained face
if
len
(
containedFaces
)
!=
1
:
status
.
increaseRejectedCounter
();
continue
;
# if too little of the candidate is inside the domain, reject it
containedFace
=
containedFaces
[
0
]
;
containedFaceArea
=
computeMagnitude
(
containedFace
)
;
containedFace
=
containedFaces
[
0
]
containedFaceArea
=
computeMagnitude
(
containedFace
)
if
containedFaceArea
<
0.25
*
candidate
.
area
():
status
.
increaseRejectedCounter
();
continue
;
...
...
@@ -210,26 +220,27 @@ while not status.finished():
print
(
'
\n\n\n
## Step 2: for each entity so far, try to generate an intersecting entity with orientation 2
\n\n\n
'
)
# after rejecting 500 candidates for an entity, we move to the next one
maxTriesForEntity
=
500
;
maxTriesForEntity
=
500
def
printSkipMessage
():
print
(
'
\t\t
-- Skipped search for this entity after '
+
\
str
(
maxTriesForEntity
)
+
' unsuccessful tries. --
\n
'
)
# function to generate a point sampler within the vicinity of the given geometry
def
getVicinityPointSampler
(
geometry
,
orientationIdx
):
rawBBox
=
getBoundingBox
(
geometry
)
charLength
=
math
.
sqrt
(
computeMagnitude
(
geometry
))
sampleBox
=
Box
(
xmin
=
rawBBox
.
xMin
()
-
charLength
,
ymin
=
rawBBox
.
yMin
()
-
charLength
,
zmin
=
rawBBox
.
zMin
()
-
charLength
,
xmax
=
rawBBox
.
xMax
()
+
charLength
,
ymax
=
rawBBox
.
yMax
()
+
charLength
,
zmax
=
rawBBox
.
zMax
()
+
charLength
)
return
getSampler
(
makeUniformPointSampler
(
sampleBox
),
orientationIdx
)
from
math
import
sqrt
status
.
reset
()
status
.
setTargetCount
(
setId2
,
len
(
entitySets
[
setId1
]))
for
e
in
entitySets
[
setId1
]:
# sample an entity of orientation 2 in the vicinity of e
rawBBox
=
getBoundingBox
(
e
);
charLength
=
sqrt
(
computeMagnitude
(
e
));
sampleBox
=
Box
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
tryCount
=
0
;
accepted
=
False
;
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
2
);
tryCount
=
0
accepted
=
False
sampler
=
getVicinityPointSampler
(
e
,
2
)
while
not
accepted
and
tryCount
<=
maxTriesForEntity
:
tryCount
+=
1
candidate
=
sampler
()
...
...
@@ -253,8 +264,8 @@ for e in entitySets[setId1]:
if
len
(
containedFaces
)
!=
1
:
status
.
increaseRejectedCounter
();
continue
;
# if too little of the candidate is inside the domain, reject it
containedFace
=
containedFaces
[
0
]
;
containedFaceArea
=
computeMagnitude
(
containedFace
)
;
containedFace
=
containedFaces
[
0
]
containedFaceArea
=
computeMagnitude
(
containedFace
)
if
containedFaceArea
<
0.25
*
candidate
.
area
():
status
.
increaseRejectedCounter
();
continue
;
...
...
@@ -268,7 +279,7 @@ for e in entitySets[setId1]:
if
not
constraintsOnBoundary
.
evaluate
(
tunnelSurfaceShape
,
containedFace
):
status
.
increaseRejectedCounter
();
continue
;
accepted
=
True
;
accepted
=
True
area
+=
containedFaceArea
entitySets
[
setId2
].
append
(
containedFace
)
status
.
increaseCounter
(
setId2
)
...
...
@@ -283,16 +294,9 @@ status.setTargetCount(setId3, len(entitySets[setId2]))
for
e
in
entitySets
[
setId2
]:
# sample an entity of orientation 2 in the vicinity of e
rawBBox
=
getBoundingBox
(
e
);
charLength
=
sqrt
(
computeMagnitude
(
e
));
sampleBox
=
Box
(
rawBBox
.
xMin
()
-
charLength
,
rawBBox
.
yMin
()
-
charLength
,
rawBBox
.
zMin
()
-
charLength
,
rawBBox
.
xMax
()
+
charLength
,
rawBBox
.
yMax
()
+
charLength
,
rawBBox
.
zMax
()
+
charLength
);
tryCount
=
0
;
accepted
=
False
;
sampler
=
getSampler
(
makeUniformPointSampler
(
sampleBox
),
3
);
tryCount
=
0
accepted
=
False
sampler
=
getVicinityPointSampler
(
e
,
3
)
while
not
accepted
and
tryCount
<=
maxTriesForEntity
:
tryCount
+=
1
candidate
=
sampler
()
...
...
@@ -316,8 +320,8 @@ for e in entitySets[setId2]:
if
len
(
containedFaces
)
!=
1
:
status
.
increaseRejectedCounter
();
continue
;
# if too little of the candidate is inside the domain, reject it
containedFace
=
containedFaces
[
0
]
;
containedFaceArea
=
computeMagnitude
(
containedFace
)
;
containedFace
=
containedFaces
[
0
]
containedFaceArea
=
computeMagnitude
(
containedFace
)
if
containedFaceArea
<
0.25
*
candidate
.
area
():
status
.
increaseRejectedCounter
();
continue
;
...
...
@@ -333,7 +337,7 @@ for e in entitySets[setId2]:
if
not
constraintsOnBoundary
.
evaluate
(
tunnelSurfaceShape
,
containedFace
):
status
.
increaseRejectedCounter
();
continue
;
accepted
=
True
;
accepted
=
True
area
+=
containedFaceArea
entitySets
[
setId3
].
append
(
containedFace
)
status
.
increaseCounter
(
setId3
)
...
...
@@ -342,7 +346,7 @@ for e in entitySets[setId2]:
if
tryCount
>=
maxTriesForEntity
:
printSkipMessage
()
# print the final entity density
domainVolume
=
computeMagnitude
(
domain
)
;
domainVolume
=
computeMagnitude
(
domain
)
numEntities
=
len
(
entitySets
[
setId1
])
+
len
(
entitySets
[
setId2
])
+
len
(
entitySets
[
setId3
])
print
(
"
\n
Created a network consisting of {:d} entities."
.
format
(
numEntities
))
print
(
"Volume of the domain: {:.2f}."
.
format
(
domainVolume
))
...
...
@@ -354,7 +358,7 @@ print("Area of the network: {:.2f} m², which corresponds to a density of {:.2f}
print
(
"Building and writing network"
)
from
frackit.entitynetwork
import
ContainedEntityNetworkBuilder
builder
=
ContainedEntityNetworkBuilder
()
;
builder
=
ContainedEntityNetworkBuilder
()
# add sub-domains
builder
.
addConfiningSubDomain
(
domain
,
Id
(
1
))
...
...
@@ -362,7 +366,7 @@ for id in entitySets: builder.addSubDomainEntities(entitySets[id], Id(1))
# now we can build and write out the network in Gmsh file format
from
frackit.io
import
GmshWriter
gmshWriter
=
GmshWriter
(
builder
.
build
())
;
gmshWriter
=
GmshWriter
(
builder
.
build
())
gmshWriter
.
setMeshSize
(
GmshWriter
.
GeometryTag
.
entity
,
0.5
)
gmshWriter
.
setMeshSize
(
GmshWriter
.
GeometryTag
.
subDomain
,
5.0
)
gmshWriter
.
write
(
"network"
)
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment