Skip to content

GitLab

Commit 68ee90dc authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[examples] add first example

parent ab7606ef
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ 68ee90dc
......@@ -110,5 +110,6 @@ include_directories( ${CMAKE_SOURCE_DIR} )
# include subdirectories
add_subdirectory(cmake)
add_subdirectory(doc)
add_subdirectory(examples)
add_subdirectory(frackit)
add_subdirectory(test)
examples/CMakeLists.txt 0 → 100644
View file @ 68ee90dc
add_subdirectory(layers)
examples/layers/CMakeLists.txt 0 → 100644
View file @ 68ee90dc
frackit_add_test(NAME example_layers
SOURCES example_layers.cc
COMPILE_DEFINITIONS BREPFILE="${CMAKE_SOURCE_DIR}/examples/layers/layers.brep"
LABELS examples)
set(CMAKE_BUILD_TYPE Debug)
examples/layers/example_layers.cc 0 → 100644
View file @ 68ee90dc
#include <iostream>
#include <fstream>
#include <string>
#include <stdexcept>
#include <random>
// Contains functionality to read in shapes
#include <BRepTools.hxx>
// OpenCascade shape class
#include <TopoDS_Shape.hxx>
// basic geometry types
#include <frackit/geometry/box.hh>
#include <frackit/geometry/disk.hh>
// headers containing functions to compute lengths/areas/volumes
#include <frackit/magnitude/magnitude.hh>
#include <frackit/magnitude/containedmagnitude.hh>
// parser functions between internal & brep data structures
#include <frackit/occ/breputilities.hh>
// sampler for points and disks
#include <frackit/sampling/pointsampling.hh>
#include <frackit/sampling/geometrysampling.hh>
// constraints to be enforced on the network (distance, angles, etc.)
#include <frackit/entitynetwork/constraints.hh>
// builder class for creating networks of entities confined in (sub-)domains
#include <frackit/entitynetwork/networkbuilder.hh>
// writes an entity network to a meshable Gmsh .geo file format
#include <frackit/io/gmshwriter.hh>
int main(int argc, char** argv)
{
//! print welcome message
std::cout << "\n\n"
<< "###################################################################\n"
<< "## Example: Creation of disk-shaped entities in a layered medium ##\n"
<< "###################################################################\n"
<< "\n\n";
/////////////////////////////////////////////////////
// 1. Read in the domain geometry from .brep file. //
// The file name is defined in CMakeLists.txt //
/////////////////////////////////////////////////////
using namespace Frackit;
const auto domainFileName = BREPFILE;
const auto domainShape = OCCUtilities::readShape(domainFileName);
// obtain the three solids contained in the file
const auto solids = OCCUtilities::getSolids(domainShape);
if (solids.size() != 3)
throw std::runtime_error("Unexpected result read from .brep file");
// the domain we want to create a network in is the center one
const auto& networkDomain = solids[1];
// compute domain volume and get its boundaries
const auto domainVolume = computeMagnitude(networkDomain);
const auto domainShells = OCCUtilities::getShells(networkDomain);
if (domainShells.size() > 1)
throw std::runtime_error(std::string("More than one shell read from solid"));
else if (domainShells.empty())
throw std::runtime_error(std::string("Could not read shell from solid"));
// get the boundary faces of our domain of intererst
const auto domainShellFaces = OCCUtilities::getFaces(domainShells[0]);
//////////////////////////////////////////////////////////////////////
// 2. Make sampler classes to randomly sample points (disk centers) //
// and disks (using the sampled center poitns) //
//////////////////////////////////////////////////////////////////////
using ctype = double;
using Disk = Disk<ctype>;
using DiskSampler = GeometrySampler<Disk>;
// sample points within bounding box of domain
const auto domainBBox = OCCUtilities::getBoundingBox(networkDomain);
// use default sampler: samples uniformly along each coordinate direction in the box
auto pointSampler = makeUniformPointSampler(domainBBox);
// mean major/minor axis length to be used in the disk samples
const ctype meanMajAxisLength = 20.0;
const ctype meanMinAxisLength = 15.0;
// sampler for disks of orientation 1
DiskSampler diskSampler1(std::normal_distribution<ctype>(meanMajAxisLength, 6.5), // major axis length: mean value & standard deviation
std::normal_distribution<ctype>(meanMinAxisLength, 4.5), // minor axis length: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)), // rotation around x-axis: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)), // rotation around y-axis: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5))); // rotation around z-axis: mean value & standard deviation
// sampler for disks of orientation 2
DiskSampler diskSampler2(std::normal_distribution<ctype>(meanMajAxisLength, 6.5), // major axis length: mean value & standard deviation
std::normal_distribution<ctype>(meanMinAxisLength, 4.5), // minor axis length: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(90.0), toRadians(7.5)), // rotation around x-axis: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5)), // rotation around y-axis: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(0.0), toRadians(7.5))); // rotation around z-axis: mean value & standard deviation
//! containers to store created entities
std::vector<Disk> diskSet1;
std::vector<Disk> diskSet2;
//! enforce some constraints on the network
auto constraintsOnSelf = makeDefaultConstraints<ctype>();
auto constraintsOnOther = makeDefaultConstraints<ctype>();
auto constraintsOnDomain = makeDefaultConstraints<ctype>();
// constraints among disks of the same set
constraintsOnSelf.setMinDistance(2.5);
constraintsOnSelf.setMinIntersectingAngle(toRadians(45.0));
constraintsOnSelf.setMinIntersectionMagnitude(2.5);
constraintsOnSelf.setMinIntersectionDistance(2.0);
// constraints with the other disk set
constraintsOnOther.setMinDistance(1.5);
constraintsOnOther.setMinIntersectingAngle(toRadians(10.0));
constraintsOnOther.setMinIntersectionMagnitude(2.5);
constraintsOnOther.setMinIntersectionDistance(2.0);
// constraints with the domain boundary
constraintsOnDomain.setMinDistance(1.5);
constraintsOnDomain.setMinIntersectingAngle(toRadians(5.0));
constraintsOnDomain.setMinIntersectionMagnitude(20.0);
constraintsOnDomain.setMinIntersectionDistance(2.0);
//! create pseudo-random number generator to select set during creation
std::default_random_engine randomNumber{std::random_device{}()};
//! keep track of number of created & accepted entities
std::size_t total = 0;
std::size_t accepted = 0;
std::size_t accepted_1 = 0;
std::size_t accepted_2 = 0;
//! Define how many entities should be created
const std::size_t numTargetEntities_1 = 4;
const std::size_t numTargetEntities_2 = 4;
//! print header of status output
std::cout << "##############################################################################################################################\n"
<< "# No. of entities | No. rejected entities | Acceptance Ratio | Current entity density [m²/m³] | Progress [%] #\n"
<< "#----------------------------------------------------------------------------------------------------------------------------#\n";
ctype currentDensity = 0.0;
ctype currentDiskArea = 0.0;
while (accepted_1 != numTargetEntities_1 || accepted_2 != numTargetEntities_2)
{
bool createSecondary = randomNumber()%2;
if (!createSecondary && accepted_1 == numTargetEntities_1)
createSecondary = true;
else if (createSecondary && accepted_2 == numTargetEntities_2)
createSecondary = false;
auto disk = createSecondary ? diskSampler1(pointSampler)
: diskSampler2(pointSampler);
total++;
// We don't want ellipses of too large aspect ratio
if (disk.majorAxisLength() > disk.minorAxisLength()*3.0) continue;
auto& diskSetSelf = createSecondary ? diskSet2 : diskSet1;
const auto& diskSetOther = createSecondary ? diskSet1 : diskSet2;
// enforce constraints w.r.t. to the own disk set
if (!constraintsOnSelf.evaluate(diskSetSelf, disk)) continue;
// enforce constraints w.r.t. to the other disk set
if (!constraintsOnOther.evaluate(diskSetOther, disk)) continue;
// enforce constraints w.r.t. the domain boundaries
if (!constraintsOnDomain.evaluate(domainShellFaces, disk)) continue;
// reject if disk area contained in the domain is too small
const auto containedArea = computeContainedMagnitude(disk, networkDomain);
if (containedArea < meanMajAxisLength*meanMinAxisLength*M_PI/5.0) continue;
// if we get here, the disk is admissible
diskSetSelf.push_back(std::move(disk));
accepted++;
if (createSecondary) accepted_2++;
else accepted_1++;
// compute new density (use minimum w.r.t. domain/network volume)
currentDiskArea += containedArea;
currentDensity = currentDiskArea/domainVolume;
std::cout << " " << std::setw(18) << std::setfill(' ')
<< std::to_string(diskSet1.size() + diskSet2.size()) + std::string(9, ' ')
<< "| " << std::setprecision(6) << std::setw(24) << std::setfill(' ')
<< std::to_string(total-accepted) + std::string(12, ' ')
<< "| " << std::setprecision(6) << std::setw(19) << std::setfill(' ')
<< std::to_string( double(double(accepted)/double(total)) ) + std::string(7, ' ')
<< "| " << std::setprecision(6) << std::setw(33) << std::setfill(' ')
<< std::to_string(currentDensity) + std::string(17, ' ')
<< "| " << std::string(3, ' ') + std::to_string( double(accepted_1+accepted_2)/double(numTargetEntities_1+numTargetEntities_2) ) << std::endl;
}
//////////////////////////////////////////////////////////////////////////
// 3. The entities of the network have been created. We can now //
// construct different types of networks (contained, confined, etc.) //
//////////////////////////////////////////////////////////////////////////
std::cout << "Constructing contained, confined network" << std::endl;
// 3.1 Construct a network such that the entities are confined in the
// sub-domain and which contains information about the containing domain.
ContainedEntityNetworkBuilder containedConfinedBuilder;
// define sub-domains
containedConfinedBuilder.addSubDomain(solids[0], /*subDomainIndex*/1);
containedConfinedBuilder.addSubDomain(networkDomain, /*subDomainIndex*/2);
containedConfinedBuilder.addSubDomain(solids[2], /*subDomainIndex*/3);
// define entity network for sub-domain 2
containedConfinedBuilder.addSubDomainEntities(diskSet1, 2);
containedConfinedBuilder.addSubDomainEntities(diskSet2, 2);
// build network
const auto containedConfinedNetwork = containedConfinedBuilder.build();
std::cout << "Constructing contained, unconfined network" << std::endl;
// 3.2 Construct a network such that the entities are NOT confined in the
// sub-domain and which contains information about the containing domain.
ContainedEntityNetworkBuilder containedUnconfinedBuilder;
// define sub-domains
containedUnconfinedBuilder.addSubDomain(solids[0], /*subDomainIndex*/1);
containedUnconfinedBuilder.addSubDomain(solids[2], /*subDomainIndex*/3);
containedUnconfinedBuilder.addSubDomain(networkDomain,
/*subDomainIndex*/2,
/*confineEmbeddedNetwork*/false);
// define entity network for sub-domain 2
containedUnconfinedBuilder.addSubDomainEntities(diskSet1, 2);
containedUnconfinedBuilder.addSubDomainEntities(diskSet2, 2);
// build network
const auto containedUnconfinedNetwork = containedUnconfinedBuilder.build();
std::cout << "Constructing uncontained, confined network" << std::endl;
// 3.3 Construct a network such that the entities are confined in the
// sub-domain, but, which does not contain information about the domains.
// This is computationally more efficient if only the 2d network is of interest.
// Moreover, the files written to disk are smaller as the domain is not contained.
EntityNetworkBuilder uncontainedConfinedBuilder;
// define sub-domains
uncontainedConfinedBuilder.addSubDomain(networkDomain,
/*subDomainIndex*/1,
/*confineEmbeddedNetwork*/true);
// define entity network for sub-domain 2
uncontainedConfinedBuilder.addSubDomainEntities(diskSet1, 1);
uncontainedConfinedBuilder.addSubDomainEntities(diskSet2, 1);
// build network
const auto uncontainedConfinedNetwork = uncontainedConfinedBuilder.build();
std::cout << "Constructing uncontained, unconfined network" << std::endl;
// 3.4 Construct the network only, independent of any domains
EntityNetworkBuilder uncontainedUnconfinedBuilder;
uncontainedUnconfinedBuilder.addEntities(diskSet1);
uncontainedUnconfinedBuilder.addEntities(diskSet2);
// build network
const auto uncontainedUnconfinedNetwork = uncontainedUnconfinedBuilder.build();
///////////////////////////////////////////////////////
// 4. Write the networks to disk in Gmsh .geo format //
///////////////////////////////////////////////////////
std::cout << "Write networks to disk" << std::endl;
GmshWriter containedConfinedWriter(containedConfinedNetwork);
containedConfinedWriter.write("contained_confined", // body of the filename to be used (will add .geo)
2.5, // mesh size to be used on entities
5.0); // mesh size to be used on domain boundaries
GmshWriter containedUnconfinedWriter(containedUnconfinedNetwork);
containedUnconfinedWriter.write("contained_unconfined", 2.5, 3.0);
GmshWriter uncontainedConfinedWriter(uncontainedConfinedNetwork);
uncontainedConfinedWriter.write("uncontained_confined", 2.5);
GmshWriter uncontainedUnconfinedWriter(uncontainedUnconfinedNetwork);
uncontainedUnconfinedWriter.write("uncontained_unconfined", 2.5);
return 0;
}
examples/layers/layers.brep 0 → 100644
View file @ 68ee90dc
DBRep_DrawableShape
CASCADE Topology V1, (c) Matra-Datavision
Locations 10
1
1 0 0 0
0 1 0 100
0 0 1 0
1
1 0 0 0
0 1 0 0
0 0 1 25
2 1 1 2 1 0
1
1 0 0 0
0 1 0 100
0 0 1 0
1
1 0 0 0
0 1 0 100
0 0 1 0
2 1 -1 0
2 2 -1 0
2 2 -1 1 -1 0
2 4 -1 0
2 5 -1 0
Curve2ds 20
1 0 0 1 0
1 0 -100 1 0
1 0 0 1 0
1 0 -25 1 0
1 0 0 1 0
1 0 -25 1 0
1 0 0 0 -1
1 1 0 0 -1
1 0 0 0 -1
1 0 0 0 -1
1 1 0 0 -1
1 1 0 0 -1
1 0 0 1 0
1 0 -100 1 0
1 0 0 0 -1
1 1 0 0 -1
1 0 0 1 0
1 0 -100 1 0
1 0 0 0 -1
1 1 0 0 -1
Curves 21
6 0 3 0 0 0 50 0 25 100 0 0 150 0 0
1 0 0 0 0 1 0
1 150 0 0 0 1 0
1 0 0 0 0 0 1
1 0 100 0 0 0 1
1 150 0 0 0 0 1
1 150 100 0 0 0 1
6 0 3 0 0 55 50 0 55 120 0 40 150 0 55
1 0 0 55 0 1 0
1 150 0 55 0 1 0
1 0 0 55 0 0 -1
1 150 0 55 0 0 -1
1 150 100 55 0 0 -1
1 0 100 55 0 0 -1
6 0 3 0 0 75 50 0 75 120 0 60 150 0 65
1 0 0 75 0 1 0
1 150 0 65 0 1 0
1 0 0 75 0 0 -1
1 150 0 65 0 0 -1
1 150 100 65 0 0 -1
1 0 100 75 0 0 -1
Polygon3D 0
PolygonOnTriangulations 0
Surfaces 15
6 -0 -1 -0
6 0 3 0 0 0 50 0 25 100 0 0 150 0 0
6 -0 -0 -1
6 0 3 0 0 0 50 0 25 100 0 0 150 0 0
6 -0 -0 -1
6 0 3 0 100 0 50 100 25 100 100 0 150 100 0
6 -0 -1 -0
6 0 3 0 0 55 50 0 55 120 0 40 150 0 55
6 -0 -1 -0
6 0 3 0 0 75 50 0 75 120 0 60 150 0 65
1 0 0 0 -1 0 0 0 1 0 0 0 -1
1 150 0 0 -1 0 0 0 1 0 0 0 -1
1 0 50 40 1 0 0 -0 0 1 0 -1 0
1 74.822958864791801 0 41.350284355050981 -0 -1 -0 0 0 -1 1 -0 0
1 150 50 40 -1 -0 -0 0 0 -1 0 -1 0
1 74.822958864791815 100 41.350284355050981 0 1 0 0 -0 1 1 0 -0
1 0 50 65 -1 -0 -0 0 0 -1 0 -1 0
1 73.573477252097632 0 59.699899660234827 -0 -1 -0 0 0 -1 1 -0 0
1 150 50 60 1 0 0 -0 0 1 0 -1 0
1 73.573477252097632 100.00000000000001 59.699899660234827 -2.0170901931128908e-31 1 -1.8564303712652899e-30 0 1.8564303712652899e-30 1 1 2.0170901931128908e-31 -3.7445874960761392e-61
Triangulations 0
TShapes 78
Ve
1e-07
0 0 0
0 0
0101101
*
Ve
1e-07
50 0 25
0 0
0101101
*
Ve
1e-07
100 0 0
0 0
0101101
*
Ve
1e-07
150 0 0
0 0
0101101
*
Ve
1e-07
0 0 55
0 0
0101101
*
Ve
1e-07
50 0 55
0 0
0101101
*
Ve
1e-07
120 0 40
0 0
0101101
*
Ve
1e-07
150 0 55
0 0
0101101
*
Ve
1e-07
0 0 75
0 0
0101101
*
Ve
1e-07
50 0 75
0 0
0101101
*
Ve
1e-07
120 0 60
0 0
0101101
*
Ve
1e-07
150 0 65
0 0
0101101
*
Ed
1e-07 1 1 0
1 1 0 0 1
2 1 1 0 0 1
2 2 1 6 0 1
2 3 2 0 0 1
2 4 2 7 0 1
2 5 3 6 0 1
2 6 3 8 0 1
0
0101000
+78 0 -75 0 *
Ed
1e-07 1 1 0
1 2 0 0 100
2 7 1 0 0 100
0
0101000
+78 0 -78 1 *
Ed
1e-07 1 1 0
1 3 0 0 100
2 8 1 0 0 100
0
0101000
+75 0 -75 1 *
Ed
1e-07 1 1 0
1 4 0 0 25
2 9 2 0 0 25
0
0101000
+78 0 -78 2 *
Ed
1e-07 1 1 0
1 5 0 0 25
2 10 3 0 0 25
0
0101000
+78 1 -78 3 *
Ed
1e-07 1 1 0
1 6 0 0 25
2 11 2 0 0 25
0
0101000
+75 0 -75 2 *
Ed
1e-07 1 1 0
1 7 0 0 25
2 12 3 0 0 25
0
0101000
+75 1 -75 3 *
Ed
1e-07 1 1 0
1 8 0 0 1
2 13 4 0 0 1
2 14 4 9 0 1
0
0101000
+74 0 -71 0 *
Ed
1e-07 1 1 0
1 9 0 0 100
2 15 4 0 0 100
0
0101000
+74 0 -74 4 *
Ed
1e-07 1 1 0
1 10 0 0 100
2 16 4 0 0 100
0
0101000
+71 0 -71 4 *
Ed
1e-07 1 1 0
1 11 0 0 30
0
0101000
+74 0 -78 2 *
Ed
1e-07 1 1 0
1 12 0 0 30
0
0101000
+71 0 -75 2 *
Ed
1e-07 1 1 0
1 13 0 0 30
0
0101000
+71 4 -75 3 *
Ed
1e-07 1 1 0
1 14 0 0 30
0
0101000
+74 4 -78 3 *
Ed
1e-07 1 1 0
1 15 0 0 1
2 17 5 0 0 1
2 18 5 10 0 1
0
0101000
+70 0 -67 0 *
Ed
1e-07 1 1 0
1 16 0 0 100
2 19 5 0 0 100
0
0101000
+70 0 -70 5 *
Ed
1e-07 1 1 0
1 17 0 0 100
2 20 5 0 0 100
0
0101000