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

[draft] example5

parent e09964ad
......@@ -2,3 +2,4 @@ add_subdirectory(example1)
add_subdirectory(example2)
add_subdirectory(example3)
add_subdirectory(example4)
add_subdirectory(example5)
frackit_add_test(NAME example5
SOURCES example5.cc
LABELS example
COMMAND ./example5
CMD_ARGS 5)
frackit_symlink_or_copy(FILES example5.py)
#include <algorithm>
#include <iostream>
#include <stdexcept>
#include <random>
#include <vector>
#include <cmath>
#include <chrono>
#include <string>
// utility functions
#include <frackit/common/math.hh>
#include <frackit/common/id.hh>
// basic geometry types
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/cylinder.hh>
#include <frackit/geometry/boundingbox.hh>
#include <frackit/geometryutilities/getboundingbox.hh>
// header containing functions to compute contained lengths/areas/volumes
#include <frackit/magnitude/containedmagnitude.hh>
// sampler for points and disks
#include <frackit/sampling/makeuniformpointsampler.hh>
#include <frackit/sampling/disksampler.hh>
#include <frackit/sampling/status.hh>
// constraints to be enforced on the network (distance, angles, etc.)
#include <frackit/entitynetwork/constraints.hh>
// brep utilities of OpenCascade, for boolean operations on shapes
#include <frackit/occ/breputilities.hh>
// intersection detection for internal geometry types
#include <frackit/intersection/intersect.hh>
#include <frackit/intersection/emptyintersection.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)
{
// get start time to compute CPU time later
const auto start = std::chrono::steady_clock::now();
//! print welcome message
std::cout << "\n\n"
<< "##################################################################################\n"
<< "## Example 5: Creation of elliptical entities between injection/extraction wells #\n"
<< "##################################################################################\n"
<< "\n\n\n";
// Define some types used here
using namespace Frackit;
using ctype = double;
using Disk = Frackit::Disk<ctype>;
using BoundingBox = Frackit::BoundingBox<ctype>;
using Cylinder = Frackit::Cylinder<ctype>;
using Point = typename Disk::Point;
using Direction = typename Disk::Direction;
using Vector = typename Direction::Vector;
using Circle = typename Cylinder::Circle;
/////////////////////////////
// 1. Construct the domain //
/////////////////////////////
const ctype wellRadius = 0.25;
const ctype wellDepth = 15.0;
const ctype wellDistance = 35.0;
const ctype injectionLength = 2.5;
const ctype domainWidth = wellDistance*2.0;
const ctype domainDepth = wellDepth*2.0;
// define injection/extraction wells
const Direction wellNormal(Vector(0.0, 0.0, 1.0));
const Point injWellCenter(-wellDistance/2.0, 0.0, -wellDepth);
const Point extWellCenter( wellDistance/2.0, 0.0, -wellDepth);
const Circle injWellBase(injWellCenter, wellNormal, wellRadius);
const Circle extWellBase(extWellCenter, wellNormal, wellRadius);
const Cylinder injWell(injWellBase, wellDepth);
const Cylinder extWell(extWellBase, wellDepth);
// define domain
const Box rawDomain(-domainWidth/2.0, -domainWidth/2.0, -domainDepth,
domainWidth/2.0, domainWidth/2.0, 0.0);
using Frackit::OCCUtilities::getShape;
using Frackit::OCCUtilities::getSolids;
using Frackit::OCCUtilities::fuse;
using Frackit::OCCUtilities::cut;
const auto wellFusion = fuse({getShape(injWell), getShape(extWell)}, 1e-6);
const auto domainCut = cut(getShape(rawDomain), wellFusion, 1e-6);
const auto domainSolids = getSolids(domainCut);
if (domainSolids.size() != 1)
throw std::runtime_error("Expected domain to consist of a single solid");
const auto& domain = domainSolids[0];
//////////////////////////////////////////////////////////////////////////////
// 2. Make sampler classes to randomly sample points (entity centers) //
// and disk as entities (using the sampled center points) //
//////////////////////////////////////////////////////////////////////////////
// lambda to make disk sampler for given bounding box
const ctype diskCharLength = 3.0;
auto makeDiskSampler = [diskCharLength] (const auto& bbox)
{
using NormalDistro = std::normal_distribution<ctype>;
return DiskSampler(makeUniformPointSampler(bbox), // sampler for disk center points
NormalDistro(diskCharLength, diskCharLength/4.0), // major axis length: mean value & standard deviation
NormalDistro(diskCharLength*0.8, diskCharLength*0.8/0.4), // minor axis length: mean value & standard deviation
NormalDistro(toRadians(0.0), toRadians(10.0)), // rotation around x-axis: mean value & standard deviation
NormalDistro(toRadians(0.0), toRadians(35.0)), // rotation around y-axis: mean value & standard deviation
NormalDistro(toRadians(0.0), toRadians(5.0))); // rotation around z-axis: mean value & standard deviation
};
///////////////////////////////////////////////////////////////////////
// 3. Define constraints that should be fulfilled among the entities //
///////////////////////////////////////////////////////////////////////
// Define constraints between entities of orientation 1
using Constraints = EntityNetworkConstraints<ctype>;
Constraints constraints;
constraints.setMinDistance(1.0);
constraints.setMinIntersectingAngle(toRadians(15.0));
constraints.setMinIntersectionMagnitude(0.1);
constraints.setMinIntersectionDistance(0.1);
///////////////////////////
// 4. Network generation //
///////////////////////////
// allow passing number of entities via command line. This is used in the
// test suite in order to speed up testing time. We default to 0, which
// represents the case of the actual example -> construction until extraction
// well is reached
const std::size_t numTargetEntities = argc > 1 ? std::stoi(argv[1]) : 0;
const bool useTargetNumber = numTargetEntities > 0;
const auto setId = Id(1);
// we want the network to not grow too much outside a predefined
// bounding box in order to enforce a predominantly lateral extension
const ctype xMin = injWellBase.center().x() - wellDistance*0.15;
const ctype xMax = extWellBase.center().x() + wellDistance*0.15;
const ctype yMin = injWellBase.center().y() - wellDistance*0.2;
const ctype yMax = extWellBase.center().y() + wellDistance*0.2;
const ctype zMin = injWellBase.center().z() - injectionLength*1.0;
const ctype zMax = extWellBase.center().z() + injectionLength*1.0;
BoundingBox networkBBox(xMin, yMin, zMin, xMax, yMax, zMax);
// returns bbox for sampling for the bbox of the current network
// - increases the network bounding box to make space for the next entity
// - takes the minimum between the result and networkBBox
auto getSampleBBox = [&] (const auto& curNetBbox) -> BoundingBox
{
using std::min; using std::max;
return BoundingBox{max(xMin, curNetBbox.xMin() - 0.5*diskCharLength),
max(yMin, curNetBbox.yMin() - 0.5*diskCharLength),
max(zMin, curNetBbox.zMin() - 0.3*diskCharLength),
min(xMax, curNetBbox.xMax() + 0.5*diskCharLength),
min(yMax, curNetBbox.yMax() + 0.5*diskCharLength),
min(zMax, curNetBbox.zMax() + 0.3*diskCharLength)};
};
// sample an entity within in the injection/extraction wells
// these mark the start and end of the network to be created
const auto injStageCenter = injWellBase.center() + wellNormal*(0.5*injectionLength);
const auto extStageCenter = extWellBase.center() + wellNormal*(0.5*injectionLength);
const auto injEntity = makeDiskSampler(getSampleBBox(getBoundingBox(injStageCenter)))();
const auto extEntity = makeDiskSampler(getSampleBBox(getBoundingBox(extStageCenter)))();
std::vector<Disk> entitySet;
entitySet.reserve(100);
entitySet.push_back(injEntity);
// stopping criterion
// -> for test: numTargetEntities
// -> as example: when extraction well is reached
bool reachedExtWell = false;
auto finished = [&] ()
{
return useTargetNumber ? entitySet.size() >= numTargetEntities
: reachedExtWell;
};
ctype area = 0.0;
SamplingStatus status;
auto sampleBox = getSampleBBox(getBoundingBox(injEntity));
auto diskSampler = makeDiskSampler(sampleBox);
while (!finished())
{
auto candidate = diskSampler();
// we only want entities that interect at least one other
if (std::all_of(entitySet.begin(), entitySet.end(),
[&candidate] (const auto& e)
{ return isEmptyIntersection(intersect(e, candidate)); }))
{ status.increaseRejectedCounter("no intersection"); continue; }
// if too little of the candidate is inside the network bounding box, reject it
const auto containedArea = computeContainedMagnitude(candidate, networkBBox);
if (containedArea < 0.25*candidate.area())
{ status.increaseRejectedCounter("small contained area"); continue; }
// check constraints to other entities
const auto eval = constraints.evaluate(entitySet, candidate);
if (eval.violationDetected())
{ status.increaseRejectedCounter(eval.violationLabel()); continue; }
// increase sampling box for next candidate
sampleBox = getSampleBBox(sampleBox + getBoundingBox(candidate));
diskSampler = makeDiskSampler(sampleBox);
// check is extraction well has been reached
if (!isEmptyIntersection(intersect(candidate, extEntity)))
reachedExtWell = true;
area += containedArea;
entitySet.emplace_back(std::move(candidate));
status.increaseCounter(setId);
status.print();
}
entitySet.push_back(extEntity);
status.printRejectionData();
// print the final entity density
const auto domainVolume = computeMagnitude(domain);
const auto numEntities = entitySet.size();
std::cout << "\nCreated a network consisting of " << numEntities << " entities\n";
std::cout << "Volume of the domain: " << domainVolume << std::endl;
std::cout << "Area of the network: " << area << " m²"
<< ", which corresponds to a density of " << area/domainVolume << " m²/m³" << std::endl;
/////////////////////////////////////////////////////////////////////////////////////////
// 5. The entities of the network have been created. We can now construct the network. //
/////////////////////////////////////////////////////////////////////////////////////////
// construct and write a contained network, i.e. write out both network and domain.
std::cout << "Building and writing network" << std::endl;
ContainedEntityNetworkBuilder<ctype> builder;
// add sub-domains (TODO: Ghost sub-domains etc)
builder.addConfiningSubDomain(domain, Id(1));
builder.addConfiningSubDomain(injWell, Id(2));
builder.addConfiningSubDomain(extWell, Id(3));
builder.addSubDomainEntities(entitySet, Id(1));
// now we can build and write out the network in Gmsh file format
GmshWriter gmshWriter(builder.build());
gmshWriter.write("network",
/*sizeAtEntities*/diskCharLength/4.0,
/*sizeInDomain*/diskCharLength*5.0);
// print time it took to generate the network
const auto end = std::chrono::steady_clock::now();
const auto duration = std::chrono::duration<double>(end-start).count();
std::cout << "Overall CPU time was " << duration << " seconds." << std::endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment