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

Merge branch 'feature/applications' into 'master'

Feature/applications

See merge request DennisGlaeser/frackit!43
parents 36fb9081 1a5e4ad8
......@@ -108,8 +108,8 @@ endif (DOXYGEN_FOUND)
include_directories( ${CMAKE_SOURCE_DIR} )
# include subdirectories
add_subdirectory(appl)
add_subdirectory(cmake)
add_subdirectory(doc)
add_subdirectory(examples)
add_subdirectory(frackit)
add_subdirectory(test)
add_subdirectory(example1)
add_subdirectory(example2)
add_subdirectory(example3)
frackit_add_application(NAME example1
SOURCES example1.cc
LABELS examples)
set(CMAKE_BUILD_TYPE Debug)
#include <random>
#include <iostream>
#include <frackit/common/id.hh>
#include <frackit/io/gmshwriter.hh>
#include <frackit/sampling/makeuniformpointsampler.hh>
#include <frackit/sampling/quadrilateralsampler.hh>
#include <frackit/sampling/status.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/box.hh>
#include <frackit/entitynetwork/constraints.hh>
#include <frackit/entitynetwork/networkbuilder.hh>
//! create a network of 3d quadrilaterals
int main()
{
using namespace Frackit;
// We consider 3d space here
static constexpr int worldDimension = 3;
// Define the type used for coordinates
using ctype = double;
// Internal geometry type for 3d quadrilaterals
using Quad = Quadrilateral<ctype, worldDimension>;
// Define a domain (here: unit cube) in which the quadrilaterals should be created.
// Boxes are created by providing xmin, ymin, zmin and xmax, ymax and zmax in constructor.
Box<ctype> domain(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
// We now create a sampler instance that uniformly samples points within this box.
// These points will be used as the quadrilateral centers.
auto pointSampler = makeUniformPointSampler(domain);
// Sampler class for quadrilaterals. Per default, this uses
// uniform distributions for all parameters defining the quadrilaterals.
// Quadrilateral samplers require distributions for strike angle, dip angle,
// edge length (see class description for more details). Moreover, we define
// a minimum edge length.
using QuadSampler = QuadrilateralSampler<worldDimension>;
using Distro = std::normal_distribution<ctype>;
QuadSampler quadSampler1(pointSampler,
Distro(toRadians(45.0), toRadians(5.0)), // strike angle: mean value & standard deviation
Distro(toRadians(90.0), toRadians(5.0)), // dip angle: mean value & standard deviation
Distro(0.5, 0.1), // edge length: mean value & standard deviation
0.05); // threshold for minimum edge length
// We use this sampler to create quadrilaterals based on the distributions.
// However, we want to create another set of quadrilaterals, whose entities
// are approximately orthogonal to those defined with the above sampler.
QuadSampler quadSampler2(makeUniformPointSampler(domain), // use a new point sampler instance!
Distro(toRadians(0.0), toRadians(5.0)), // strike angle: mean value & standard deviation
Distro(toRadians(0.0), toRadians(5.0)), // dip angle: mean value & standard deviation
Distro(0.5, 0.1), // edge length: mean value & standard deviation
0.05); // threshold for minimum edge length
// We want to enforce some constraints on the set of quadrilaterals.
// In particular, for entities of the same set we want a minimum spacing
// distance of 5cm, and the quadrilaterals must not intersect in angles
// less than 30°. Moreover, if they intersect, we don't want intersection
// edges whose length is smaller than 5cm, and, the intersection should not
// be too close to the boundary of one of two intersecting quadrilaterals. Here: 5cm.
EntityNetworkConstraints constraintsOnSelf;
constraintsOnSelf.setMinDistance(0.05);
constraintsOnSelf.setMinIntersectingAngle(toRadians(30.0));
constraintsOnSelf.setMinIntersectionMagnitude(0.05);
constraintsOnSelf.setMinIntersectionDistance(0.05);
// with respect to entities of the other set, we want to have larger intersection angles
auto constraintsOnOther = constraintsOnSelf;
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
// container to store created entities
std::vector<Quad> entitySet1;
std::vector<Quad> entitySet2;
// we give unique identifiers to both entity sets
const Id idSet1(1);
const Id idSet2(2);
// use the status class to define when to stop sampling
SamplingStatus status;
status.setTargetCount(idSet1, 15); // we want 15 entities in set 1
status.setTargetCount(idSet2, 15); // we want 15 entities in set 2
// start sampling into set 1 and keep alternating
bool sampleIntoSet1 = true;
std::cout << "\n --- Start entity sampling ---\n" << std::endl;
while (!status.finished())
{
// sample a quadrilateral, alternating between sampler 1 and sampler 2
auto quad = sampleIntoSet1 ? quadSampler1() : quadSampler2();
auto& entitySet = sampleIntoSet1 ? entitySet1 : entitySet2;
const auto& otherEntitySet = sampleIntoSet1 ? entitySet2 : entitySet1;
// sample again if constraints w.r.t. other
// entities of this set are not fulfilled
if (!constraintsOnSelf.evaluate(entitySet, quad))
{ status.increaseRejectedCounter(); continue; }
// sample again if constraints w.r.t. other
// entities of the other set are not fulfilled
if (!constraintsOnOther.evaluate(otherEntitySet, quad))
{ status.increaseRejectedCounter(); continue; }
// the quadrilateral is admissible
entitySet.push_back(quad);
// tell the status we have a new entity in this set
const auto& setId = sampleIntoSet1 ? idSet1 : idSet2;
status.increaseCounter(setId);
status.print();
// sample into the other set the next time
sampleIntoSet1 = !sampleIntoSet1;
}
std::cout << "\n --- Finished entity sampling ---\n" << std::endl;
// We can now create an entity network from the two sets
EntityNetworkBuilder builder;
builder.addEntities(entitySet1);
builder.addEntities(entitySet2);
const auto network = builder.build();
// This can be written out in Gmsh (.geo) format to be
// meshed by a two-dimensional surface mesh
std::cout << "\n --- Writing .geo file ---\n" << std::endl;
GmshWriter writer(network);
writer.write("network", // filename of the .geo files (will add extension .geo automatically)
0.1); // element size to be used
std::cout << "\n --- Finished writing .geo file ---\n" << std::endl;
return 0;
}
frackit_add_application(NAME example2
SOURCES example2.cc
LABELS examples)
set(CMAKE_BUILD_TYPE Debug)
#include <random>
#include <iostream>
#include <frackit/common/id.hh>
#include <frackit/io/gmshwriter.hh>
#include <frackit/sampling/makeuniformpointsampler.hh>
#include <frackit/sampling/quadrilateralsampler.hh>
#include <frackit/sampling/status.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/box.hh>
#include <frackit/geometry/cylinder.hh>
#include <frackit/magnitude/containedmagnitude.hh>
#include <frackit/entitynetwork/constraints.hh>
#include <frackit/entitynetwork/networkbuilder.hh>
//! create a network of 3d quadrilaterals contained in a cylinder
int main()
{
using namespace Frackit;
// We use the same samplers for quadrilaterals as in example 1.
// Therefore, we do not repeat the comments here. For more details
// please see the main file of example1. We will only add comments
// here in those places that deviate from the first example.
static constexpr int worldDimension = 3;
using ctype = double;
using Quad = Quadrilateral<ctype, worldDimension>;
// In this example we consider a cylindrical domain
Cylinder<ctype> domain(/*radius*/0.5, /*height*/1.0);
// However, we still sample the quadrilateral sample points
// within the unit cube (a sampler on cylinders is also available)
// Note that we shift the unit cube to have the origin in
// the center of the bottom boundary
Box<ctype> unitCube(-0.5, -0.5, 0.0, 0.5, 0.5, 1.0);
using QuadSampler = QuadrilateralSampler<worldDimension>;
using Distro = std::normal_distribution<ctype>;
QuadSampler quadSampler1(makeUniformPointSampler(unitCube),
Distro(toRadians(45.0), toRadians(5.0)),
Distro(toRadians(90.0), toRadians(5.0)),
Distro(0.5, 0.1),
0.05);
QuadSampler quadSampler2(makeUniformPointSampler(unitCube),
Distro(toRadians(0.0), toRadians(5.0)),
Distro(toRadians(0.0), toRadians(5.0)),
Distro(0.5, 0.1),
0.05);
EntityNetworkConstraints constraintsOnSelf;
constraintsOnSelf.setMinDistance(0.05);
constraintsOnSelf.setMinIntersectingAngle(toRadians(30.0));
constraintsOnSelf.setMinIntersectionMagnitude(0.05);
constraintsOnSelf.setMinIntersectionDistance(0.05);
auto constraintsOnOther = constraintsOnSelf;
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
// We can now also enforce constraints w.r.t. the domain boundary
auto constraintsOnBoundary = constraintsOnSelf;
std::vector<Quad> entitySet1;
std::vector<Quad> entitySet2;
const Id idSet1(1);
const Id idSet2(2);
SamplingStatus status;
status.setTargetCount(idSet1, 15);
status.setTargetCount(idSet2, 15);
bool sampleIntoSet1 = true;
std::cout << "\n --- Start entity sampling ---\n" << std::endl;
while (!status.finished())
{
auto quad = sampleIntoSet1 ? quadSampler1() : quadSampler2();
auto& entitySet = sampleIntoSet1 ? entitySet1 : entitySet2;
const auto& otherEntitySet = sampleIntoSet1 ? entitySet2 : entitySet1;
if (!constraintsOnSelf.evaluate(entitySet, quad))
{ status.increaseRejectedCounter(); continue; }
if (!constraintsOnOther.evaluate(otherEntitySet, quad))
{ status.increaseRejectedCounter(); continue; }
// 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; }
if (!constraintsOnBoundary.evaluate(domain.lateralFace(), quad))
{ status.increaseRejectedCounter(); continue; }
// we also want to neglect quadrilaterals of which only
// small fragments remain after confining them to the
// domain. This computes the area of the quad that is
// inside the domain
const auto containedArea = computeContainedMagnitude(quad, domain);
// reject if this is too small (< 0.01 m^2)
if (containedArea < 0.01)
{ status.increaseRejectedCounter(); continue; }
entitySet.push_back(quad);
status.increaseCounter(sampleIntoSet1 ? idSet1 : idSet2);
status.print();
sampleIntoSet1 = !sampleIntoSet1;
}
std::cout << "\n --- Finished entity sampling ---\n" << std::endl;
// We can now create a contained entity network from the two sets,
// which has information on both the entities and the domain.
ContainedEntityNetworkBuilder builder;
// define the domain (single sub-domain) and give it a unique id
builder.addConfiningSubDomain(domain, Id(1));
// define entities to be embedded in this domain
builder.addSubDomainEntities(entitySet1, Id(1));
builder.addSubDomainEntities(entitySet2, Id(1));
const auto network = builder.build();
// This can be written out in Gmsh (.geo) format to be
// meshed by a three-dimensional mesh that is conforming
// to the two-dimensional quadrilaterals
std::cout << "\n --- Writing .geo file ---\n" << std::endl;
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
std::cout << "\n --- Finished writing .geo file ---\n" << std::endl;
return 0;
}
frackit_add_application(NAME example3
SOURCES example3.cc
COMPILE_DEFINITIONS BREPFILE="${CMAKE_SOURCE_DIR}/appl/example3/layers.brep"
LABELS examples)
set(CMAKE_BUILD_TYPE Debug)
#include <iostream>
#include <fstream>
#include <string>
#include <stdexcept>
#include <random>
......@@ -18,7 +16,6 @@
#include <frackit/sampling/disksampler.hh>
#include <frackit/sampling/quadrilateralsampler.hh>
#include <frackit/sampling/multigeometrysampler.hh>
#include <frackit/sampling/sequentialsamplingstrategy.hh>
#include <frackit/sampling/status.hh>
// constraints to be enforced on the network (distance, angles, etc.)
......@@ -74,22 +71,25 @@ int main(int argc, char** argv)
// and disk/quadrilaterals as entities (using the sampled center points) //
//////////////////////////////////////////////////////////////////////////////
// we use the default sampler types -> normal distributions for all parameters
using Distro = std::normal_distribution<ctype>;
// Bounding box of the domain in which we want to place the entities
const auto domainBBox = OCCUtilities::getBoundingBox(networkDomain);
// sampler for disks (orientation 1)
DiskSampler diskSampler(makeUniformPointSampler(domainBBox), // sampler for disk center points
std::normal_distribution<ctype>(30.0, 6.5), // major axis length: mean value & standard deviation
std::normal_distribution<ctype>(24.0, 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
Distro(30.0, 6.5), // major axis length: mean value & standard deviation
Distro(24.0, 4.5), // minor axis length: mean value & standard deviation
Distro(toRadians(0.0), toRadians(7.5)), // rotation around x-axis: mean value & standard deviation
Distro(toRadians(0.0), toRadians(7.5)), // rotation around y-axis: mean value & standard deviation
Distro(toRadians(0.0), toRadians(7.5))); // rotation around z-axis: mean value & standard deviation
// sampler for quadrilaterals (orientation 2)
QuadrilateralSampler<3> quadSampler(makeUniformPointSampler(domainBBox), // sampler for quadrilateral center points
std::normal_distribution<ctype>(toRadians(0.0), toRadians(5.0)), // strike angle: mean value & standard deviation
std::normal_distribution<ctype>(toRadians(45.0), toRadians(5.0)), // dip angle: mean value & standard deviation
std::normal_distribution<ctype>(35.0, 5.0), // edge length: mean value & standard deviation
Distro(toRadians(45.0), toRadians(5.0)), // strike angle: mean value & standard deviation
Distro(toRadians(90.0), toRadians(5.0)), // dip angle: mean value & standard deviation
Distro(45.0, 5.0), // edge length: mean value & standard deviation
5.0); // threshold for minimum edge length
// Define ids for the two entity sets
......@@ -114,26 +114,22 @@ int main(int argc, char** argv)
// Define constraints between entities of orientation 1
using Constraints = EntityNetworkConstraints<ctype>;
Constraints constraints1;
constraints1.setMinDistance(1.5);
constraints1.setMinDistance(2.5);
constraints1.setMinIntersectingAngle(toRadians(25.0));
constraints1.setMinIntersectionMagnitude(2.5);
constraints1.setMinIntersectionDistance(2.0);
constraints1.setMinIntersectionMagnitude(5.0);
constraints1.setMinIntersectionDistance(2.5);
// Define constraints between entities of orientation 2
Constraints constraints2;
constraints2.setMinDistance(10.0);
constraints2.setMinIntersectingAngle(toRadians(25.0));
constraints2.setMinIntersectionMagnitude(2.5);
constraints2.setMinIntersectionDistance(2.0);
// we want to enforce larger spacing between those entities
auto constraints2 = constraints1;
constraints2.setMinDistance(5.0);
// Define constraints between entities of different sets
Constraints constraintsOnOther;
auto constraintsOnOther = constraints1;
constraintsOnOther.setMinDistance(2.5);
constraintsOnOther.setMinIntersectingAngle(toRadians(40.0));
constraintsOnOther.setMinIntersectionMagnitude(2.5);
constraintsOnOther.setMinIntersectionDistance(2.0);
// We can use a constraint matrix to facilitate constraint evaluation
// We can use the constraints matrix to facilitate constraint evaluation
EntityNetworkConstraintsMatrix<Constraints> constraintsMatrix;
constraintsMatrix.addConstraints(constraints1, // constraint instance
IdPair(diskSetId, diskSetId)); // set between which to use these constraints
......@@ -146,11 +142,8 @@ int main(int argc, char** argv)
IdPair(quadSetId, diskSetId)}); // sets between which to use these constraints
// Moreover, we define constraints w.r.t. the domain boundary
EntityNetworkConstraints constraintsOnDomain;
constraintsOnDomain.setMinDistance(1.5);
constraintsOnDomain.setMinIntersectingAngle(toRadians(5.0));
constraintsOnDomain.setMinIntersectionMagnitude(20.0);
constraintsOnDomain.setMinIntersectionDistance(2.0);
auto constraintsOnDomain = constraints1;
constraintsOnDomain.setMinIntersectingAngle(toRadians(15.0));
......@@ -166,8 +159,8 @@ int main(int argc, char** argv)
// Helper class for terminal output of the creation
// progress and definition of stop criterion etc
SamplingStatus status;
status.setTargetCount(diskSetId, 10); // we want 10 entities of orientation 1
status.setTargetCount(quadSetId, 10); // we want 10 entities of orientation 2
status.setTargetCount(diskSetId, 12); // we want 11 entities of orientation 1
status.setTargetCount(quadSetId, 16); // we want 13 entities of orientation 2
// The actual network generation loop
ctype containedNetworkArea = 0.0;
......@@ -184,7 +177,7 @@ int main(int argc, char** argv)
// Moreover, we want to avoid small fragments (< 250 m²)
const auto containedArea = computeContainedMagnitude(geom, networkDomain);
if (containedArea < 250.0)
if (containedArea < 350.0)
{ status.increaseRejectedCounter(); continue; }
// enforce constraints w.r.t. to the other entities
......
install(FILES
FrackitApplicationMacros.cmake
FrackitMacros.cmake
FrackitTestMacros.cmake
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/cmake/modules)
# Introduce a target that triggers the building of all applications
add_custom_target(build_applications)
# processes labels passed to frackit_add_application
function(frackit_declare_application_label)
include(CMakeParseArguments)
set(OPTIONS)
set(SINGLEARGS)
set(MULTIARGS LABELS)
cmake_parse_arguments(arg "${OPTIONS}" "${SINGLEARGS}" "${MULTIARGS}" ${ARGN})
if( (DEFINED arg_UNPARSED_ARGUMENTS) AND NOT ( arg_UNPARSED_ARGUMENTS STREQUAL "" ) )
message(FATAL_ERROR "Unhandled extra arguments given to frackit_declare_application_label(): "
"<${arg_UNPARSED_ARGUMENTS}>")
endif()
foreach(label IN LISTS arg_LABELS)
# Make sure the label is not empty, and does not contain any funny
# characters, in particular regex characters
if(NOT (label MATCHES "[-_0-9a-zA-Z]+"))
message(FATAL_ERROR "Refusing to add label \"${label}\" since it is "
"empty or contains funny characters (characters other than "
"alphanumeric ones and \"-\" or \"_\"; the intent of this restriction "
"is to make construction of the argument to \"ctest -L\" easier")
endif()
set(target "build_${label}_applications")
if(NOT TARGET "${target}")
add_custom_target("${target}")
endif()
endforeach()
endfunction(frackit_declare_application_label)
# function to add an application
function(frackit_add_application)
include(CMakeParseArguments)
set(SINGLEARGS NAME TARGET TIMEOUT)
set(MULTIARGS SOURCES COMPILE_DEFINITIONS COMPILE_FLAGS CMD_ARGS COMMAND CMAKE_GUARD LABELS)
cmake_parse_arguments(ADDAPPL "${OPTIONS}" "${SINGLEARGS}" "${MULTIARGS}" ${ARGN})
# Check whether the parser produced any errors
if(ADDAPPL_UNPARSED_ARGUMENTS)
message(WARNING "Unrecognized arguments ('${ADDAPPL_UNPARSED_ARGUMENTS}') for frackit_add_application!")
endif()
# Check input for validity and apply defaults
if(NOT ADDAPPL_SOURCES AND NOT ADDAPPL_TARGET)
message(FATAL_ERROR "You need to specify either the SOURCES or the TARGET option for frackit_add_application!")
endif()
if(ADDAPPL_SOURCES AND ADDAPPL_TARGET)
message(FATAL_ERROR "You cannot specify both SOURCES and TARGET for frackit_add_application")
endif()
if(NOT ADDAPPL_NAME)
message(FATAL_ERROR "You have to define a name for the application!")
endif()
if(NOT ADDAPPL_COMMAND)
set(ADDAPPL_COMMAND ${ADDAPPL_NAME})
endif()
if(NOT ADDAPPL_TIMEOUT)
set(ADDAPPL_TIMEOUT 3600)
endif()
# Add the executable if it is not already present
if(ADDAPPL_SOURCES)
add_executable(${ADDAPPL_NAME} ${ADDAPPL_SOURCES})
target_compile_definitions(${ADDAPPL_NAME} PUBLIC ${ADDAPPL_COMPILE_DEFINITIONS})
target_compile_options(${ADDAPPL_NAME} PUBLIC ${ADDAPPL_COMPILE_FLAGS})
# link OCC library to each application
target_link_libraries(${ADDAPPL_NAME} ${ADDAPPL_LINK_LIBRARIES} ${OCC_LIBS})
set(ADDAPPL_TARGET ${ADDAPPL_NAME})
endif()
# make sure each label exists and its name is acceptable
frackit_declare_application_label(LABELS ${ADDAPPL_LABELS})
# Make sure to exclude the target from all
set_property(TARGET ${ADDAPPL_TARGET} PROPERTY EXCLUDE_FROM_ALL 1)
# Have build_applications and build_${label}_applications depend on the given target in
# order to trigger the build correctly
add_dependencies(build_applications ${ADDAPPL_TARGET})
foreach(label IN LISTS ADDAPPL_LABELS)
add_dependencies(build_${label}_applications ${ADDAPPL_TARGET})
endforeach()
endfunction(frackit_add_application)
# includes all relevant macros
include(FrackitTestMacros)
include(FrackitApplicationMacros)
......@@ -73,6 +73,7 @@ function(frackit_add_test)
set(ADDTEST_TARGET ${ADDTEST_NAME})
endif()
# Now add the actual test
add_test(NAME ${ADDTEST_NAME} COMMAND "${ADDTEST_COMMAND}" ${ADDTEST_CMD_ARGS})
......
add_subdirectory(layers)
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)
......@@ -143,7 +143,8 @@ intersect_cylinderSurface_planarGeometry(const CylinderSurface<ctype>& cylSurfac
majAxis.invert();
using std::cos;
const ctype majAxisLength = cylSurface.radius()/cos(dn.Angle(ca));
using std::abs;
const ctype majAxisLength = abs(cylSurface.radius()/cos(dn.Angle(ca)));
const auto& cylAxisLine = cylSurface.centerSegment().supportingLine();
const auto center = std::get<Point>(intersect(faceGeomPlane, cylAxisLine, eps));
infEllipse = Ellipse(center, majAxis, minAxis, majAxisLength, cylSurface.radius());
......@@ -174,7 +175,7 @@ intersect_cylinderSurface_planarGeometry(const CylinderSurface<ctype>& cylSurfac
const auto p1 = OCCUtilities::point(v1);
const auto p2 = OCCUtilities::point(v2);
if (!faceGeomIsParallel)
if (faceGeomIsOrthogonal || !faceGeomIsParallel)
{
// select the arc whose center is on the set of given edges
EllipseArc arc1(infEllipse, p1, p2);
......
Supports Markdown
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