manolas
/
MySources
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Backup commit

This commit is contained in:
iasonmanolas 2022-09-14 13:04:05 +02:00
parent 01dff95f9d
commit 266aca08c2
14 changed files with 6146 additions and 5829 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
CMakeLists.txt
View File

@ -45,7 +45,8 @@ target_include_directories(${PROJECT_NAME} PUBLIC ${DRMSimulationModelDir})
#add_compile_definitions(DLIB_DEFINED)
## polyscope
#if(NOT TARGET polyscope AND ${USE_POLYSCOPE})
if(#[[NOT TARGET polyscope AND]] ${USE_POLYSCOPE})
message("Using polyscope")
FetchContent_Declare(polyscope
GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git
GIT_TAG master
@ -55,7 +56,7 @@ target_include_directories(${PROJECT_NAME} PUBLIC ${DRMSimulationModelDir})
# target_include_directories(${PROJECT_NAME} PUBLIC ${polyscope_SOURCE_DIR}/deps/imgui)
target_sources(${PROJECT_NAME} PUBLIC ${polyscope_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.h ${polyscope_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
#endif()
endif()
#if(NOT TARGET polyscope AND ${USE_POLYSCOPE})
# set(POLYSCOPE_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/polyscope)
# download_project(PROJ POLYSCOPE
@ -149,8 +150,9 @@ add_subdirectory(${ARMADILLO_SOURCE_DIR} ${ARMADILLO_BIN_DIR})
target_include_directories(${PROJECT_NAME} PUBLIC ${ARMADILLO_SOURCE_DIR}/include)
add_compile_definitions(ARMA_DONT_USE_WRAPPER)
target_link_libraries(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/armadillo/build/libarmadillo.a" #[[blas lapack]])
find_package(Armadillo REQUIRED)
target_link_libraries(${PROJECT_NAME} PUBLIC ${ARMADILLO_LIBRARIES})
#find_package(Armadillo REQUIRED)
#target_link_libraries(${PROJECT_NAME} PUBLIC ${ARMADILLO_LIBRARIES})
#if(NOT TARGET ThreedBeamFEA)
#FetchContent_Declare(armadillo
# GIT_REPOSITORY https://gitlab.com/conradsnicta/armadillo-code.git
@ -203,7 +205,7 @@ add_compile_definitions(USE_ENSMALLEN)
# LINK_FLAGS "${CHRONO_LINKER_FLAGS}")
#include_directories(${CHRONO_INCLUDE_DIRS})
#target_link_libraries(${PROJECT_NAME} PUBLIC ${CHRONO_LIBRARIES})
#target_include_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/chrono-7.0.3/src" #[[${CHRONO_INCLUDE_DIRS}]] #[["/home/iason/Coding/Libraries/chrono-7.0.3/src"]] #[["/usr/include/irrlicht"]] )
target_include_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/chrono-7.0.3/src" #[[${CHRONO_INCLUDE_DIRS}]] #[["/home/iason/Coding/Libraries/chrono-7.0.3/src"]] #[["/usr/include/irrlicht"]] )
#target_link_directories(${PROJECT_NAME} PRIVATE "/home/iason/Coding/build/external dependencies/CHRONO-src/build/RelWithDebInfo/lib")
#target_link_libraries(${PROJECT_NAME} PUBLIC "/home/iason/Coding/build/external dependencies/CHRONO-src/build/RelWithDebInfo/lib/libChronoEngine.a")
target_link_libraries(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/chrono-7.0.3/build/Release/lib/libChronoEngine.a")
@ -252,11 +254,12 @@ if(NOT TARGET tbb_static AND NOT TARGET tbb)
# option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
# add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
#link_directories(${TBB_BINARY_DIR})
message("Using tbb")
FetchContent_Declare(tbb
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
)
FetchContent_MakeAvailable(tbb)
# target_link_libraries(${PROJECT_NAME} PRIVATE tbb_static)
target_link_libraries(${PROJECT_NAME} PRIVATE tbb)
endif()
target_link_libraries(${PROJECT_NAME} PRIVATE tbb_static)

68
chronoseulersimulationmodel.cpp
View File

@ -106,6 +106,33 @@ void ChronosEulerSimulationModel::parseForces(
}
}
void ChronosEulerSimulationModel::parseForcedDisplacements(
const std::shared_ptr<const SimulationJob>& pJob,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
edgeMeshVertsToChronoNodes,
chrono::ChSystemSMC& my_system) {
assert(!edgeMeshVertsToChronoNodes.empty());
ChState x;
ChStateDelta v;
double t;
for (const std::pair<VertexIndex, Eigen::Vector3d>& forcedDisplacement :
pJob->nodalForcedDisplacements) {
assert(false);
std::cerr << "Forced displacements dont work" << std::endl;
// std::terminate();
const int& constrainedVi = forcedDisplacement.first;
ChVector<double> displacementVector(
pJob->nodalForcedDisplacements.at(constrainedVi));
const std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>& constrainedChronoNode =
edgeMeshVertsToChronoNodes[constrainedVi];
constrainedChronoNode->NodeIntStateGather(0, x, 0, v, t);
std::cout << "state rows:" << x.rows() << std::endl;
std::cout << "state cols:" << x.cols() << std::endl;
// x = x + displacementVector;
constrainedChronoNode->NodeIntStateScatter(0, x, 0, v, t);
}
}
void ChronosEulerSimulationModel::parseConstrainedVertices(
const std::shared_ptr<const SimulationJob>& pJob,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
@ -120,30 +147,17 @@ void ChronosEulerSimulationModel::parseConstrainedVertices(
auto truss = chrono_types::make_shared<ChBody>();
truss->SetBodyFixed(true);
my_system.Add(truss);
const auto& constrainedChronoNode =
const std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>& constrainedChronoNode =
edgeMeshVertsToChronoNodes[constrainedVi];
// Create a constraint at the end of the beam
auto constr_a = chrono_types::make_shared<ChLinkMateGeneric>();
constr_a->SetConstrainedCoords(
auto constraint = chrono_types::make_shared<ChLinkMateGeneric>();
constraint->SetConstrainedCoords(
constrainedDoF.contains(0), constrainedDoF.contains(1),
constrainedDoF.contains(2), constrainedDoF.contains(3),
constrainedDoF.contains(4), constrainedDoF.contains(5));
constr_a->Initialize(constrainedChronoNode, truss, false,
constraint->Initialize(constrainedChronoNode, truss, false,
constrainedChronoNode->Frame(),
constrainedChronoNode->Frame());
// const auto frameNode = constrainedChronoNode->Frame();
my_system.Add(constr_a);
// edgeMeshVertsToChronoNodes[constrainedVi]->SetFixed(true);
// if (vertexIsFullyConstrained) {
// } else {
// std::cerr << "Currently only rigid vertices are handled." <<
// std::endl;
// // SimulationResults simulationResults;
// // simulationResults.converged = false;
// // assert(false);
// // return simulationResults;
// }
my_system.Add(constraint);
}
}
@ -177,6 +191,7 @@ SimulationResults ChronosEulerSimulationModel::executeSimulation(
// -1));
// my_system.SetTimestepperType(chrono::ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
// parse forces
parseForcedDisplacements(pJob, edgeMeshVertsToChronoNodes, my_system);
parseForces(mesh_chronos, edgeMeshVertsToChronoNodes,
pJob->nodalExternalForces);
// parse constrained vertices
@ -198,20 +213,29 @@ SimulationResults ChronosEulerSimulationModel::executeSimulation(
auto solver = chrono_types::make_shared<ChSolverMINRES>();
my_system.SetSolver(solver);
// solver->SetMaxIterations(1e5);
// solver->SetTolerance(1e-12);
solver->SetMaxIterations(1e5);
solver->EnableWarmStart(
true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED
solver->EnableDiagonalPreconditioner(true);
my_system.SetSolverForceTolerance(1e-9);
solver->SetVerbose(false);
// solver->SetTolerance(1e-15);
my_system.SetSolverForceTolerance(1e-15);
SimulationResults simulationResults;
//#ifdef POLYSCOPE_DEFINED
// solver->SetVerbose(true);
// // edgeMeshVertsToChronoNodes[10]->Frame().Move({0, 0, 1e-1});
// simulationResults.converged = my_system.DoEntireKinematics(5e2);
//// edgeMeshVertsToChronoNodes[10]->SetForce({0, 0, 1e6});
//#else
solver->SetVerbose(false);
if (settings.analysisType == Settings::AnalysisType::Linear) {
simulationResults.converged = my_system.DoStaticLinear();
// simulationResults.converged = my_system.DoStaticRelaxing();
} else {
simulationResults.converged = my_system.DoStaticNonlinear();
// simulationResults.converged = my_system.DoStaticRelaxing();
}
//#endif
if (!simulationResults.converged) {
std::cerr << "Simulation failed" << std::endl;
assert(false);

41
chronoseulersimulationmodel.hpp
View File

@ -11,32 +11,41 @@ class ChNodeFEAxyzrot;
} // namespace fea
} // namespace chrono
class ChronosEulerSimulationModel : public SimulationModel
{
class ChronosEulerSimulationModel : public SimulationModel {
std::shared_ptr<chrono::fea::ChMesh> mesh_chronos;
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> edgeMeshVertsToChronoNodes;
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>
edgeMeshVertsToChronoNodes;
static void parseConstrainedVertices(
const std::shared_ptr<const SimulationJob> &pJob,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
chrono::ChSystemSMC &my_system);
const std::shared_ptr<const SimulationJob>& pJob,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
edgeMeshVertsToChronoNodes,
chrono::ChSystemSMC& my_system);
static void parseForces(
const std::shared_ptr<chrono::fea::ChMesh> &mesh_chronos,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
const std::unordered_map<VertexIndex, Vector6d> &nodalExternalForces);
const std::shared_ptr<chrono::fea::ChMesh>& mesh_chronos,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
edgeMeshVertsToChronoNodes,
const std::unordered_map<VertexIndex, Vector6d>& nodalExternalForces);
public:
static void parseForcedDisplacements(
const std::shared_ptr<const SimulationJob>& pJob,
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
edgeMeshVertsToChronoNodes,
chrono::ChSystemSMC& my_system);
public:
ChronosEulerSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob) override;
void setStructure(const std::shared_ptr<SimulationEdgeMesh> &pMesh) override;
SimulationResults executeSimulation(
const std::shared_ptr<SimulationJob>& pJob) override;
void setStructure(const std::shared_ptr<SimulationEdgeMesh>& pMesh) override;
static std::shared_ptr<chrono::fea::ChMesh> convertToChronosMesh_Euler(
const std::shared_ptr<SimulationEdgeMesh> &pMesh,
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronosNodes);
const std::shared_ptr<SimulationEdgeMesh>& pMesh,
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>>&
edgeMeshVertsToChronosNodes);
inline const static std::string label{"Chronos_Euler"};
struct Settings
{
struct Settings {
enum AnalysisType { Linear = 0, NonLinear };
AnalysisType analysisType{NonLinear};
};

297
polymesh.hpp
View File

@ -1,297 +0,0 @@
#ifndef POLYMESH_HPP
#define POLYMESH_HPP
#include "mesh.hpp"
#include "utilities.hpp"
#include "vcg/complex/complex.h"
#include <filesystem>
#include <wrap/io_trimesh/export_obj.h>
#include <wrap/io_trimesh/export_ply.h>
#include <wrap/io_trimesh/import.h>
#ifdef POLYSCOPE_DEFINED
#include <polyscope/surface_mesh.h>
#endif
class PFace;
class PVertex;
class PEdge;
struct PUsedTypes : public vcg::UsedTypes<vcg::Use<PVertex>::AsVertexType,
vcg::Use<PFace>::AsFaceType,
vcg::Use<PEdge>::AsEdgeType>
{};
class PVertex : public vcg::Vertex<PUsedTypes,
vcg::vertex::Coord3d,
vcg::vertex::Normal3d,
vcg::vertex::Mark,
vcg::vertex::Qualityf,
vcg::vertex::BitFlags,
vcg::vertex::VFAdj,
vcg::vertex::VEAdj>
{};
class PEdge : public vcg::Edge<PUsedTypes,
vcg::edge::VertexRef,
vcg::edge::BitFlags,
vcg::edge::EEAdj,
vcg::edge::EFAdj,
vcg::edge::VEAdj,
vcg::edge::EVAdj>
{};
class PFace : public vcg::Face<PUsedTypes,
vcg::face::PolyInfo // this is necessary if you use component in
// vcg/simplex/face/component_polygon.h
// It says "this class is a polygon and the memory for its components
// (e.g. pointer to its vertices will be allocated dynamically")
,
// vcg::face::FHAdj,
vcg::face::PVFAdj,
vcg::face::PFEAdj,
vcg::face::PFVAdj,
vcg::face::PVFAdj,
// vcg::face::PVFAdj,
vcg::face::PFFAdj // Pointer to edge-adjacent face (just like FFAdj )
,
vcg::face::BitFlags // bit flags
,
vcg::face::Qualityf // quality
,
vcg::face::Normal3f // normal
,
vcg::face::Color4b>
{};
class VCGPolyMesh
: public vcg::tri::TriMesh<std::vector<PVertex>, std::vector<PFace>, std::vector<PEdge>>,
public Mesh
{
public:
virtual bool load(const std::filesystem::__cxx11::path &meshFilePath) override
{
int mask;
vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(meshFilePath.c_str(), mask);
int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, meshFilePath.c_str(), mask);
if (error != 0) {
std::cerr << "Could not load polygonal mesh:" << meshFilePath << std::endl;
return false;
}
vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
vcg::tri::Clean<VCGPolyMesh>::RemoveUnreferencedVertex(*this);
//finally remove valence 1 vertices on the border
// vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
return true;
}
// // vcg::tri::io::ImporterOBJ<VCGPolyMesh>::Open();
// // unsigned int mask = 0;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOORD;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTNORMAL;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_VERTCOLOR;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_EDGEINDEX;
// // mask |= nanoply::NanoPlyWrapper<VCGPolyMesh>::IO_FACEINDEX;
// // if (nanoply::NanoPlyWrapper<VCGPolyMesh>::LoadModel(
// // std::filesystem::absolute(filename).c_str(), *this, mask) !=
// // 0) {
// // std::cout << "Could not load tri mesh" << std::endl;
// // }
// }
Eigen::MatrixX3d getVertices() const
{
Eigen::MatrixX3d vertices(VN(), 3);
for (size_t vi = 0; vi < VN(); vi++) {
VCGPolyMesh::CoordType vertexCoordinates = vert[vi].cP();
vertices.row(vi) = vertexCoordinates.ToEigenVector<Eigen::Vector3d>();
}
return vertices;
}
std::vector<std::vector<int>> getFaces() const
{
std::vector<std::vector<int>> faces(FN());
for (const VCGPolyMesh::FaceType &f : this->face) {
const int fi = vcg::tri::Index<VCGPolyMesh>(*this, f);
for (size_t vi = 0; vi < f.VN(); vi++) {
const size_t viGlobal = vcg::tri::Index<VCGPolyMesh>(*this, f.cV(vi));
faces[fi].push_back(viGlobal);
}
}
return faces;
}
// bool load(const std::filesystem::__cxx11::path &meshFilePath)
// {
// const std::string extension = ".ply";
// std::filesystem::path filePath = meshFilePath;
// assert(std::filesystem::path(filePath).extension().string() == extension);
// unsigned int mask = 0;
// mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
// mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
// mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
// mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
// if (vcg::tri::io::Importer<VCGPolyMesh>::Open(*this, filePath.c_str()) != 0) {
// return false;
// }
// label = meshFilePath.filename();
// return true;
// }
bool save(const std::filesystem::__cxx11::path &meshFilePath = std::filesystem::path()) override
{
if (meshFilePath.extension() == ".obj") {
return saveOBJ(meshFilePath);
} else if (meshFilePath.extension() == ".ply") {
return savePLY(meshFilePath);
}
return false;
}
bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path())
{
const std::string extension = ".obj";
std::filesystem::path filePath = objFilePath;
if (filePath.empty()) {
filePath = std::filesystem::current_path().append(getLabel() + extension).string();
} else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) {
filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string();
}
assert(std::filesystem::path(filePath).extension().string() == extension);
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask) != 0) {
return false;
}
return true;
}
bool savePLY(const std::filesystem::path &objFilePath = std::filesystem::path())
{
const std::string extension = ".ply";
std::filesystem::path filePath = objFilePath;
if (filePath.empty()) {
filePath = std::filesystem::current_path().append(getLabel() + extension).string();
} else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) {
filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string();
}
assert(std::filesystem::path(filePath).extension().string() == extension);
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
mask |= vcg::tri::io::Mask::IOM_FACEINDEX;
mask |= vcg::tri::io::Mask::IOM_FACECOLOR;
if (vcg::tri::io::ExporterPLY<VCGPolyMesh>::Save(*this, filePath.string().c_str(), mask, false) != 0) {
return false;
}
return true;
}
#ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true)
{
auto vertices = getVertices();
auto faces = getFaces();
PolyscopeInterface::init();
polyscope::SurfaceMesh *polyscopeHandle_mesh = polyscope::registerSurfaceMesh(label,
vertices,
faces);
const double drawingRadius = 0.002;
polyscopeHandle_mesh->setEnabled(shouldEnable);
polyscopeHandle_mesh->setEdgeWidth(drawingRadius);
if (desiredColor.has_value()) {
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_mesh->setSurfaceColor(desiredColor_glm);
}
return polyscopeHandle_mesh;
}
#endif
void moveToCenter()
{
CoordType centerOfMass(0, 0, 0);
for (int vi = 0; vi < VN(); vi++) {
centerOfMass += vert[vi].cP();
}
centerOfMass /= VN();
vcg::tri::UpdatePosition<VCGPolyMesh>::Translate(*this, -centerOfMass);
}
/*
* Returns the average distance from the center of each edge to the center of its face over the whole mesh
* */
double getAverageFaceRadius() const
{
double averageFaceRadius = 0;
for (int fi = 0; fi < FN(); fi++) {
const VCGPolyMesh::FaceType &f = face[fi];
CoordType centerOfFace(0, 0, 0);
for (int vi = 0; vi < f.VN(); vi++) {
centerOfFace = centerOfFace + f.cP(vi);
}
centerOfFace /= f.VN();
double faceRadius = 0;
// for (int face_ei = 0; face_ei < f.EN(); face_ei++) {
// std::cout << "fi:" << getIndex(f) << std::endl;
// auto vps = f.FVp(0);
// auto vpe = vps;
for (int i = 0; i < f.VN(); i++) {
faceRadius += vcg::Distance(centerOfFace, (f.cP0(i) + f.cP1(i)) / 2);
}
// }
const int faceEdges = f.VN(); //NOTE: When does this not hold?
faceRadius /= faceEdges;
averageFaceRadius += faceRadius;
}
averageFaceRadius /= FN();
return averageFaceRadius;
}
bool copy(VCGPolyMesh &copyFrom)
{
vcg::tri::Append<VCGPolyMesh, VCGPolyMesh>::MeshCopy(*this, copyFrom);
label = copyFrom.getLabel();
// eigenEdges = mesh.getEigenEdges();
// if (eigenEdges.rows() == 0) {
// getEdges(eigenEdges);
// }
// eigenVertices = mesh.getEigenVertices();
// if (eigenVertices.rows() == 0) {
// getVertices();
// }
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::AllocateEdge(*this);
vcg::tri::UpdateTopology<VCGPolyMesh>::EdgeEdge(*this);
// vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
return true;
}
VCGPolyMesh(VCGPolyMesh &copyFrom) { copy(copyFrom); }
VCGPolyMesh() {}
template<typename MeshElement>
size_t getIndex(const MeshElement &meshElement) const
{
return vcg::tri::Index<VCGPolyMesh>(*this, meshElement);
}
};
using ConstVCGPolyMesh = VCGPolyMesh;
#endif // POLYMESH_HPP

642
reducedmodelevaluator.cpp
View File

@ -1,121 +1,45 @@
#include "reducedmodelevaluator.hpp"
#include <execution>
#include <filesystem>
#include "hexagonremesher.hpp"
#include "reducedmodel.hpp"
#include "reducedmodeloptimizer.hpp"
#include "simulationmodelfactory.hpp"
#include "trianglepatterngeometry.hpp"
#include <execution>
#include <filesystem>
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
using PatternVertexIndex = VertexIndex;
using ReducedModelVertexIndex = VertexIndex;
ReducedModelEvaluator::ReducedModelEvaluator()
{
ReducedModelEvaluator::ReducedModelEvaluator() {
pTileIntoSurface = [&]() {
std::istringstream inputStream_tileIntoTriSurface(tileIntoSurfaceFileContent);
std::istringstream inputStream_tileIntoTriSurface(
tileIntoSurfaceFileContent);
VCGTriMesh tileInto_triMesh;
const bool surfaceLoadSuccessfull = tileInto_triMesh.load(inputStream_tileIntoTriSurface);
const bool surfaceLoadSuccessfull =
tileInto_triMesh.load(inputStream_tileIntoTriSurface);
tileInto_triMesh.setLabel("instantMeshes_plane_34");
assert(surfaceLoadSuccessfull);
return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
}();
}
//double ReducedModelEvaluator::evaluateOptimizationSettings(
// const ReducedModelOptimization::Settings &optimizationSettings,
// const std::vector<std::shared_ptr<PatternGeometry>> &pPatterns,
// std::vector<ReducedModelEvaluator::Results> &patternEvaluationResults)
//{
// assert(!pPatterns.empty());
// double optimizationError = 0;
// auto start = std::chrono::high_resolution_clock::now();
// std::vector<double> averageNormalizedError(pPatterns.size(), 0);
// patternEvaluationResults.clear();
// patternEvaluationResults.resize(pPatterns.size());
// std::for_each(
// // std::execution::par_unseq,
// pPatterns.begin(),
// pPatterns.end(),
// [&](const std::shared_ptr<ConstPatternGeometry> &pPattern) {
// // std::cout << "Optimizing " << pPattern->getLabel() << std::endl;
// ReducedModelOptimization::Results optimizationResults;
// ReducedModelOptimizer optimizer;
// optimizer.optimize(*pPattern, optimizationSettings, optimizationResults);
// const auto evaluationResults
// = ReducedModelEvaluator::evaluateReducedModel(optimizationResults,
// tileIntoSurfaceFilePath,
// scenariosDirPath,
// fullPatternTessellatedResultsDirPath);
// const double averageNormalizedErrorOfPattern
// = std::reduce(evaluationResults.distances_normalizedDrm2reduced.begin(),
// evaluationResults.distances_normalizedDrm2reduced.end())
// / evaluationResults.distances_normalizedDrm2reduced.size();
// const int patternIndex = &pPattern - &patterns[0];
// averageNormalizedError[patternIndex] = averageNormalizedErrorOfPattern;
// patternsEvaluationResults[patternIndex] = evaluationResults;
// });
// const double strategyAverageNormalizedError = std::reduce(std::execution::par_unseq,
// averageNormalizedError.begin(),
// averageNormalizedError.end())
// / pPointers.size();
// const auto totalDuration_min = std::chrono::duration_cast<std::chrono::seconds>(
// std::chrono::high_resolution_clock::now() - start)
// .count()
// / 60.0;
// std::cout << "Optimized pattern(s) in:" << totalDuration_min << " minutes." << std::endl;
// std::cout << "Average time per pattern:" << totalDuration_min / patternsPointers.size()
// << " minutes." << std::endl;
// std::cout << "Objective value:" << strategyAverageNormalizedError << std::endl;
// return strategyAverageNormalizedError;
// // std::cout << "After:" << ++numberOfOptimizationRoundsExecuted << " iterations." << std::endl;
//}
ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResult)
{
const std::filesystem::path scenariosDirectoryPath
= "/home/iason/Coding/Projects/Approximating shapes with flat "
ReducedModelOptimization::Results& optimizationResult) {
const std::filesystem::path scenariosDirectoryPath =
"/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/ReducedModelEvaluator/Scenarios";
const std::filesystem::path fullPatternTessellatedResultsDirectoryPath
= "/home/iason/Coding/Projects/Approximating shapes with flat "
const std::filesystem::path fullPatternTessellatedResultsDirectoryPath =
"/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/ReducedModelEvaluator/TessellatedResults";
return evaluateReducedModel(optimizationResult,
scenariosDirectoryPath,
return evaluateReducedModel(optimizationResult, scenariosDirectoryPath,
fullPatternTessellatedResultsDirectoryPath);
}
//void ReducedModelEvaluator::printResults(const Results &evaluationResults,
// const std::string &resultsLabel,
// const std::filesystem::path &resultsOutputPath)
//{
// const bool outputPathIsDirectory = !resultsOutputPath.empty()
// && !resultsOutputPath.has_extension();
// const bool outputPathIsFile = !resultsOutputPath.empty() && resultsOutputPath.has_extension();
// assert(outputPathIsDirectory && outputPathIsFile);
// if (outputPathIsDirectory) {
// std::filesystem::create_directories(resultsOutputPath);
// }
//#else
// std::filesystem::path csvOutputFilePath;
// bool shouldOverwrite = false;
// if (outputPathIsDirectory) {
// csvOutputFilePath = std::filesystem::path(resultsOutputPath)
// .append("distances_" + resultsLabel + ".csv")
// .string();
// shouldOverwrite = true;
// } else if (outputPathIsFile) {
// csvOutputFilePath = resultsOutputPath;
// }
// csvFile csvOutput(csvOutputFilePath, shouldOverwrite);
// printResults(evaluationResults, resultsLabel, csvOutput);
//}
void ReducedModelEvaluator::printResults(const Results &evaluationResults,
const std::string &resultsLabel)
{
void ReducedModelEvaluator::printResults(const Results& evaluationResults,
const std::string& resultsLabel) {
csvFile csvOutputToCout({}, true);
Settings exportSettings;
exportSettings.exportingDirection = Vertical;
@ -124,15 +48,12 @@ void ReducedModelEvaluator::printResults(const Results &evaluationResults,
printResults(evaluationResults, exportSettings, csvOutputToCout);
}
void ReducedModelEvaluator::printResults(const Results &evaluationResults,
const Settings &settings,
csvFile &csvOutput)
{
void ReducedModelEvaluator::printResults(const Results& evaluationResults,
const Settings& settings,
csvFile& csvOutput) {
if (settings.shouldWriteHeader) {
csvOutput << csvExportingDataStrings[settings.exportingData];
printHeader(settings, csvOutput);
// csvOutput << "Average error";
// csvOutput<<"Cumulative error";
csvOutput << endrow;
}
if (!settings.resultsLabel.empty()) {
@ -145,13 +66,15 @@ void ReducedModelEvaluator::printResults(const Results &evaluationResults,
}
}
void ReducedModelEvaluator::printHeader(const Settings &settings, csvFile &csvOutput)
{
void ReducedModelEvaluator::printHeader(const Settings& settings,
csvFile& csvOutput) {
if (settings.exportingDirection == Horizontal) {
// csvOutput << "Job label";
for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
for (int jobIndex = 0;
jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
jobIndex++) {
const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
const std::string& jobLabel =
ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
csvOutput << jobLabel;
}
} else {
@ -161,20 +84,20 @@ void ReducedModelEvaluator::printHeader(const Settings &settings, csvFile &csvOu
}
}
void ReducedModelEvaluator::printResultsHorizontally(const Results &evaluationResults,
csvFile &csvOutput)
{
//print header
//print raw error
void ReducedModelEvaluator::printResultsHorizontally(
const Results& evaluationResults,
csvFile& csvOutput) {
// print header
// print raw error
constexpr bool shouldPrintRawError = false;
if (shouldPrintRawError) {
// csvOutput << "drm2Reduced";
double sumOfFull2Reduced = 0;
int numOfNonEvaluatedScenarios = 0;
for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
for (int jobIndex = 0;
jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
jobIndex++) {
const double &distance_fullDrmToReduced = evaluationResults
.distances_drm2reduced[jobIndex];
const double& distance_fullDrmToReduced =
evaluationResults.distances_drm2reduced[jobIndex];
if (distance_fullDrmToReduced == -1) {
csvOutput << "notEvaluated";
numOfNonEvaluatedScenarios++;
@ -183,19 +106,14 @@ void ReducedModelEvaluator::printResultsHorizontally(const Results &evaluationRe
sumOfFull2Reduced += distance_fullDrmToReduced;
}
}
// const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
// - numOfNonEvaluatedScenarios;
// const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
// csvOutput << averageDistance_full2Reduced;
// csvOutput << endrow;
}
//print normalized error
// csvOutput << "norm_drm2Reduced";
// print normalized error
double sumOfNormalizedFull2Reduced = 0;
for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
for (int jobIndex = 0;
jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
jobIndex++) {
const double &distance_normalizedFullDrmToReduced
= evaluationResults.distances_normalizedDrm2reduced[jobIndex];
const double& distance_normalizedFullDrmToReduced =
evaluationResults.distances_normalizedDrm2reduced[jobIndex];
if (distance_normalizedFullDrmToReduced == -1) {
csvOutput << "notEvaluated";
} else {
@ -203,15 +121,11 @@ void ReducedModelEvaluator::printResultsHorizontally(const Results &evaluationRe
sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
}
}
// const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
// / numOfEvaluatedScenarios;
// csvOutput << averageDistance_normalizedFull2Reduced;
// csvOutput << endrow;
}
void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResults,
csvFile &csvOutput)
{
void ReducedModelEvaluator::printResultsVertically(
const Results& evaluationResults,
csvFile& csvOutput) {
#ifdef POLYSCOPE_DEFINED
csvOutput << "pattern2Reduced"
<< "norm_pattern2Reduced";
@ -225,13 +139,16 @@ void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResu
double sumOfFull2Reduced = 0;
double sumOfNormalizedFull2Reduced = 0;
int numOfNonEvaluatedScenarios = 0;
for (int jobIndex = 0; jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
for (int jobIndex = 0;
jobIndex < ReducedModelEvaluator::scenariosTestSetLabels.size();
jobIndex++) {
const double &distance_fullDrmToReduced = evaluationResults.distances_drm2reduced[jobIndex];
const double &distance_normalizedFullDrmToReduced
= evaluationResults.distances_normalizedDrm2reduced[jobIndex];
const double& distance_fullDrmToReduced =
evaluationResults.distances_drm2reduced[jobIndex];
const double& distance_normalizedFullDrmToReduced =
evaluationResults.distances_normalizedDrm2reduced[jobIndex];
#ifndef POLYSCOPE_DEFINED
const std::string &jobLabel = ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
const std::string& jobLabel =
ReducedModelEvaluator::scenariosTestSetLabels[jobIndex];
csvOutput << jobLabel;
#endif
if (distance_fullDrmToReduced == -1) {
@ -239,22 +156,25 @@ void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResu
<< "notEvaluated";
numOfNonEvaluatedScenarios++;
} else {
csvOutput << distance_fullDrmToReduced << distance_normalizedFullDrmToReduced;
csvOutput << distance_fullDrmToReduced
<< distance_normalizedFullDrmToReduced;
sumOfFull2Reduced += distance_fullDrmToReduced;
sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
}
csvOutput << endrow;
// sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
}
const int numOfEvaluatedScenarios = ReducedModelEvaluator::scenariosTestSetLabels.size()
- numOfNonEvaluatedScenarios;
const double averageDistance_full2Reduced = sumOfFull2Reduced / numOfEvaluatedScenarios;
const double averageDistance_normalizedFull2Reduced = sumOfNormalizedFull2Reduced
/ numOfEvaluatedScenarios;
const int numOfEvaluatedScenarios =
ReducedModelEvaluator::scenariosTestSetLabels.size() -
numOfNonEvaluatedScenarios;
const double averageDistance_full2Reduced =
sumOfFull2Reduced / numOfEvaluatedScenarios;
const double averageDistance_normalizedFull2Reduced =
sumOfNormalizedFull2Reduced / numOfEvaluatedScenarios;
#ifndef POLYSCOPE_DEFINED
csvOutput << "Average error";
#endif
csvOutput << averageDistance_full2Reduced << averageDistance_normalizedFull2Reduced;
csvOutput << averageDistance_full2Reduced
<< averageDistance_normalizedFull2Reduced;
csvOutput << endrow;
#ifndef POLYSCOPE_DEFINED
csvOutput << "Cumulative error";
@ -264,106 +184,67 @@ void ReducedModelEvaluator::printResultsVertically(const Results &evaluationResu
csvOutput << endrow;
}
//void ReducedModelEvaluator::createFullAndReducedPatternTessellations(){
//}
//ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
// std::vector<ReducedModelOptimization::Results> &optimizationResults,
// const std::filesystem::path &tileInto_triMesh_filePath,
// const std::filesystem::path &scenariosDirectoryPath,
// // const std::filesystem::path &reducedPatternFilePath,
// const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath)
//{
// //Load surface
// std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
// VCGTriMesh tileInto_triMesh;
// const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filePath);
// assert(surfaceLoadSuccessfull);
// return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
// }();
// const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0),
// (optimizationResult.baseTriangle.cP(1)
// + optimizationResult.baseTriangle.cP(
// 2))
// / 2);
// pTileIntoSurface->moveToCenter();
// const double scaleFactor = optimizedBaseTriangleHeight
// / pTileIntoSurface->getAverageFaceRadius();
// vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
// //Tile full pattern into surface
//}
ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResult,
const std::filesystem::path &scenariosDirectoryPath,
const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath)
{
// const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResult.baseTriangle.cP(0),
// (optimizationResult.baseTriangle.cP(1)
// + optimizationResult.baseTriangle.cP(
// 2))
// / 2);
pTileIntoSurface->moveToCenter();
const double scaleFactor = optimizationResult.settings.targetBaseTriangleSize
/ pTileIntoSurface->getAverageFaceRadius();
vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
// tileIntoSurface.registerForDrawing();
// polyscope::show();
ReducedModelOptimization::Results& optimizationResult,
const std::filesystem::path& scenariosDirectoryPath,
const std::filesystem::path& patternTessellatedResultsDirectoryPath) {
// Apply optimization results to the reduced model
ReducedModel reducedModel;
reducedModel.deleteDanglingVertices();
std::unordered_map<std::string, double> optimalXVariables_set(
optimizationResult.optimalXNameValuePairs.begin(),
optimizationResult.optimalXNameValuePairs.end());
reducedModel.updateBaseTriangleGeometry_R(optimalXVariables_set.at("R"));
reducedModel.updateBaseTriangleGeometry_theta(
optimalXVariables_set.at("Theta"));
//Tile full pattern into surface
std::vector<PatternGeometry> fullPatterns(1);
fullPatterns[0].copy(optimizationResult.baseTriangleFullPattern);
//// Base triangle pattern might contain dangling vertices.Remove those
fullPatterns[0].interfaceNodeIndex = 3;
fullPatterns[0].deleteDanglingVertices();
// reducedModel.registerForDrawing();
// Scale tile-into surface
pTileIntoSurface->moveToCenter();
const double scaleFactor =
optimizationResult.settings.targetBaseTriangleSize /
pTileIntoSurface->getAverageFaceRadius();
vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
#ifdef POLYSCOPE_DEFINED
pTileIntoSurface->registerForDrawing(color_tileIntoSurface);
#endif
// Tile full pattern into surface
std::vector<PatternGeometry> patterns(1);
patterns[0].copy(optimizationResult.baseTrianglePattern);
patterns[0].interfaceNodeIndex = 3;
patterns[0].deleteDanglingVertices();
std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledFullVi;
std::shared_ptr<PatternGeometry> pTiledFullPattern
= PatternGeometry::tilePattern(fullPatterns,
{},
*pTileIntoSurface,
std::shared_ptr<PatternGeometry> pTilled_pattern =
PatternGeometry::tilePattern(patterns, {}, *pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledFullVi,
perPatternIndexTiledFullPatternEdgeIndices);
pTiledFullPattern->setLabel("Tiled_full_patterns");
// pTiledFullPattern->registerForDrawing();
//Tile reduced pattern into surface
ReducedModel reducedModel;
reducedModel.deleteDanglingVertices();
std::unordered_map<std::string, double>
optimalXVariables_set(optimizationResult.optimalXNameValuePairs.begin(),
optimizationResult.optimalXNameValuePairs.end());
reducedModel.updateBaseTriangleGeometry_R(optimalXVariables_set.at("R"));
reducedModel.updateBaseTriangleGeometry_theta(optimalXVariables_set.at("Theta"));
// const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle();
// ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned(
// optimizationResult, reducedPatternBaseTriangle, reducedModel);
pTilled_pattern->setLabel("Tilled_pattern");
std::vector<PatternGeometry> reducedPatterns(1);
reducedPatterns[0].copy(reducedModel);
std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
// Tile reduced pattern into surface
std::vector<PatternGeometry> reducedModels(1);
reducedModels[0].copy(reducedModel);
std::vector<std::vector<size_t>>
perPatternIndexTiledReducedPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledReducedVi;
std::shared_ptr<PatternGeometry> pTiledReducedPattern
= PatternGeometry::tilePattern(reducedPatterns,
{0},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
std::shared_ptr<PatternGeometry> pTilled_reducedModel =
PatternGeometry::tilePattern(
reducedModels, {0}, *pTileIntoSurface, perSurfaceFacePatternIndices,
tileIntoEdgeToTiledReducedVi,
perPatternIndexTiledReducedPatternEdgeIndices);
pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
#ifdef POLYSCOPE_DEFINED
pTiledReducedPattern->registerForDrawing();
#endif
pTilled_reducedModel->setLabel("Tilled_reduced_model");
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
fullToReducedViMap; //of only the common vertices
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
reducedToFullViMap; //of only the common vertices
std::unordered_map<PatternVertexIndex, ReducedModelVertexIndex>
fullToReducedViMap; // of only the common vertices
std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>
reducedToFullViMap; // of only the common vertices
for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] =
tileIntoEdgeToTiledReducedVi[ei];
}
constructInverseMap(fullToReducedViMap, reducedToFullViMap);
@ -375,179 +256,208 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
}
//Create simulation meshes
// Create simulation meshes
////Tessellated full pattern simulation mesh
std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
//NOTE: Those should be derived from the optimization results instead of hardcoding them
// const double beamWidth = optimizationResult.settings.patternBeamDimensions.getWidth();
// const double beamHeight = optimizationResult.settings.patternBeamDimensions.getHeight();
pTiledFullPattern_simulationMesh->setBeamCrossSection(
optimizationResult.settings.patternBeamDimensions);
if (optimizationResult.fullPatternYoungsModulus == 0) {
std::shared_ptr<SimulationEdgeMesh> pSimulationEdgeMesh_tilledPattern =
std::make_shared<SimulationEdgeMesh>(*pTilled_pattern);
pSimulationEdgeMesh_tilledPattern->setBeamCrossSection(
optimizationResult.settings.beamDimensions_pattern);
if (optimizationResult.settings.youngsModulus_pattern == 0) {
std::cerr << "Full pattern's young modulus not found." << std::endl;
std::terminate();
}
pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
optimizationResult.fullPatternYoungsModulus);
pTiledFullPattern_simulationMesh->reset();
pSimulationEdgeMesh_tilledPattern->setBeamMaterial(
0.3, optimizationResult.settings.youngsModulus_pattern);
pSimulationEdgeMesh_tilledPattern->reset();
// optimizationResult.draw();
#ifdef POLYSCOPE_DEFINED
pSimulationEdgeMesh_tilledPattern->registerForDrawing(
color_tesselatedPatterns);
#endif
////Tessellated reduced pattern simulation mesh
std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
= perPatternIndexTiledReducedPatternEdgeIndices[0];
std::shared_ptr<SimulationEdgeMesh> pSimulationEdgeMesh_tilledReducedModel;
pSimulationEdgeMesh_tilledReducedModel =
std::make_shared<SimulationEdgeMesh>(*pTilled_reducedModel);
ReducedModelOptimization::Results::applyOptimizationResults_elements(
optimizationResult, pTiledReducedPattern_simulationMesh);
// pTiledReducedPattern_simulationMesh->reset();
optimizationResult, pSimulationEdgeMesh_tilledReducedModel);
pSimulationEdgeMesh_tilledReducedModel->reset();
#ifdef POLYSCOPE_DEFINED
pSimulationEdgeMesh_tilledReducedModel->registerForDrawing(
color_tesselatedReducedModels);
polyscope::show();
#endif
Results evaluationResults;
evaluationResults.distances_drm2reduced.fill(-1);
evaluationResults.distances_normalizedDrm2reduced.fill(-1);
DRMSimulationModel::Settings drmSimulationSettings;
drmSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// drmSimulationSettings.threshold_residualToExternalForcesNorm = 1e-3;
// drmSimulationSettings.load(drmSettingsFilePath);
drmSimulationSettings.beVerbose = true;
drmSimulationSettings.maxDRMIterations = 5e6;
drmSimulationSettings.debugModeStep = 100000;
drmSimulationSettings.translationalKineticEnergyThreshold = 1e-15;
drmSimulationSettings.linearGuessForceScaleFactor = 0.8;
drmSimulationSettings.viscousDampingFactor = 7e-3;
drmSimulationSettings.xi = 0.9999;
// drmSimulationSettings.Dtini = 5.86;
drmSimulationSettings.gamma = 0.25;
drmSimulationSettings.threshold_totalTranslationalKineticEnergy = 1e-14;
// drmSimulationSettings.threshold_currentToFirstPeakTranslationalKineticEnergy
// =
// 1e-10;
drmSimulationSettings.threshold_averageResidualToExternalForcesNorm = 1e-5;
// drmSimulationSettings.linearGuessForceScaleFactor = 0.8;
// drmSimulationSettings.viscousDampingFactor = 7e-3;
// drmSimulationSettings.xi = 0.9999;
// drmSimulationSettings.gamma = 0.25;
#ifdef POLYSCOPE_DEFINED
// drmSimulationSettings.shouldDraw = true;
drmSimulationSettings.shouldCreatePlots = true;
drmSimulationSettings.shouldCreatePlots = false;
constexpr bool shouldDrawScenarioResults = true;
if (shouldDrawScenarioResults) {
pSimulationEdgeMesh_tilledPattern->registerForDrawing(
ReducedModelOptimization::Colors::patternInitial);
}
#endif
constexpr bool shouldRerunFullPatternSimulation = false;
const std::string &simulationModelLabel = optimizationResult.settings.simulationModelLabel;
// for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
const std::string& simulationModelLabel_pattern =
optimizationResult.settings.simulationModelLabel_groundTruth;
const std::string& simulationModelLabel_reducedModel =
optimizationResult.settings.simulationModelLabel_reducedModel;
const bool shouldRerunFullPatternSimulation = [&]() {
// if (simulationModelLabel_pattern == DRMSimulationModel::label) {
return false;
// }
// return true;
}();
std::for_each(
//#ifndef POLYSCOPE_DEFINED
#ifndef POLYSCOPE_DEFINED
std::execution::par_unseq,
//#endif
scenariosTestSetLabels.begin(),
scenariosTestSetLabels.end(),
[&](const std::string &jobLabel) {
//check if reduced model scenario exists
// const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
const std::filesystem::path tiledReducedPatternJobFilePath
= std::filesystem::path(scenariosDirectoryPath)
#endif
scenariosTestSetLabels.begin(), scenariosTestSetLabels.end(),
[&](const std::string& jobLabel) {
// check if reduced model scenario exists
// const std::string &jobLabel =
// scenariosTestSetLabels[jobIndex];
const std::filesystem::path tiledReducedPatternJobFilePath =
std::filesystem::path(scenariosDirectoryPath)
.append(pTileIntoSurface->getLabel())
.append(jobLabel)
.append("ReducedJob")
.append(SimulationJob::jsonDefaultFileName);
if (!std::filesystem::exists(tiledReducedPatternJobFilePath)) {
std::cerr << "Scenario " << jobLabel
<< " not found in:" << tiledReducedPatternJobFilePath << std::endl;
<< " not found in:" << tiledReducedPatternJobFilePath
<< std::endl;
// continue; //if not move on to the next scenario
return;
}
//Map the reduced job to the job on the pattern tessellation
//set jobs
std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
// Map the reduced job to the job on the pattern tessellation
// set jobs
std::shared_ptr<SimulationJob> pJob_tiledReducedModel;
pJob_tiledReducedModel =
std::make_shared<SimulationJob>(SimulationJob());
pJob_tiledReducedModel->load(tiledReducedPatternJobFilePath, false);
pJob_tiledReducedModel->pMesh = pSimulationEdgeMesh_tilledReducedModel;
std::shared_ptr<SimulationJob> pJob_tilledPattern;
pJob_tilledPattern = std::make_shared<SimulationJob>(SimulationJob());
pJob_tilledPattern->pMesh = pSimulationEdgeMesh_tilledPattern;
pJob_tiledReducedModel->remap(reducedToFullViMap, *pJob_tilledPattern);
// pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
// pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
// polyscope::show();
const std::filesystem::path surfaceFolderPath
= std::filesystem::path(fullPatternTessellatedResultsDirectoryPath)
.append(simulationModelLabel + "_" + pTileIntoSurface->getLabel());
const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
const std::filesystem::path scenarioDirectoryPath
= std::filesystem::path(surfaceFolderPath).append(scenarioLabel);
//Save reduced job
const std::filesystem::path surfaceFolderPath =
std::filesystem::path(patternTessellatedResultsDirectoryPath)
.append(simulationModelLabel_pattern + "_" +
pTileIntoSurface->getLabel());
const std::string scenarioLabel = pJob_tilledPattern->getLabel();
const std::filesystem::path scenarioDirectoryPath =
std::filesystem::path(surfaceFolderPath).append(scenarioLabel);
// Save reduced job
constexpr bool exportReducedJob = false;
if (exportReducedJob) {
const std::filesystem::path reducedJobDirectoryPath
= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
const std::filesystem::path reducedJobDirectoryPath =
std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
std::filesystem::create_directories(reducedJobDirectoryPath);
pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
pJob_tiledReducedModel->save(reducedJobDirectoryPath);
}
//Check if the drm simulation of the full pattern has already been computed
// Check if the drm simulation of the full pattern has already been
// computed
////Full
const std::string &patternLabel = [&]() {
const std::string patternLabel = optimizationResult.baseTriangleFullPattern
.getLabel();
const std::string& patternLabel = [&]() {
const std::string patternLabel =
optimizationResult.baseTrianglePattern.getLabel();
if (patternLabel.find("_") == std::string::npos) {
return std::to_string(optimizationResult.baseTriangleFullPattern.EN()) + "_"
+ patternLabel;
return std::to_string(optimizationResult.baseTrianglePattern.EN()) +
"_" + patternLabel;
} else {
return patternLabel;
}
}();
const auto tilledPatternResultsFolderPath = std::filesystem::path(scenarioDirectoryPath)
const auto tilledPatternResultsFolderPath =
std::filesystem::path(scenarioDirectoryPath)
.append(patternLabel)
.append("Results");
if (shouldRerunFullPatternSimulation
&& std::filesystem::exists(tilledPatternResultsFolderPath)) {
if (shouldRerunFullPatternSimulation &&
std::filesystem::exists(tilledPatternResultsFolderPath)) {
std::filesystem::remove_all(tilledPatternResultsFolderPath);
}
const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
scenarioDirectoryPath)
const std::filesystem::path fullPatternJobFolderPath =
std::filesystem::path(scenarioDirectoryPath)
.append(patternLabel)
.append("SimulationJob");
SimulationResults simulationResults_tilledPattern;
if (std::filesystem::exists(tilledPatternResultsFolderPath)) {
//Load full pattern results
// Load full pattern results
assert(std::filesystem::exists(fullPatternJobFolderPath));
simulationResults_tilledPattern.load(tilledPatternResultsFolderPath,
fullPatternJobFolderPath);
//#ifdef POLYSCOPE_DEFINED
// std::array<double, 3> resultsColor({28.0, 99.0, 227.0});
// simulationResults_tiledFullPattern_drm.registerForDrawing(resultsColor);
#ifdef POLYSCOPE_DEFINED
simulationResults_tilledPattern.registerForDrawing(
color_tesselatedPatterns);
// std::ifstream ifs("CameraSettings.json");
// nlohmann::json json;
// ifs >> json;
// polyscope::view::setCameraFromJson(json.dump(), false);
// // polyscope::show();
// const std::string cameraJson = polyscope::view::getCameraJson();
// std::filesystem::path jsonFilePath("CameraSettings.json");
// std::ofstream jsonFile_cameraSettings(jsonFilePath.string());
// jsonFile_cameraSettings << cameraJson;
// jsonFile_cameraSettings.close();
// std::filesystem::create_directories("screenshots");
// const std::string screenshotOutputFilePath
// = (std::filesystem::current_path()
// .append("screenshots")
// .append(optimizationResult.label + "_" + pJob_tiledFullPattern->getLabel()))
// .string()
// + ".png";
// // std::cout << "Saving image to:" << screenshotOutputFilePath << std::endl;
// polyscope::screenshot(screenshotOutputFilePath, false);
// simulationResults_tiledFullPattern_drm.unregister();
//#endif
polyscope::show();
const std::string cameraJson = polyscope::view::getCameraJson();
std::filesystem::path jsonFilePath("CameraSettings.json");
std::ofstream jsonFile_cameraSettings(jsonFilePath.string());
jsonFile_cameraSettings << cameraJson;
jsonFile_cameraSettings.close();
std::filesystem::create_directories("screenshots");
const std::string screenshotOutputFilePath =
(std::filesystem::current_path()
.append("screenshots")
.append(optimizationResult.label + "_" +
pJob_tilledPattern->getLabel()))
.string() +
".png";
std::cout << "Saving image to:" << screenshotOutputFilePath
<< std::endl;
polyscope::screenshot(screenshotOutputFilePath, false);
simulationResults_tilledPattern.unregister();
#endif
simulationResults_tilledPattern.converged = true;
} else {
std::cout << "Tilled pattern simulation results not found in:"
<< tilledPatternResultsFolderPath << std::endl;
//Full
// Full
std::cout << "Executing:" << jobLabel << std::endl;
SimulationModelFactory factory;
std::unique_ptr<SimulationModel> pTilledPatternSimulationModel = factory.create(
simulationModelLabel);
//TODO: since the drm simulation model does not have a common interface with the rest of simulation models I need to cast it in order to pass simulation settings. Fix it by removing the SimulationSettings argument
if (simulationModelLabel == DRMSimulationModel::label) {
simulationResults_tilledPattern = static_cast<DRMSimulationModel *>(
std::unique_ptr<SimulationModel> pTilledPatternSimulationModel =
factory.create(simulationModelLabel_pattern);
// TODO: since the drm simulation model does not have a common
// interface with the rest of simulation models I need to cast it in
// order to pass simulation settings. Fix it by removing the
// SimulationSettings argument
if (simulationModelLabel_pattern == DRMSimulationModel::label) {
simulationResults_tilledPattern =
static_cast<DRMSimulationModel*>(
pTilledPatternSimulationModel.get())
->executeSimulation(pJob_tiledFullPattern,
->executeSimulation(pJob_tilledPattern,
drmSimulationSettings);
} else if (simulationModelLabel == ChronosEulerSimulationModel::label) {
simulationResults_tilledPattern
= pTilledPatternSimulationModel->executeSimulation(pJob_tiledFullPattern);
} else {
std::cerr << "Simulation model used for computing the optimization results was "
"not recognized"
<< std::endl;
std::terminate();
simulationResults_tilledPattern =
pTilledPatternSimulationModel->executeSimulation(
pJob_tilledPattern);
}
}
if (!simulationResults_tilledPattern.converged) {
@ -557,37 +467,65 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
return;
}
std::filesystem::create_directories(tilledPatternResultsFolderPath);
const std::filesystem::path drmResultsOutputPath
= std::filesystem::path(scenarioDirectoryPath).append(patternLabel);
const std::filesystem::path drmResultsOutputPath =
std::filesystem::path(scenarioDirectoryPath).append(patternLabel);
simulationResults_tilledPattern.save(drmResultsOutputPath);
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_tiledReducedModel
= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
// LinearSimulationModel linearSimulationModel;
SimulationModelFactory factory;
std::unique_ptr<SimulationModel> pSimulationModel_tilledReducedModel =
factory.create(simulationModelLabel_reducedModel);
SimulationResults simulationResults_tiledReducedModel =
pSimulationModel_tilledReducedModel->executeSimulation(
pJob_tiledReducedModel);
// ChronosEulerNonLinearSimulationModel
// debug_chronosNonLinearSimulationModel;
// const auto debug_chronosResults =
// debug_chronosNonLinearSimulationModel.executeSimulation(
// pJob_tilledPattern);
// LinearSimulationModel debug_linearSimulationModel;
// const auto debug_linearSimResults =
// debug_linearSimulationModel.executeSimulation(pJob_tilledPattern);
#ifdef POLYSCOPE_DEFINED
simulationResults_tiledReducedModel.registerForDrawing();
simulationResults_tilledPattern.registerForDrawing();
if (shouldDrawScenarioResults) {
simulationResults_tiledReducedModel.registerForDrawing(
ReducedModelOptimization::Colors::reducedDeformed, true,
simulationResults_tilledPattern.pJob->pMesh
->getBeamDimensions()[0]
.getDrawingRadius());
simulationResults_tilledPattern.registerForDrawing(
ReducedModelOptimization::Colors::patternDeformed);
// debug_chronosResults.registerForDrawing();
// debug_linearSimResults.registerForDrawing();
polyscope::show();
// debug_linearSimResults.unregister();
simulationResults_tiledReducedModel.unregister();
simulationResults_tilledPattern.unregister();
// debug_chronosResults.unregister();
}
#endif
//compute the full2reduced distance
const double distance_patternToReduced
= simulationResults_tilledPattern
.computeDistance(simulationResults_tiledReducedModel, fullToReducedViMap);
// compute the full2reduced distance
const double distance_patternToReduced =
simulationResults_tilledPattern.computeDistance(
simulationResults_tiledReducedModel, fullToReducedViMap);
double distance_patternSumOfAllVerts = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
distance_patternSumOfAllVerts += simulationResults_tilledPattern
distance_patternSumOfAllVerts +=
simulationResults_tilledPattern
.displacements[fullToReducedPair.first]
.getTranslation()
.norm();
}
const double distance_normalizedPatternToReduced = distance_patternToReduced
/ distance_patternSumOfAllVerts;
const double distance_normalizedPatternToReduced =
distance_patternToReduced / distance_patternSumOfAllVerts;
const int jobIndex = &jobLabel - &scenariosTestSetLabels[0];
evaluationResults.distances_drm2reduced[jobIndex] = distance_patternToReduced;
evaluationResults.distances_normalizedDrm2reduced[jobIndex]
= distance_normalizedPatternToReduced;
evaluationResults.distances_drm2reduced[jobIndex] =
distance_patternToReduced;
evaluationResults.distances_normalizedDrm2reduced[jobIndex] =
distance_normalizedPatternToReduced;
});
return evaluationResults;

96
reducedmodelevaluator.hpp
View File

@ -2,10 +2,10 @@
#define REDUCEDMODELEVALUATOR_HPP
#include "reducedmodeloptimizer_structs.hpp"
#include "utilities.hpp"
class ReducedModelEvaluator
{
public:
class ReducedModelEvaluator {
public:
enum CSVExportingDirection { Vertical = 0, Horizontal };
enum CSVExportingData {
raw_drm2Reduced = 0,
@ -14,9 +14,9 @@ public:
NumberOfDataTypes
};
inline static std::array<std::string, NumberOfDataTypes>
csvExportingDataStrings{"raw_drm2Reduced", "norm_drm2Reduced", "raw_and_norm_drm2Reduced"};
struct Settings
{
csvExportingDataStrings{"raw_drm2Reduced", "norm_drm2Reduced",
"raw_and_norm_drm2Reduced"};
struct Settings {
CSVExportingDirection exportingDirection{Horizontal};
CSVExportingData exportingData{norm_drm2Reduced};
bool shouldWriteHeader{true};
@ -24,52 +24,48 @@ public:
};
inline static constexpr int NumberOfEvaluationScenarios{22};
struct Results
{
struct Results {
std::array<double, NumberOfEvaluationScenarios> distances_drm2reduced;
std::array<double, NumberOfEvaluationScenarios> distances_normalizedDrm2reduced;
std::array<std::string, NumberOfEvaluationScenarios> evaluationScenarioLabels;
std::array<double, NumberOfEvaluationScenarios>
distances_normalizedDrm2reduced;
std::array<std::string, NumberOfEvaluationScenarios>
evaluationScenarioLabels;
};
ReducedModelEvaluator();
Results evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResult,
const std::filesystem::path &scenariosDirectoryPath,
ReducedModelOptimization::Results& optimizationResult,
const std::filesystem::path& scenariosDirectoryPath,
// const std::filesystem::path &reducedPatternFilePath,
const std::filesystem::path &fullPatternTessellatedResultsDirectoryPath);
Results evaluateReducedModel(ReducedModelOptimization::Results &optimizationResult);
static void printResultsVertically(const ReducedModelEvaluator::Results &evaluationResults,
csvFile &csvOutput);
static void printResults(const ReducedModelEvaluator::Results &evaluationResults,
const std::string &resultsLabel);
const std::filesystem::path& fullPatternTessellatedResultsDirectoryPath);
Results evaluateReducedModel(
ReducedModelOptimization::Results& optimizationResult);
static void printResultsVertically(
const ReducedModelEvaluator::Results& evaluationResults,
csvFile& csvOutput);
static void printResults(
const ReducedModelEvaluator::Results& evaluationResults,
const std::string& resultsLabel);
inline static std::array<std::string, NumberOfEvaluationScenarios> scenariosTestSetLabels{
"22Hex_randomBending0",
"22Hex_randomBending1",
inline static std::array<std::string, NumberOfEvaluationScenarios>
scenariosTestSetLabels{
"22Hex_randomBending0", "22Hex_randomBending1",
"22Hex_randomBending2",
// "22Hex_randomBending3",
"22Hex_randomBending4",
"22Hex_randomBending5",
"22Hex_randomBending4", "22Hex_randomBending5",
// "22Hex_randomBending6",
// "22Hex_randomBending7",
"22Hex_randomBending8",
"22Hex_randomBending9",
"22Hex_randomBending10",
"22Hex_randomBending11",
"22Hex_randomBending8", "22Hex_randomBending9",
"22Hex_randomBending10", "22Hex_randomBending11",
"22Hex_randomBending12",
// "22Hex_randomBending13",
// "22Hex_randomBending14",
// "22Hex_randomBending15",
"22Hex_randomBending16",
"22Hex_randomBending17",
"22Hex_randomBending18",
"22Hex_randomBending19",
"22Hex_randomBending16", "22Hex_randomBending17",
"22Hex_randomBending18", "22Hex_randomBending19",
// "22Hex_randomBending20",
"22Hex_bending_0.005N",
"22Hex_bending_0.01N",
"22Hex_bending_0.03N",
"22Hex_bending_0.005N", "22Hex_bending_0.01N", "22Hex_bending_0.03N",
// "22Hex_bending_0.05N",
"22Hex_pullOppositeVerts_0.05N",
"22Hex_pullOppositeVerts_0.1N",
"22Hex_pullOppositeVerts_0.05N", "22Hex_pullOppositeVerts_0.1N",
// "22Hex_pullOppositeVerts_0.3N",
//#ifdef POLYSCOPE_DEFINED
// "22Hex_shear_2N",
@ -82,22 +78,30 @@ public:
// "notUsed_22Hex_axial_10N",
// "notUsed_22Hex_axial_20N",
//#endif
"22Hex_cylinder_0.05N",
"22Hex_cylinder_0.1N",
"22Hex_s_0.05N",
"22Hex_cylinder_0.05N", "22Hex_cylinder_0.1N", "22Hex_s_0.05N",
// "22Hex_s_0.1N"
};
static void printResultsHorizontally(const Results &evaluationResults, csvFile &csvOutput);
static void printResults(const Results &evaluationResults,
const Settings &settings,
csvFile &csvOutput);
static void printHeader(const Settings &settings, csvFile &csvOutput);
static void printResultsHorizontally(const Results& evaluationResults,
csvFile& csvOutput);
static void printResults(const Results& evaluationResults,
const Settings& settings,
csvFile& csvOutput);
static void printHeader(const Settings& settings, csvFile& csvOutput);
// static double evaluateOptimizationSettings(
// const ReducedModelOptimization::Settings &optimizationSettings,
// const std::vector<std::shared_ptr<PatternGeometry>> &pPatterns,
// std::vector<ReducedModelEvaluator::Results> &patternEvaluationResults);
// std::vector<ReducedModelEvaluator::Results>
// &patternEvaluationResults);
std::shared_ptr<VCGPolyMesh> pTileIntoSurface;
inline static constexpr char *tileIntoSurfaceFileContent = R"~(OFF
ReducedModelOptimization::Colors::RGBColor color_tesselatedPatterns{
24.0 / 255, 23.0 / 255, 23.0 / 255};
ReducedModelOptimization::Colors::RGBColor color_tesselatedReducedModels{
67.0 / 255, 160.00 / 255, 232.0 / 255};
ReducedModelOptimization::Colors::RGBColor color_tileIntoSurface{
222 / 255.0, 235 / 255.0, 255 / 255.0};
ReducedModelOptimization::Colors::RGBColor interfaceNodes_color{
63.0 / 255, 85.0 / 255, 42.0 / 255};
inline static constexpr char* tileIntoSurfaceFileContent = R"~(OFF
46 66 0
-0.0745923 0.03573945 0
-0.07464622 0.02191801 0

2223
reducedmodeloptimizer.cpp
View File

File diff suppressed because it is too large Load Diff

443
reducedmodeloptimizer.hpp
View File

@ -1,18 +1,17 @@
#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP
#include "chronoseulersimulationmodel.hpp"
#include "csvfile.hpp"
//#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "edgemesh.hpp"
#include "linearsimulationmodel.hpp"
#ifdef POLYSCOPE_DEFINED
#include "matplot/matplot.h"
#endif
#include <Eigen/Dense>
#include "reducedmodel.hpp"
#include "reducedmodeloptimizer_structs.hpp"
#include "simulationmesh.hpp"
#include <Eigen/Dense>
#ifdef DLIB_DEFINED
#include <dlib/global_optimization.h>
#include <dlib/optimization.h>
@ -21,85 +20,86 @@
#ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h"
#endif // POLYSCOPE_DEFINED
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
using PatternVertexIndex = VertexIndex;
using ReducedModelVertexIndex = VertexIndex;
class ReducedModelOptimizer
{
public:
struct OptimizationState
{
class ReducedModelOptimizer {
public:
struct OptimizationState {
std::vector<SimulationResults> fullPatternResults;
std::vector<double> translationalDisplacementNormalizationValues;
std::vector<double> rotationalDisplacementNormalizationValues;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
std::vector<std::shared_ptr<SimulationJob>> pSimulationJobs_pattern;
std::vector<std::shared_ptr<SimulationJob>> simulationJobs_reducedModel;
std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>
reducedToFullInterfaceViMap;
std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>
fullPatternOppositeInterfaceViPairs;
matplot::line_handle gPlotHandle;
std::vector<size_t> objectiveValueHistory_iteration;
std::vector<double> objectiveValueHistory;
std::vector<double> plotColors;
std::array<double,
ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
ReducedModelOptimization::OptimizationParameterIndex::
NumberOfOptimizationVariables>
parametersInitialValue;
std::array<double,
ReducedModelOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
ReducedModelOptimization::OptimizationParameterIndex::
NumberOfOptimizationVariables>
optimizationInitialValue;
std::vector<int> simulationScenarioIndices;
double minY{DBL_MAX};
std::vector<double> minX;
int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0};
//Variables for finding the full pattern simulation forces
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
std::array<std::function<void(const double &newValue,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>,
// Variables for finding the full pattern simulation forces
std::shared_ptr<SimulationEdgeMesh> pFullPatternSimulationEdgeMesh;
std::array<std::function<void(const double& newValue,
std::shared_ptr<SimulationEdgeMesh>&
pReducedPatternSimulationEdgeMesh)>,
7>
functions_updateReducedPatternParameter;
std::vector<double> xMin;
std::vector<double> xMax;
std::vector<double> scenarioWeights;
std::vector<ReducedModelOptimization::Settings::ObjectiveWeights> objectiveWeights;
std::vector<ReducedModelOptimization::Settings::ObjectiveWeights>
objectiveWeights;
std::string simulationModelLabel_reducedModel;
};
private:
private:
OptimizationState optimizationState;
vcg::Triangle3<double> baseTriangle;
std::function<void(const double &,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
SimulationJob &)>
std::function<void(
const double&,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&,
SimulationJob&)>
constructScenarioFunction;
std::shared_ptr<SimulationMesh> m_pReducedModelSimulationMesh;
std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
std::shared_ptr<SimulationEdgeMesh> m_pReducedModelSimulationEdgeMesh;
std::shared_ptr<SimulationEdgeMesh> m_pSimulationEdgeMesh_pattern;
std::unordered_map<PatternVertexIndex, ReducedModelVertexIndex>
m_fullToReducedInterfaceViMap;
std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>
m_fullPatternOppositeInterfaceViPairs;
std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
std::string optimizationNotes;
std::array<std::function<void(
const double &,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
SimulationJob &)>,
std::array<
std::function<void(
const double&,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&,
SimulationJob&)>,
ReducedModelOptimization::NumberOfBaseSimulationScenarios>
constructBaseScenarioFunctions;
std::vector<bool> scenarioIsSymmetrical;
int fullPatternNumberOfEdges;
constexpr static double youngsModulus{1 * 1e9};
constexpr static double defaultBeamWidth{0.002};
constexpr static double defaultBeamHeight{0.002};
std::string fullPatternLabel;
// ReducedModelOptimization::Settings optimizationSettings;
public:
struct FunctionEvaluation
{
public:
struct FunctionEvaluation {
FunctionEvaluation() = default;
FunctionEvaluation(const std::vector<double> &x, double y) : x(x), y(y) {}
FunctionEvaluation(const std::vector<double>& x, double y) : x(x), y(y) {}
std::vector<double> x;
double y = std::numeric_limits<double>::quiet_NaN();
@ -117,170 +117,195 @@ public:
// };
// inline constexpr static ParameterLabels parameterLabels();
inline static std::array<std::string, ReducedModelOptimization::NumberOfOptimizationVariables>
inline static std::array<
std::string,
ReducedModelOptimization::NumberOfOptimizationVariables>
parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"};
constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
constexpr static std::
array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
simulationScenariosResolution = {12, 12, 22, 22, 22, 22};
constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
// simulationScenariosResolution = {2, 2, 2, 22, 22, 22};
constexpr static std::
array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
baseScenarioWeights = {1, 1, 2, 2, 2};
inline static int totalNumberOfSimulationScenarios
= std::accumulate(simulationScenariosResolution.begin(),
inline static int totalNumberOfSimulationScenarios =
std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.end(),
0);
inline static int fanSize{6};
inline static double initialHexagonSize{0.3};
inline static int fanCardinality{6};
inline static double initialValue_R{0.3};
inline static double initialValue_theta{0};
int interfaceNodeIndex;
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
double operator()(const Eigen::VectorXd& x, Eigen::VectorXd&) const;
ReducedModelOptimizer();
static void computeReducedModelSimulationJob(
const SimulationJob &simulationJobOfFullModel,
const std::unordered_map<size_t, size_t> &fullToReducedMap,
SimulationJob &simulationJobOfReducedModel);
const SimulationJob& simulationJobOfFullModel,
const std::unordered_map<size_t, size_t>& fullToReducedMap,
SimulationJob& simulationJobOfReducedModel);
SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
SimulationJob getReducedSimulationJob(
const SimulationJob& fullModelSimulationJob);
static void runSimulation(const std::string &filename, std::vector<double> &x);
static void runSimulation(const std::string& filename,
std::vector<double>& x);
static std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);
static std::vector<std::shared_ptr<SimulationJob>>
createFullPatternSimulationJobs(
const std::shared_ptr<SimulationEdgeMesh>& pMesh,
const std::unordered_map<size_t, size_t>&
fullPatternOppositeInterfaceViMap);
static void createSimulationMeshes(
PatternGeometry &pattern,
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions,
std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &fullPattern,
ReducedModel &reducedPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
&fullToReducedInterfaceViMap,
std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
&fullPatternOppositeInterfaceViMap);
static void createSimulationEdgeMeshes(
PatternGeometry& pattern,
PatternGeometry& reducedModel,
const ReducedModelOptimization::Settings& optimizationSettings,
std::shared_ptr<SimulationEdgeMesh>& pFullPatternSimulationEdgeMesh,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh);
void computeMaps(
const std::unordered_map<size_t, std::unordered_set<size_t>>& slotToNode,
PatternGeometry& fullPattern,
ReducedModel& reducedPattern,
std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap,
std::unordered_map<PatternVertexIndex, ReducedModelVertexIndex>&
fullToReducedInterfaceViMap,
std::vector<std::pair<PatternVertexIndex, ReducedModelVertexIndex>>&
fullPatternOppositeInterfaceViMap);
static void visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
const std::vector<std::shared_ptr<SimulationJob>>&
fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>>&
reducedPatternSimulationJobs,
const std::vector<ReducedModelOptimization::BaseSimulationScenario>&
simulationScenarios,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap);
static void registerResultsForDrawing(
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
const std::shared_ptr<SimulationJob>& pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob>& pReducedPatternSimulationJob,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap);
static double computeRawTranslationalError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
const std::vector<Vector6d>& fullPatternDisplacements,
const std::vector<Vector6d>& reducedPatternDisplacements,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap);
static double computeDisplacementError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
const std::vector<Vector6d>& fullPatternDisplacements,
const std::vector<Vector6d>& reducedPatternDisplacements,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap,
const double& normalizationFactor);
static double computeRawRotationalError(
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
const std::vector<Eigen::Quaterniond>& rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond>& rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap);
static double computeRotationalError(const std::vector<Eigen::Quaterniond> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
static double computeRotationalError(
const std::vector<Eigen::Quaterniond>& rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaterniond>& rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap,
const double& normalizationFactor);
static double computeError(
const SimulationResults &simulationResults_fullPattern,
const SimulationResults &simulationResults_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor_translationalDisplacement,
const double &normalizationFactor_rotationalDisplacement,
const double &scenarioWeight,
const ReducedModelOptimization::Settings::ObjectiveWeights &objectiveWeights);
const SimulationResults& simulationResults_fullPattern,
const SimulationResults& simulationResults_reducedPattern,
const std::unordered_map<ReducedModelVertexIndex, PatternVertexIndex>&
reducedToFullInterfaceViMap,
const double& normalizationFactor_translationalDisplacement,
const double& normalizationFactor_rotationalDisplacement,
const double& scenarioWeight,
const ReducedModelOptimization::Settings::ObjectiveWeights&
objectiveWeights);
static void constructAxialSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static void constructShearSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static void constructBendingSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static void constructDomeSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static void constructSaddleSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static void constructSSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job);
static std::function<void(const std::vector<double> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
const double& forceMagnitude,
const std::vector<std::pair<PatternVertexIndex, PatternVertexIndex>>&
oppositeInterfaceViPairs,
SimulationJob& job);
static std::function<void(
const std::vector<double>& x,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern;
static std::function<void(const double &newE,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newE,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_E;
static std::function<void(const double &newA,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newA,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_A;
static std::function<void(const double &newI,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newI,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_I;
static std::function<void(const double &newI2,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newI2,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_I2;
static std::function<void(const double &newI3,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newI3,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_I3;
static std::function<void(const double &newJ,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
static std::function<void(
const double& newJ,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
function_updateReducedPattern_material_J;
static double objective(const std::vector<double> &x);
static double objective(const std::vector<double>& x);
void initializeUpdateReducedPatternFunctions();
// static double objective(const double &xValue);
ReducedModelOptimization::Results optimize(
ConstPatternGeometry &fullPattern,
const ReducedModelOptimization::Settings &optimizationSettings);
ConstPatternGeometry& fullPattern,
const ReducedModelOptimization::Settings& optimizationSettings);
void optimize(ConstPatternGeometry &fullPattern,
ReducedModelOptimization::Settings &optimizationSettings,
ReducedModelOptimization::Results &optimizationResults);
static double objective(const std::vector<double> &x,
ReducedModelOptimizer::OptimizationState &optimizationState);
void optimize(ConstPatternGeometry& fullPattern,
ReducedModelOptimization::Settings& optimizationSettings,
ReducedModelOptimization::Results& optimizationResults);
static double objective(
const std::vector<double>& x,
ReducedModelOptimizer::OptimizationState& optimizationState);
private:
private:
void optimize(
ReducedModelOptimization::Settings &optimizationSettings,
ReducedModelOptimization::Results &results,
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios
= std::vector<ReducedModelOptimization::BaseSimulationScenario>(
ReducedModelOptimization::Settings& optimizationSettings,
ReducedModelOptimization::Results& results,
const std::vector<ReducedModelOptimization::BaseSimulationScenario>&
simulationScenarios =
std::vector<ReducedModelOptimization::BaseSimulationScenario>(
{ReducedModelOptimization::Axial,
ReducedModelOptimization::Shear,
ReducedModelOptimization::Bending,
@ -288,78 +313,88 @@ private:
ReducedModelOptimization::Saddle,
ReducedModelOptimization::S}));
void initializePatterns(PatternGeometry &fullPattern,
ReducedModel &reducedModel,
const ReducedModelOptimization::Settings &optimizationSettings);
void initializePatterns(
PatternGeometry& fullPattern,
ReducedModel& reducedModel,
const ReducedModelOptimization::Settings& optimizationSettings);
static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
void runOptimization(const ReducedModelOptimization::Settings &settings,
ReducedModelOptimization::Results &results);
void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel);
void createSimulationMeshes(PatternGeometry &fullModel,
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions);
const SimulationResults& fullModelResults,
const std::unordered_map<size_t, size_t>& displacementsReducedToFullMap,
Eigen::MatrixX3d& optimalDisplacementsOfReducedModel);
void runOptimization(const ReducedModelOptimization::Settings& settings,
ReducedModelOptimization::Results& results);
void computeMaps(PatternGeometry& fullModel, ReducedModel& reducedModel);
void createSimulationEdgeMeshes(
PatternGeometry& fullModel,
PatternGeometry& reducedModel,
const ReducedModelOptimization::Settings& optimizationSettings);
void initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh,
const std::shared_ptr<SimulationEdgeMesh>& mesh,
const std::array<ReducedModelOptimization::xRange,
ReducedModelOptimization::NumberOfOptimizationVariables>
&optimizationParamters);
ReducedModelOptimization::NumberOfOptimizationVariables>&
optimizationParamters);
ChronosEulerSimulationModel patternSimulator;
LinearSimulationModel reducedModelSimulator;
void computeObjectiveValueNormalizationFactors(
const ReducedModelOptimization::Settings &optimizationSettings);
const ReducedModelOptimization::Settings& optimizationSettings);
void getResults(const FunctionEvaluation &optimalObjective,
const ReducedModelOptimization::Settings &settings,
ReducedModelOptimization::Results &results);
void getResults(const FunctionEvaluation& optimalObjective,
const ReducedModelOptimization::Settings& settings,
ReducedModelOptimization::Results& results);
double computeFullPatternMaxSimulationForce(
const ReducedModelOptimization::BaseSimulationScenario &scenario) const;
const ReducedModelOptimization::BaseSimulationScenario& scenario) const;
#ifdef DLIB_DEFINED
static double objective(const dlib::matrix<double, 0, 1> &x);
static double objective(const dlib::matrix<double, 0, 1>& x);
#endif
std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
computeFullPatternMaxSimulationForces(
const std::vector<ReducedModelOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenario) const;
const std::vector<ReducedModelOptimization::BaseSimulationScenario>&
desiredBaseSimulationScenario) const;
std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
&baseScenarioMaxForceMagnitudes) const;
const std::shared_ptr<SimulationEdgeMesh>& pMesh,
const std::array<
double,
ReducedModelOptimization::NumberOfBaseSimulationScenarios>&
baseScenarioMaxForceMagnitudes) const;
std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
getFullPatternMaxSimulationForces(
const std::vector<ReducedModelOptimization::BaseSimulationScenario>
&desiredBaseSimulationScenarioIndices,
const std::filesystem::path &intermediateResultsDirectoryPath,
const bool &recomputeForceMagnitudes);
const std::vector<ReducedModelOptimization::BaseSimulationScenario>&
desiredBaseSimulationScenarioIndices,
const std::filesystem::path& intermediateResultsDirectoryPath,
const bool& recomputeForceMagnitudes);
std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
getFullPatternMaxSimulationForces();
void computeScenarioWeights(
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &baseSimulationScenarios,
const ReducedModelOptimization::Settings &optimizationSettings);
const std::vector<ReducedModelOptimization::BaseSimulationScenario>&
baseSimulationScenarios,
const ReducedModelOptimization::Settings& optimizationSettings);
};
inline std::function<void(const double &newE,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newE,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_E;
inline std::function<void(const double &newA,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newA,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_A;
inline std::function<void(const double &newI,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newI,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_I;
inline std::function<void(const double &newI2,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newI2,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_I2;
inline std::function<void(const double &newI3,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newI3,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_I3;
inline std::function<void(const double &newJ,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>
inline std::function<void(
const double& newJ,
std::shared_ptr<SimulationEdgeMesh>& pReducedPatternSimulationEdgeMesh)>
ReducedModelOptimizer::function_updateReducedPattern_material_J;
inline std::function<void(const std::vector<double> &x, std::shared_ptr<SimulationMesh> &m)>
inline std::function<void(const std::vector<double>& x,
std::shared_ptr<SimulationEdgeMesh>& m)>
ReducedModelOptimizer::function_updateReducedPattern;
extern ReducedModelOptimizer::OptimizationState global;

1009
reducedmodeloptimizer_structs.hpp
View File

File diff suppressed because it is too large Load Diff

591
topologyenumerator.cpp
View File

@ -1,12 +1,13 @@
#include "topologyenumerator.hpp"
#include <math.h>
#include <algorithm>
#include <boost/graph/biconnected_components.hpp>
#include <iostream>
#include <math.h>
#include <numeric>
#include <thread>
#include <unordered_set>
const bool debugIsOn{false};
const bool debugIsOn{true};
const bool savePlyFiles{true};
// size_t binomialCoefficient(size_t n, size_t m) {
@ -34,42 +35,46 @@ const bool savePlyFiles{true};
// .string());
//}
size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const
{
size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const {
if (ni1 <= ni0) {
std::swap(ni0, ni1);
}
assert(ni1 > ni0);
const size_t &n = numberOfNodes;
const size_t& n = numberOfNodes;
return (n * (n - 1) / 2) - (n - ni0) * ((n - ni0) - 1) / 2 + ni1 - ni0 - 1;
}
TopologyEnumerator::TopologyEnumerator() {}
void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reducedNumberOfNodesPerSlot,
const std::string &desiredResultsPath,
const int &numberOfDesiredEdges)
{
void TopologyEnumerator::computeValidPatterns(
const std::vector<size_t>& reducedNumberOfNodesPerSlot,
const std::string& desiredResultsPath,
const int& numberOfDesiredEdges) {
assert(reducedNumberOfNodesPerSlot.size() == 5);
assert(reducedNumberOfNodesPerSlot[0] == 0 || reducedNumberOfNodesPerSlot[0] == 1);
assert(reducedNumberOfNodesPerSlot[1] == 0 || reducedNumberOfNodesPerSlot[1] == 1);
std::vector<size_t> numberOfNodesPerSlot{reducedNumberOfNodesPerSlot[0],
reducedNumberOfNodesPerSlot[1],
reducedNumberOfNodesPerSlot[1],
reducedNumberOfNodesPerSlot[2],
reducedNumberOfNodesPerSlot[3],
reducedNumberOfNodesPerSlot[2],
assert(reducedNumberOfNodesPerSlot[0] == 0 ||
reducedNumberOfNodesPerSlot[0] == 1);
assert(reducedNumberOfNodesPerSlot[1] == 0 ||
reducedNumberOfNodesPerSlot[1] == 1);
std::vector<size_t> numberOfNodesPerSlot{
reducedNumberOfNodesPerSlot[0], reducedNumberOfNodesPerSlot[1],
reducedNumberOfNodesPerSlot[1], reducedNumberOfNodesPerSlot[2],
reducedNumberOfNodesPerSlot[3], reducedNumberOfNodesPerSlot[2],
reducedNumberOfNodesPerSlot[4]};
// Generate an edge mesh wih all possible edges
numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(), numberOfNodesPerSlot.end(), 0);
const size_t numberOfAllPossibleEdges = numberOfNodes * (numberOfNodes - 1) / 2;
numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(),
numberOfNodesPerSlot.end(), 0);
const size_t numberOfAllPossibleEdges =
numberOfNodes * (numberOfNodes - 1) / 2;
std::vector<vcg::Point2i> allPossibleEdges(numberOfAllPossibleEdges);
const int &n = numberOfNodes;
for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges; edgeIndex++) {
const int ni0 = n - 2
- std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5);
const int ni1 = edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2;
const int& n = numberOfNodes;
for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges;
edgeIndex++) {
const int ni0 =
n - 2 -
std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5);
const int ni1 =
edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2;
allPossibleEdges[edgeIndex] = vcg::Point2i(ni0, ni1);
}
PatternGeometry patternGeometryAllEdges;
@ -88,13 +93,15 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
std::string elemID;
if (numberOfNodesPerSlotIndex == 0 || numberOfNodesPerSlotIndex == 1) {
elemID = "v";
} else if (numberOfNodesPerSlotIndex == 2 || numberOfNodesPerSlotIndex == 3) {
} else if (numberOfNodesPerSlotIndex == 2 ||
numberOfNodesPerSlotIndex == 3) {
elemID = "e";
} else {
elemID = "c";
}
setupString += std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex])
+ elemID + "_";
setupString +=
std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex]) +
elemID + "_";
}
setupString += std::to_string(PatternGeometry().getFanSize()) + "fan";
if (debugIsOn) {
@ -105,34 +112,39 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
std::filesystem::create_directory(resultsPath);
if (debugIsOn) {
patternGeometryAllEdges.save(
std::filesystem::path(resultsPath).append("allPossibleEdges.ply").string());
patternGeometryAllEdges.save(std::filesystem::path(resultsPath)
.append("allPossibleEdges.ply")
.string());
}
// statistics.numberOfPossibleEdges = numberOfAllPossibleEdges;
std::vector<vcg::Point2i> validEdges = getValidEdges(numberOfNodesPerSlot,
resultsPath,
patternGeometryAllEdges,
std::vector<vcg::Point2i> validEdges =
getValidEdges(numberOfNodesPerSlot, resultsPath, patternGeometryAllEdges,
allPossibleEdges);
PatternGeometry patternAllValidEdges;
patternAllValidEdges.add(patternGeometryAllEdges.computeVertices(), validEdges);
patternAllValidEdges.add(patternGeometryAllEdges.computeVertices(),
validEdges);
if (debugIsOn) {
// Export all valid edges in a ply
patternAllValidEdges.save(
std::filesystem::path(resultsPath).append("allValidEdges.ply").string());
patternAllValidEdges.save(std::filesystem::path(resultsPath)
.append("allValidEdges.ply")
.string());
}
// statistics.numberOfValidEdges = validEdges.size();
// Find pairs of intersecting edges
const std::unordered_map<size_t, std::unordered_set<size_t>> intersectingEdges
= patternAllValidEdges.getIntersectingEdges(statistics.numberOfIntersectingEdgePairs);
const std::unordered_map<size_t, std::unordered_set<size_t>>
intersectingEdges = patternAllValidEdges.getIntersectingEdges(
statistics.numberOfIntersectingEdgePairs);
if (debugIsOn) {
auto intersectingEdgesPath = std::filesystem::path(resultsPath)
.append("All_intersecting_edge_pairs");
std::filesystem::create_directory(intersectingEdgesPath);
// Export intersecting pairs in ply files
for (auto mapIt = intersectingEdges.begin(); mapIt != intersectingEdges.end(); mapIt++) {
for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end(); setIt++) {
for (auto mapIt = intersectingEdges.begin();
mapIt != intersectingEdges.end(); mapIt++) {
for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end();
setIt++) {
PatternGeometry intersectingEdgePair;
const size_t ei0 = mapIt->first;
const size_t ei1 = *setIt;
@ -145,15 +157,16 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
patternGeometryAllEdges.computeVertices()[validEdges[ei1][0]],
patternGeometryAllEdges.computeVertices()[validEdges[ei1][1]]);
intersectingEdgePair.save(std::filesystem::path(intersectingEdgesPath)
.append(std::to_string(mapIt->first) + "_"
+ std::to_string(*setIt) + ".ply")
.append(std::to_string(mapIt->first) +
"_" + std::to_string(*setIt) +
".ply")
.string());
}
}
}
const std::unordered_set<VertexIndex> interfaceNodes = patternGeometryAllEdges.getInterfaceNodes(
numberOfNodesPerSlot);
const std::unordered_set<VertexIndex> interfaceNodes =
patternGeometryAllEdges.getInterfaceNodes(numberOfNodesPerSlot);
// assert(validEdges.size() == allPossibleEdges.size() -
// coincideEdges.size() -
@ -174,9 +187,10 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
// std::filesystem::path(resultsPath).append("patterns.patt"));
// }
if (numberOfDesiredEdges == -1) {
for (size_t numberOfEdges = 2; numberOfEdges <= validEdges.size(); numberOfEdges++) {
std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges."
<< std::endl;
for (size_t numberOfEdges = 2; numberOfEdges <= validEdges.size();
numberOfEdges++) {
std::cout << "Computing " + setupString << " with " << numberOfEdges
<< " edges." << std::endl;
auto perEdgeResultPath = std::filesystem::path(resultsPath)
.append(std::to_string(numberOfEdges));
@ -189,19 +203,16 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
// }
}
std::filesystem::create_directory(perEdgeResultPath);
computeValidPatterns(numberOfNodesPerSlot,
numberOfEdges,
computeValidPatterns(numberOfNodesPerSlot, numberOfEdges,
perEdgeResultPath,
patternGeometryAllEdges.computeVertices(),
intersectingEdges,
validEdges,
interfaceNodes);
intersectingEdges, validEdges, interfaceNodes);
statistics.print(setupString, perEdgeResultPath);
statistics.reset();
}
} else {
std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges << " edges."
<< std::endl;
std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges
<< " edges." << std::endl;
auto perEdgeResultPath = std::filesystem::path(resultsPath)
.append(std::to_string(numberOfDesiredEdges));
@ -210,22 +221,19 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
std::filesystem::remove_all(perEdgeResultPath);
}
std::filesystem::create_directory(perEdgeResultPath);
computeValidPatterns(numberOfNodesPerSlot,
numberOfDesiredEdges,
computeValidPatterns(numberOfNodesPerSlot, numberOfDesiredEdges,
perEdgeResultPath,
patternGeometryAllEdges.computeVertices(),
intersectingEdges,
validEdges,
interfaceNodes);
intersectingEdges, validEdges, interfaceNodes);
statistics.print(setupString, perEdgeResultPath);
}
}
void TopologyEnumerator::computeEdgeNodes(const std::vector<size_t> &numberOfNodesPerSlot,
std::vector<size_t> &nodesEdge0,
std::vector<size_t> &nodesEdge1,
std::vector<size_t> &nodesEdge2)
{
void TopologyEnumerator::computeEdgeNodes(
const std::vector<size_t>& numberOfNodesPerSlot,
std::vector<size_t>& nodesEdge0,
std::vector<size_t>& nodesEdge1,
std::vector<size_t>& nodesEdge2) {
// Create vectors holding the node indices of each pattern node of each
// triangle edge
size_t nodeIndex = 0;
@ -238,18 +246,21 @@ void TopologyEnumerator::computeEdgeNodes(const std::vector<size_t> &numberOfNod
nodesEdge2.push_back(nodeIndex++);
if (numberOfNodesPerSlot[3] != 0) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3]; edgeNodeIndex++) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3];
edgeNodeIndex++) {
nodesEdge0.push_back(nodeIndex++);
}
}
if (numberOfNodesPerSlot[4] != 0) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4]; edgeNodeIndex++) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4];
edgeNodeIndex++) {
nodesEdge1.push_back(nodeIndex++);
}
}
if (numberOfNodesPerSlot[5] != 0) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5]; edgeNodeIndex++) {
for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5];
edgeNodeIndex++) {
nodesEdge2.push_back(nodeIndex++);
}
}
@ -265,8 +276,7 @@ void TopologyEnumerator::computeEdgeNodes(const std::vector<size_t> &numberOfNod
}
std::unordered_set<size_t> TopologyEnumerator::computeCoincideEdges(
const std::vector<size_t> &numberOfNodesPerSlot)
{
const std::vector<size_t>& numberOfNodesPerSlot) {
/*
* A coincide edge is defined as an edge connection between two nodes that lay
* on a triangle edge and which have another node in between
@ -279,12 +289,11 @@ std::unordered_set<size_t> TopologyEnumerator::computeCoincideEdges(
std::vector<size_t> coincideEdges0 = getCoincideEdges(nodesEdge0);
std::vector<size_t> coincideEdges1 = getCoincideEdges(nodesEdge1);
std::vector<size_t> coincideEdges2 = getCoincideEdges(nodesEdge2);
std::unordered_set<size_t> coincideEdges{coincideEdges0.begin(), coincideEdges0.end()};
std::copy(coincideEdges1.begin(),
coincideEdges1.end(),
std::unordered_set<size_t> coincideEdges{coincideEdges0.begin(),
coincideEdges0.end()};
std::copy(coincideEdges1.begin(), coincideEdges1.end(),
std::inserter(coincideEdges, coincideEdges.end()));
std::copy(coincideEdges2.begin(),
coincideEdges2.end(),
std::copy(coincideEdges2.begin(), coincideEdges2.end(),
std::inserter(coincideEdges, coincideEdges.end()));
if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[1]) {
@ -300,8 +309,7 @@ std::unordered_set<size_t> TopologyEnumerator::computeCoincideEdges(
}
std::unordered_set<size_t> TopologyEnumerator::computeDuplicateEdges(
const std::vector<size_t> &numberOfNodesPerSlot)
{
const std::vector<size_t>& numberOfNodesPerSlot) {
/*
* A duplicate edges are all edges the "right" edge since due to rotational
* symmetry "left" edge=="right" edge
@ -312,7 +320,8 @@ std::unordered_set<size_t> TopologyEnumerator::computeDuplicateEdges(
std::vector<size_t> nodesEdge2; // right edge
computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2);
if (numberOfNodesPerSlot[5]) {
for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1; edge2NodeIndex++) {
for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1;
edge2NodeIndex++) {
const size_t nodeIndex = nodesEdge2[edge2NodeIndex];
const size_t nextNodeIndex = nodesEdge2[edge2NodeIndex + 1];
duplicateEdges.insert(getEdgeIndex(nodeIndex, nextNodeIndex));
@ -323,18 +332,20 @@ std::unordered_set<size_t> TopologyEnumerator::computeDuplicateEdges(
}
std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
const std::vector<size_t> &numberOfNodesPerSlot,
const std::filesystem::path &resultsPath,
const PatternGeometry &patternGeometryAllEdges,
const std::vector<vcg::Point2i> &allPossibleEdges)
{
std::unordered_set<size_t> coincideEdges = computeCoincideEdges(numberOfNodesPerSlot);
const std::vector<size_t>& numberOfNodesPerSlot,
const std::filesystem::path& resultsPath,
const PatternGeometry& patternGeometryAllEdges,
const std::vector<vcg::Point2i>& allPossibleEdges) {
std::unordered_set<size_t> coincideEdges =
computeCoincideEdges(numberOfNodesPerSlot);
// Export each coincide edge into a ply file
if (!coincideEdges.empty() && debugIsOn) {
auto coincideEdgesPath = std::filesystem::path(resultsPath).append("Coincide_edges");
auto coincideEdgesPath =
std::filesystem::path(resultsPath).append("Coincide_edges");
std::filesystem::create_directories(coincideEdgesPath);
for (auto coincideEdgeIndex : coincideEdges) {
PatternGeometry::EdgeType e = patternGeometryAllEdges.edge[coincideEdgeIndex];
PatternGeometry::EdgeType e =
patternGeometryAllEdges.edge[coincideEdgeIndex];
PatternGeometry singleEdgeMesh;
vcg::Point3d p0 = e.cP(0);
vcg::Point3d p1 = e.cP(1);
@ -345,34 +356,40 @@ std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
singleEdgeMesh.add(std::vector<vcg::Point2i>{vcg::Point2i{0, 1}});
singleEdgeMesh.save(std::filesystem::path(coincideEdgesPath)
.append(std::to_string(coincideEdgeIndex))
.string()
+ ".ply");
.string() +
".ply");
}
}
statistics.numberOfCoincideEdges = coincideEdges.size();
// Compute duplicate edges
std::unordered_set<size_t> duplicateEdges = computeDuplicateEdges(numberOfNodesPerSlot);
std::unordered_set<size_t> duplicateEdges =
computeDuplicateEdges(numberOfNodesPerSlot);
if (!duplicateEdges.empty() && debugIsOn) {
// Export duplicate edges in a single ply file
auto duplicateEdgesPath = std::filesystem::path(resultsPath).append("duplicate");
auto duplicateEdgesPath =
std::filesystem::path(resultsPath).append("duplicate");
std::filesystem::create_directory(duplicateEdgesPath);
PatternGeometry patternDuplicateEdges;
for (auto duplicateEdgeIndex : duplicateEdges) {
PatternGeometry::EdgeType e = patternGeometryAllEdges.edge[duplicateEdgeIndex];
PatternGeometry::EdgeType e =
patternGeometryAllEdges.edge[duplicateEdgeIndex];
vcg::Point3d p0 = e.cP(0);
vcg::Point3d p1 = e.cP(1);
vcg::tri::Allocator<PatternGeometry>::AddEdge(patternDuplicateEdges, p0, p1);
vcg::tri::Allocator<PatternGeometry>::AddEdge(patternDuplicateEdges, p0,
p1);
}
patternDuplicateEdges.save(
std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply").string());
patternDuplicateEdges.save(std::filesystem::path(duplicateEdgesPath)
.append("duplicateEdges.ply")
.string());
}
statistics.numberOfDuplicateEdges = duplicateEdges.size();
// Create the set of all possible edges without coincide and duplicate edges
std::vector<vcg::Point2i> validEdges;
for (size_t edgeIndex = 0; edgeIndex < allPossibleEdges.size(); edgeIndex++) {
if (coincideEdges.count(edgeIndex) == 0 && duplicateEdges.count(edgeIndex) == 0) {
if (coincideEdges.count(edgeIndex) == 0 &&
duplicateEdges.count(edgeIndex) == 0) {
validEdges.push_back(allPossibleEdges[edgeIndex]);
}
}
@ -380,29 +397,32 @@ std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
return validEdges;
}
void TopologyEnumerator::exportPattern(const std::filesystem::path &saveToPath,
PatternGeometry &patternGeometry,
const bool saveTilledPattern) const
{
void TopologyEnumerator::exportPattern(const std::filesystem::path& saveToPath,
PatternGeometry& patternGeometry,
const bool saveTilledPattern) const {
const std::string patternName = patternGeometry.getLabel();
std::filesystem::create_directory(saveToPath);
patternGeometry.save(std::filesystem::path(saveToPath).append(patternName).string() + ".ply");
patternGeometry.save(
std::filesystem::path(saveToPath).append(patternName).string() + ".ply");
if (saveTilledPattern) {
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(patternGeometry);
tiledPatternGeometry.save(
std::filesystem::path(saveToPath).append(patternName + "_tiled").string() + ".ply");
PatternGeometry tiledPatternGeometry =
PatternGeometry::createTile(patternGeometry);
tiledPatternGeometry.save(std::filesystem::path(saveToPath)
.append(patternName + "_tiled")
.string() +
".ply");
}
}
void TopologyEnumerator::computeValidPatterns(
const std::vector<size_t> &numberOfNodesPerSlot,
const size_t &numberOfDesiredEdges,
const std::filesystem::path &resultsPath,
const std::vector<vcg::Point3d> &allVertices,
const std::unordered_map<size_t, std::unordered_set<size_t>> &intersectingEdges,
const std::vector<vcg::Point2i> &validEdges,
const std::unordered_set<VertexIndex> &interfaceNodes)
{
const std::vector<size_t>& numberOfNodesPerSlot,
const size_t& numberOfDesiredEdges,
const std::filesystem::path& resultsPath,
const std::vector<vcg::Point3d>& allVertices,
const std::unordered_map<size_t, std::unordered_set<size_t>>&
intersectingEdges,
const std::vector<vcg::Point2i>& validEdges,
const std::unordered_set<VertexIndex>& interfaceNodes) {
assert(numberOfNodesPerSlot.size() == 7);
// Iterate over all patterns which have numberOfDesiredEdges edges from
// from the validEdges Identify patterns that contain dangling edges
@ -414,27 +434,33 @@ void TopologyEnumerator::computeValidPatterns(
assert(enoughValidEdgesExist);
// Create pattern result paths
const auto validPatternsPath = std::filesystem::path(resultsPath).append("Valid");
const bool validPathCreatedSuccesfully = std::filesystem::create_directories(validPatternsPath);
assert(validPathCreatedSuccesfully && std::filesystem::exists(validPatternsPath));
const auto validPatternsPath =
std::filesystem::path(resultsPath).append("Valid");
const bool validPathCreatedSuccesfully =
std::filesystem::create_directories(validPatternsPath);
assert(validPathCreatedSuccesfully &&
std::filesystem::exists(validPatternsPath));
// std::ofstream validPatternsFileStream;
// validPatternsFileStream.open(
// validPatternsPath.append("patterns.patt").string());
const std::string compressedPatternsFilePath
= std::filesystem::path(validPatternsPath).append("patterns.patt").string();
const std::string compressedPatternsFilePath =
std::filesystem::path(validPatternsPath).append("patterns.patt").string();
PatternIO::PatternSet patternSet;
patternSet.nodes = allVertices;
const int patternSetBufferSize = 10000;
const size_t numberOfPatterns = PatternGeometry::binomialCoefficient(validEdges.size(),
numberOfDesiredEdges);
const size_t numberOfPatterns = PatternGeometry::binomialCoefficient(
validEdges.size(), numberOfDesiredEdges);
statistics.numberOfPatterns = numberOfPatterns;
// Initialize pattern binary representation
std::string patternBinaryRepresentation;
patternBinaryRepresentation = std::string(numberOfDesiredEdges, '1');
patternBinaryRepresentation += std::string(validEdges.size() - numberOfDesiredEdges, '0');
std::sort(patternBinaryRepresentation.begin(), patternBinaryRepresentation.end());
patternBinaryRepresentation +=
std::string(validEdges.size() - numberOfDesiredEdges, '0');
std::sort(patternBinaryRepresentation.begin(),
patternBinaryRepresentation.end());
/*TODO: Performance could be improved by changing the patternGeometry with
* respect to the previous one. Maybe I could xor the binaryRepresentation
* to the previous one.*/
@ -445,8 +471,8 @@ void TopologyEnumerator::computeValidPatterns(
constexpr bool saveCompressedFormat = false;
do {
patternIndex++;
const std::string patternName = std::to_string(numberOfDesiredEdges) + "_"
+ std::to_string(patternIndex);
const std::string patternName = std::to_string(numberOfDesiredEdges) + "_" +
std::to_string(patternIndex);
// std::cout << "Pattern name:" + patternBinaryRepresentation <<
// std::endl; isValidPattern(patternBinaryRepresentation, validEdges,
// numberOfDesiredEdges);
@ -454,7 +480,8 @@ void TopologyEnumerator::computeValidPatterns(
// Compute the pattern edges from the binary representation
std::vector<vcg::Point2i> patternEdges(numberOfDesiredEdges);
size_t patternEdgeIndex = 0;
for (size_t validEdgeIndex = 0; validEdgeIndex < patternBinaryRepresentation.size();
for (size_t validEdgeIndex = 0;
validEdgeIndex < patternBinaryRepresentation.size();
validEdgeIndex++) {
if (patternBinaryRepresentation[validEdgeIndex] == '1') {
assert(patternEdgeIndex < numberOfDesiredEdges);
@ -462,96 +489,124 @@ void TopologyEnumerator::computeValidPatterns(
}
}
// DEBUG: Export only part of the pattern set
// if (statistics.numberOfValidPatterns > 1e4) {
// break;
// }
// const bool patternContainsCentroid =
// std::find_if(patternEdges.begin(), patternEdges.end(),
// [](const vcg::Point2i& edge) {
// if (edge[0] == 0 || edge[1] == 0) {
// return true;
// }
// return false;
// }) != patternEdges.end();
// if (!patternContainsCentroid) {
// continue;
// }
PatternGeometry patternGeometry;
patternGeometry.add(allVertices, patternEdges);
patternGeometry.setLabel(patternName);
#ifdef POLYSCOPE_DEFINED
//1st example
// const bool shouldBreak = patternBinaryRepresentation == "00100000100100000"; //398
// const bool shouldBreak = patternBinaryRepresentation == "10000010101110110";//13036
// const bool shouldBreak = patternBinaryRepresentation == "00010111000010100"; //2481
// const bool shouldBreak = patternBinaryRepresentation == "10000101100110010"; //12116
// const bool shouldBreak = patternBinaryRepresentation == "10010111000000110"; //13915
//2nd example
// const bool shouldBreak = patternBinaryRepresentation == "00001011100010011"; //7_1203
// const bool shouldBreak = patternBinaryRepresentation == "00110001100100111"; //4865
// const bool shouldBreak = patternBinaryRepresentation == "00010000101000110"; //1380
// const bool shouldBreak = patternBinaryRepresentation == "00000010100010111"; //268
//3rd
// const bool shouldBreak = patternBinaryRepresentation == "10011011100000010"; //14272
// const bool shouldBreak = patternBinaryRepresentation == "10000111100110110"; //11877
// const bool shouldBreak = patternBinaryRepresentation == "00001011100010011"; //1203
// const bool shouldBreak = patternBinaryRepresentation == "00010101000110000"; //12117
// 1st example
// const bool shouldBreak =
// patternIndex == 1970494; // 398 const bool shouldBreak =
// patternBinaryRepresentation == "10000010101110110";//13036 const
// bool shouldBreak = patternBinaryRepresentation ==
// "00010111000010100"; //2481 const bool shouldBreak =
// patternBinaryRepresentation == "10000101100110010"; //12116 const
// bool shouldBreak = patternBinaryRepresentation ==
// "10010111000000110"; //13915
// 2nd example
// const bool shouldBreak = patternBinaryRepresentation ==
// "00001011100010011"; //7_1203 const bool shouldBreak =
// patternBinaryRepresentation == "00110001100100111"; //4865 const
// bool shouldBreak = patternBinaryRepresentation ==
// "00010000101000110"; //1380 const bool shouldBreak =
// patternBinaryRepresentation == "00000010100010111"; //268
// 3rd
// const bool shouldBreak = patternBinaryRepresentation ==
// "10011011100000010"; //14272 const bool shouldBreak =
// patternBinaryRepresentation == "10000111100110110"; //11877 const
// bool shouldBreak = patternBinaryRepresentation ==
// "00001011100010011"; //1203 const bool shouldBreak =
// patternBinaryRepresentation == "00010101000110000"; //12117
// const bool shouldBreak = patternBinaryRepresentation == "10000101100110100"; //12117
// const bool shouldBreak = patternBinaryRepresentation ==
// "10000101100110100"; //12117
// std::thread drawingThread([&]() {
// if (shouldBreak) {
// patternGeometry.registerForDrawing();
// polyscope::show();
// patternGeometry.unregister();
// }
// });
#endif
// Check if pattern contains intersecting edges
const bool isInterfaceConnected = patternGeometry.isInterfaceConnected(interfaceNodes);
const bool isInterfaceConnected =
patternGeometry.isInterfaceConnected(interfaceNodes);
// Export the tiled ply file if it contains intersecting edges
if (!isInterfaceConnected) {
// create the tiled geometry of the pattern
statistics.numberOfPatternViolatingInterfaceEnforcement++;
if (debugIsOn) {
if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath).append("InterfaceEnforcement"),
patternGeometry,
exportTilledPattern);
}
} else {
continue; // should be uncommented in order to improve performance
}
}
// Check if pattern contains intersecting edges
const bool patternContainsIntersectingEdges
= patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, intersectingEdges);
// Export the tiled ply file if it contains intersecting edges
if (patternContainsIntersectingEdges) {
// create the tiled geometry of the pattern
statistics.numberOfPatternsWithIntersectingEdges++;
if (debugIsOn) {
if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath).append("Intersecting"),
patternGeometry,
exportTilledPattern);
}
} else {
continue; // should be uncommented in order to improve performance
}
}
// const bool shouldBreak = numberOfDesiredEdges == 4 && patternIndex == 53;
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
= patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15);
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
statistics.numberOfPatternsViolatingAngleThreshold++;
if (debugIsOn /*|| savePlyFiles*/) {
if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath)
.append("ExceedingAngleThreshold"),
patternGeometry,
exportTilledPattern);
exportPattern(
std::filesystem::path(resultsPath).append("InterfaceEnforcement"),
patternGeometry, exportTilledPattern);
}
} else {
continue;
}
}
const bool tiledPatternHasNodeWithValenceGreaterThanDesired
= patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6);
// Check if pattern contains intersecting edges
const bool patternContainsIntersectingEdges =
patternGeometry.hasIntersectingEdges(patternBinaryRepresentation,
intersectingEdges);
// Export the tiled ply file if it contains intersecting edges
if (patternContainsIntersectingEdges) {
// create the tiled geometry of the pattern
statistics.numberOfPatternsWithIntersectingEdges++;
if (debugIsOn) {
if (savePlyFiles) {
exportPattern(
std::filesystem::path(resultsPath).append("Intersecting"),
patternGeometry, exportTilledPattern);
}
} else {
continue;
}
}
// const bool shouldBreak = numberOfDesiredEdges == 4 && patternIndex
// == 53;
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold =
patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15);
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
statistics.numberOfPatternsViolatingAngleThreshold++;
if (debugIsOn /*|| savePlyFiles*/) {
if (savePlyFiles) {
exportPattern(std::filesystem::path(resultsPath)
.append("ExceedingAngleThreshold"),
patternGeometry, exportTilledPattern);
}
} else {
continue;
}
}
const bool tiledPatternHasNodeWithValenceGreaterThanDesired =
patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6);
if (tiledPatternHasNodeWithValenceGreaterThanDesired) {
statistics.numberOfPatternsViolatingValenceThreshold++;
if (debugIsOn) {
if (savePlyFiles) {
auto highValencePath = std::filesystem::path(resultsPath)
.append("HighValencePatterns");
auto highValencePath =
std::filesystem::path(resultsPath).append("HighValencePatterns");
exportPattern(highValencePath, patternGeometry, exportTilledPattern);
}
} else {
@ -563,40 +618,47 @@ void TopologyEnumerator::computeValidPatterns(
const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges(
numberOfNodesPerSlot); // marks the nodes with valence>=1
// Create the tiled geometry of the pattern
const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenFanned();
const bool hasFloatingComponents =
!patternGeometry.isFullyConnectedWhenFanned();
PatternGeometry fanPatternGeometry = PatternGeometry::createFan(patternGeometry);
PatternGeometry fanPatternGeometry =
PatternGeometry::createFan(patternGeometry);
const int interfaceNodeVi = 3;
std::vector<PatternGeometry::EdgeType *> connectedEdges;
vcg::edge::VEStarVE(&fanPatternGeometry.vert[interfaceNodeVi], connectedEdges);
std::vector<PatternGeometry::EdgeType*> connectedEdges;
vcg::edge::VEStarVE(&fanPatternGeometry.vert[interfaceNodeVi],
connectedEdges);
if (!connectedEdges.empty()) {
for (int i = 1; i < 6; i++) {
vcg::tri::Allocator<PatternGeometry>::AddEdge(fanPatternGeometry,
interfaceNodeVi
+ (i - 1) * patternGeometry.VN(),
interfaceNodeVi
+ i * patternGeometry.VN());
vcg::tri::Allocator<PatternGeometry>::AddEdge(
fanPatternGeometry,
interfaceNodeVi + (i - 1) * patternGeometry.VN(),
interfaceNodeVi + i * patternGeometry.VN());
}
}
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanPatternGeometry, 0.0000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(fanPatternGeometry);
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(fanPatternGeometry,
0.0000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(
fanPatternGeometry);
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(fanPatternGeometry);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(fanPatternGeometry);
// for (PatternGeometry::VertexType &v : tilledPatternGeometry.vert) {
// for (PatternGeometry::VertexType &v : tilledPatternGeometry.vert)
// {
// std::vector<PatternGeometry::EdgeType *> connectedEdges;
// vcg::edge::VEStarVE(&v, connectedEdges);
// if (connectedEdges.size() == 1) {
// vcg::tri::Allocator<PatternGeometry>::DeleteVertex(tilledPatternGeometry, v);
// vcg::tri::Allocator<PatternGeometry>::DeleteVertex(tilledPatternGeometry,
// v);
// vcg::tri::Allocator<PatternGeometry>::DeleteEdge(tilledPatternGeometry,
// *connectedEdges[0]);
// }
// }
// // vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tilledPatternGeometry);
// //
// vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tilledPatternGeometry);
// fanPatternGeometry.updateEigenEdgeAndVertices();
BoostGraph fanPatternGraph(fanPatternGeometry.VN());
// std::cout << "Edges:";
for (const PatternGeometry::EdgeType &e : fanPatternGeometry.edge) {
for (const PatternGeometry::EdgeType& e : fanPatternGeometry.edge) {
if (e.IsD() || e.cV(0)->IsD() || e.cV(1)->IsD()) {
continue;
}
@ -608,21 +670,28 @@ void TopologyEnumerator::computeValidPatterns(
// std::cout << std::endl;
std::vector<vertex_t> articulationPoints;
boost::articulation_points(fanPatternGraph, std::back_inserter(articulationPoints));
boost::articulation_points(fanPatternGraph,
std::back_inserter(articulationPoints));
const bool hasArticulationPoints = !articulationPoints.empty();
// if (!hasArticulationPoints && tiledPatternHasDanglingEdges) {
// PatternGeometry tilledPatternGeometry = PatternGeometry::createTile(patternGeometry);
// PatternGeometry tilledPatternGeometry =
// PatternGeometry::createTile(patternGeometry);
// tilledPatternGeometry.updateEigenEdgeAndVertices();
// tilledPatternGeometry.registerForDrawing();
// // ->addNodeColorQuantity("de_noAp_tilled", fanVertexColors)
// // ->addNodeColorQuantity("de_noAp_tilled",
// fanVertexColors)
// // ->setEnabled(true);
// polyscope::show();
// tilledPatternGeometry.unregister();
// }
// if (hasArticulationPoints && !tiledPatternHasDanglingEdges/*&& !patternContainsIntersectingEdges
// if (hasArticulationPoints && !tiledPatternHasDanglingEdges/*&&
// !patternContainsIntersectingEdges
// && !hasFloatingComponents
// && !tiledPatternHasNodeWithValenceGreaterThanDesired
// && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) {
// &&
// !tiledPatternHasNodeWithValenceGreaterThanDesired
// &&
// !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/)
// {
// for (PatternGeometry::VertexType &v : patternGeometry.vert) {
// v.C() = vcg::Color4b::Yellow;
// }
@ -632,15 +701,20 @@ void TopologyEnumerator::computeValidPatterns(
// continue;
// }
// // std::cout << articulationPointVi << " ";
// patternGeometry.vert[articulationPointVi].C() = vcg::Color4b::Red;
// patternGeometry.vert[articulationPointVi].C() =
// vcg::Color4b::Red;
// }
// PatternGeometry tilledPatternGeometry = PatternGeometry::createTile(patternGeometry);
// PatternGeometry tilledPatternGeometry =
// PatternGeometry::createTile(patternGeometry);
// // std::cout << std::endl;
// std::vector<glm::vec3> fanVertexColors(tilledPatternGeometry.VN(), glm::vec3(0, 0, 1));
// for (const PatternGeometry::VertexType &v : tilledPatternGeometry.vert) {
// const auto vColor = glm::vec3(v.cC()[0] / 255, v.cC()[1] / 255, v.cC()[2] / 255);
// const auto vi = tilledPatternGeometry.getIndex(v);
// fanVertexColors[vi] = vColor;
// std::vector<glm::vec3>
// fanVertexColors(tilledPatternGeometry.VN(), glm::vec3(0, 0,
// 1)); for (const PatternGeometry::VertexType &v :
// tilledPatternGeometry.vert) {
// const auto vColor = glm::vec3(v.cC()[0] / 255, v.cC()[1] /
// 255, v.cC()[2] / 255); const auto vi =
// tilledPatternGeometry.getIndex(v); fanVertexColors[vi] =
// vColor;
// }
// tilledPatternGeometry.updateEigenEdgeAndVertices();
// tilledPatternGeometry.registerForDrawing()
@ -657,8 +731,10 @@ void TopologyEnumerator::computeValidPatterns(
statistics.numberOfPatternsWithADanglingEdgeOrNode++;
if (debugIsOn) {
if (savePlyFiles) {
auto danglingEdgesPath = std::filesystem::path(resultsPath).append("Dangling");
exportPattern(danglingEdgesPath, patternGeometry, exportTilledPattern);
auto danglingEdgesPath =
std::filesystem::path(resultsPath).append("Dangling");
exportPattern(danglingEdgesPath, patternGeometry,
exportTilledPattern);
}
} else {
continue;
@ -669,28 +745,29 @@ void TopologyEnumerator::computeValidPatterns(
statistics.numberOfPatternsWithMoreThanASingleCC++;
if (debugIsOn) {
if (savePlyFiles) {
auto moreThanOneCCPath = std::filesystem::path(resultsPath)
.append("MoreThanOneCC");
auto moreThanOneCCPath =
std::filesystem::path(resultsPath).append("MoreThanOneCC");
std::filesystem::create_directory(moreThanOneCCPath);
patternGeometry.save(
std::filesystem::path(moreThanOneCCPath).append(patternName).string()
+ ".ply");
patternGeometry.save(std::filesystem::path(moreThanOneCCPath)
.append(patternName)
.string() +
".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
std::vector<std::pair<int, PatternGeometry::EdgePointer>> eCC;
vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(tiledPatternGeometry,
eCC);
vcg::tri::UpdateFlags<PatternGeometry>::EdgeClear(tiledPatternGeometry);
vcg::tri::Clean<PatternGeometry>::edgeMeshConnectedComponents(
tiledPatternGeometry, eCC);
vcg::tri::UpdateFlags<PatternGeometry>::EdgeClear(
tiledPatternGeometry);
const size_t numberOfCC_edgeBased = eCC.size();
std::sort(eCC.begin(),
eCC.end(),
[](const std::pair<int, PatternGeometry::EdgePointer> &a,
const std::pair<int, PatternGeometry::EdgePointer> &b) {
std::sort(eCC.begin(), eCC.end(),
[](const std::pair<int, PatternGeometry::EdgePointer>& a,
const std::pair<int, PatternGeometry::EdgePointer>& b) {
return a.first > b.first;
});
PatternGeometry::EdgePointer &ep = eCC[0].second;
PatternGeometry::EdgePointer& ep = eCC[0].second;
size_t colorsRegistered = 0;
std::stack<EdgePointer> stack;
stack.push(ep);
@ -706,12 +783,10 @@ void TopologyEnumerator::computeValidPatterns(
stack.push(vei.E());
tiledPatternGeometry
.vert[tiledPatternGeometry.getIndex(vei.V1())]
.C()
= vcg::Color4b::Blue;
.C() = vcg::Color4b::Blue;
tiledPatternGeometry
.vert[tiledPatternGeometry.getIndex(vei.V0())]
.C()
= vcg::Color4b::Blue;
.C() = vcg::Color4b::Blue;
colorsRegistered++;
}
++vei;
@ -724,8 +799,8 @@ void TopologyEnumerator::computeValidPatterns(
if (exportTilledPattern) {
tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
.string() +
".ply");
}
}
} else {
@ -737,35 +812,40 @@ void TopologyEnumerator::computeValidPatterns(
statistics.numberOfPatternsWithArticulationPoints++;
if (debugIsOn) {
if (savePlyFiles) {
auto articulationPointsPath = std::filesystem::path(resultsPath)
.append("ArticulationPoints");
exportPattern(articulationPointsPath, patternGeometry, exportTilledPattern);
auto articulationPointsPath =
std::filesystem::path(resultsPath).append("ArticulationPoints");
exportPattern(articulationPointsPath, patternGeometry,
exportTilledPattern);
}
} else {
continue;
}
}
const bool isValidPattern = !patternContainsIntersectingEdges
&& isInterfaceConnected
const bool isValidPattern =
!patternContainsIntersectingEdges &&
isInterfaceConnected
/*&& !tiledPatternHasDanglingEdges*/
&& !hasFloatingComponents && !hasArticulationPoints
&& !tiledPatternHasNodeWithValenceGreaterThanDesired
&& !tiledPatternHasEdgesWithAngleSmallerThanThreshold;
&& !hasFloatingComponents && !hasArticulationPoints &&
!tiledPatternHasNodeWithValenceGreaterThanDesired &&
!tiledPatternHasEdgesWithAngleSmallerThanThreshold;
// if (!hasArticulationPoints && !patternContainsIntersectingEdges
// && !tiledPatternHasDanglingEdges && !hasFloatingComponents
// && !tiledPatternHasNodeWithValenceGreaterThanDesired
// && tiledPatternHasEdgesWithAngleSmallerThanThreshold && numberOfDesiredEdges > 4) {
// std::cout << "Pattern found:" << patternName << std::endl;
// && tiledPatternHasEdgesWithAngleSmallerThanThreshold &&
// numberOfDesiredEdges > 4) { std::cout << "Pattern found:" <<
// patternName << std::endl;
// patternGeometry.registerForDrawing();
// polyscope::show();
// patternGeometry.unregister();
// }
if (isValidPattern) {
// if(patternName=='2055'){
// PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
// patternGeometry); // the marked nodes of hasDanglingEdges are
// tiledPatternGeometry.registerForDrawing(std::array<double, 3>{0, 0, 1});
// PatternGeometry tiledPatternGeometry =
// PatternGeometry::createTile(
// patternGeometry); // the marked nodes of
// hasDanglingEdges are
// tiledPatternGeometry.registerForDrawing(RGBColor{0, 0, 1});
// polyscope::show();
// tiledPatternGeometry.unregister();
// }
@ -789,6 +869,10 @@ void TopologyEnumerator::computeValidPatterns(
}
}
// if (drawingThread.joinable()) {
// drawingThread.join();
// }
// assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges);
} while (std::next_permutation(patternBinaryRepresentation.begin(),
patternBinaryRepresentation.end()));
@ -805,26 +889,25 @@ void TopologyEnumerator::computeValidPatterns(
}
std::vector<size_t> TopologyEnumerator::getCoincideEdges(
const std::vector<size_t> &edgeNodeIndices) const
{
const std::vector<size_t>& edgeNodeIndices) const {
std::vector<size_t> coincideEdges;
if (edgeNodeIndices.size() < 3)
return coincideEdges;
for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2; edgeNodeIndex++) {
const size_t &firstNodeIndex = edgeNodeIndices[edgeNodeIndex];
for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2;
edgeNodeIndex++) {
const size_t& firstNodeIndex = edgeNodeIndices[edgeNodeIndex];
for (size_t secondEdgeNodeIndex = edgeNodeIndex + 2;
secondEdgeNodeIndex < edgeNodeIndices.size();
secondEdgeNodeIndex++) {
const size_t &secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex];
secondEdgeNodeIndex < edgeNodeIndices.size(); secondEdgeNodeIndex++) {
const size_t& secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex];
coincideEdges.push_back(getEdgeIndex(firstNodeIndex, secondNodeIndex));
}
}
return coincideEdges;
}
bool TopologyEnumerator::isValidPattern(const std::string &patternBinaryRepresentation,
const std::vector<vcg::Point2i> &validEdges,
const size_t &numberOfDesiredEdges) const
{
bool TopologyEnumerator::isValidPattern(
const std::string& patternBinaryRepresentation,
const std::vector<vcg::Point2i>& validEdges,
const size_t& numberOfDesiredEdges) const {
return true;
}

583
trianglepatterngeometry.cpp
View File

File diff suppressed because it is too large Load Diff

144
trianglepatterngeometry.hpp
View File

@ -1,23 +1,22 @@
#ifndef FLATPATTERNGEOMETRY_HPP
#define FLATPATTERNGEOMETRY_HPP
#include "edgemesh.hpp"
#include <string>
#include <unordered_map>
#include <unordered_set>
#include "vcgtrimesh.hpp"
#include "edgemesh.hpp"
#include "polymesh.hpp"
#include "vcgtrimesh.hpp"
class PatternGeometry;
using ConstPatternGeometry = PatternGeometry;
class PatternGeometry : public VCGEdgeMesh
{
private:
size_t
computeTiledValence(const size_t &nodeIndex,
const std::vector<size_t> &numberOfNodesPerSlot) const;
size_t getNodeValence(const size_t &nodeIndex) const;
size_t getNodeSlot(const size_t &nodeIndex) const;
class PatternGeometry : public VCGEdgeMesh {
private:
size_t computeTiledValence(
const size_t& nodeIndex,
const std::vector<size_t>& numberOfNodesPerSlot) const;
size_t getNodeValence(const size_t& nodeIndex) const;
size_t getNodeSlot(const size_t& nodeIndex) const;
void addNormals();
double baseTriangleHeight;
@ -30,13 +29,15 @@ private:
std::unordered_map<size_t, size_t> nodeTiledValence;
vcg::Triangle3<double> baseTriangle;
void
constructCorrespondingNodeMap(const std::vector<size_t> &numberOfNodesPerSlot);
void constructCorrespondingNodeMap(
const std::vector<size_t>& numberOfNodesPerSlot);
void constructNodeToTiledValenceMap(const std::vector<size_t> &numberOfNodesPerSlot);
void constructNodeToTiledValenceMap(
const std::vector<size_t>& numberOfNodesPerSlot);
std::vector<VectorType> getEdgeVectorsWithVertexAsOrigin(std::vector<EdgePointer> &edgePointers,
const int &vi);
std::vector<VectorType> getEdgeVectorsWithVertexAsOrigin(
std::vector<EdgePointer>& edgePointers,
const int& vi);
public:
inline static VectorType DefaultNormal{0.0, 0.0, 1.0};
@ -45,56 +46,58 @@ private:
* a const argument but this is impossible due to the
* vcglib interface.
* */
PatternGeometry(PatternGeometry &other);
bool load(const std::filesystem::path &meshFilePath) override;
void add(const std::vector<vcg::Point3d> &vertices);
void add(const std::vector<vcg::Point2i> &edges);
void add(const std::vector<vcg::Point3d> &vertices, const std::vector<vcg::Point2i> &edges);
void add(const std::vector<size_t> &numberOfNodesPerSlot, const std::vector<vcg::Point2i> &edges);
PatternGeometry(PatternGeometry& other);
bool load(const std::filesystem::path& meshFilePath) override;
void add(const std::vector<vcg::Point3d>& vertices);
void add(const std::vector<vcg::Point2i>& edges);
void add(const std::vector<vcg::Point3d>& vertices,
const std::vector<vcg::Point2i>& edges);
void add(const std::vector<size_t>& numberOfNodesPerSlot,
const std::vector<vcg::Point2i>& edges);
static std::vector<vcg::Point3d> constructVertexVector(
const std::vector<size_t> &numberOfNodesPerSlot,
const size_t &fanSize,
const double &triangleEdgeSize);
bool hasDanglingEdges(const std::vector<size_t> &numberOfNodesPerSlot);
bool hasValenceGreaterThan(const std::vector<size_t> &numberOfNodesPerSlot,
const size_t &valenceThreshold);
const std::vector<size_t>& numberOfNodesPerSlot,
const size_t& fanSize,
const double& triangleEdgeSize);
bool hasDanglingEdges(const std::vector<size_t>& numberOfNodesPerSlot);
bool hasValenceGreaterThan(const std::vector<size_t>& numberOfNodesPerSlot,
const size_t& valenceThreshold);
std::vector<vcg::Point3d> computeVertices() const;
static PatternGeometry createFan(PatternGeometry &pattern);
static PatternGeometry createTile(PatternGeometry &pattern);
static PatternGeometry createFan(PatternGeometry& pattern);
static PatternGeometry createTile(PatternGeometry& pattern);
double getTriangleEdgeSize() const;
bool hasUntiledDanglingEdges();
std::unordered_map<size_t, std::unordered_set<size_t>> getIntersectingEdges(
size_t &numberOfIntersectingEdgePairs) const;
size_t& numberOfIntersectingEdgePairs) const;
static size_t binomialCoefficient(size_t n, size_t m)
{
static size_t binomialCoefficient(size_t n, size_t m) {
assert(n >= m);
return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1));
}
bool isFullyConnectedWhenFanned();
bool hasIntersectingEdges(
const std::string &patternBinaryRepresentation,
const std::unordered_map<size_t, std::unordered_set<size_t>>
&intersectingEdges);
const std::string& patternBinaryRepresentation,
const std::unordered_map<size_t, std::unordered_set<size_t>>&
intersectingEdges);
bool isPatternValid();
size_t getFanSize() const;
void add(const std::vector<vcg::Point2d> &vertices,
const std::vector<vcg::Point2i> &edges);
void add(const std::vector<vcg::Point2d>& vertices,
const std::vector<vcg::Point2i>& edges);
PatternGeometry(const std::vector<size_t> &numberOfNodesPerSlot,
const std::vector<vcg::Point2i> &edges);
PatternGeometry(const std::filesystem::path &patternFilePath, bool addNormalsIfAbsent = true);
PatternGeometry(const std::vector<size_t>& numberOfNodesPerSlot,
const std::vector<vcg::Point2i>& edges);
PatternGeometry(const std::filesystem::path& patternFilePath,
bool addNormalsIfAbsent = true);
bool createHoneycombAtom();
void copy(PatternGeometry &copyFrom);
void copy(PatternGeometry& copyFrom);
void tilePattern(VCGEdgeMesh &pattern,
VCGPolyMesh &tileInto,
const int &interfaceNodeIndex,
const bool &shouldDeleteDanglingEdges);
void tilePattern(VCGEdgeMesh& pattern,
VCGPolyMesh& tileInto,
const int& interfaceNodeIndex,
const bool& shouldDeleteDanglingEdges);
void scale(const double &desiredBaseTriangleCentralEdgeSize);
void scale(const double& desiredBaseTriangleCentralEdgeSize);
double getBaseTriangleHeight() const;
vcg::Triangle3<double> computeBaseTriangle() const;
@ -102,35 +105,44 @@ private:
double computeBaseTriangleHeight() const;
void updateBaseTriangleHeight();
// PatternGeometry(const std::vector<vcg::Point2d> &vertices, const std::vector<vcg::Point2i> &edges);
// static std::shared_ptr<PatternGeometry> tilePattern(
// PatternGeometry(const std::vector<vcg::Point2d> &vertices, const
// std::vector<vcg::Point2i> &edges); static std::shared_ptr<PatternGeometry>
// tilePattern(
// std::vector<PatternGeometry> &pattern,
// const std::vector<int> &connectToNeighborsVi,
// const VCGPolyMesh &tileInto,
// std::vector<int> &tileIntoEdgesToTiledVi,
// std::vector<std::vector<size_t>> &perPatternEdges);
virtual void createFan(/*const std::vector<int> &connectToNeighborsVi = std::vector<int>(),*/
const size_t &fanSize = 6);
int interfaceNodeIndex{3}; //TODO: Fix this. This should be automatically computed
bool hasAngleSmallerThanThreshold(const std::vector<size_t> &numberOfNodesPerSlot,
const double &angleThreshold_degrees);
virtual void createFan(/*const std::vector<int> &connectToNeighborsVi =
std::vector<int>(),*/
const size_t& fanSize = 6);
int interfaceNodeIndex{
3}; // TODO: Fix this. This should be automatically computed
bool hasAngleSmallerThanThreshold(
const std::vector<size_t>& numberOfNodesPerSlot,
const double& angleThreshold_degrees,
const bool shouldBreak = false);
vcg::Triangle3<double> getBaseTriangle() const;
static std::shared_ptr<PatternGeometry> tilePattern(PatternGeometry &pattern,
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
std::vector<int> &tileIntoEdgesToTiledVi);
static std::shared_ptr<PatternGeometry> tilePattern(
std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
const std::vector<int> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
std::unordered_set<VertexIndex> getInterfaceNodes(const std::vector<size_t> &numberOfNodesPerSlot);
bool isInterfaceConnected(const std::unordered_set<VertexIndex> &interfaceNodes);
PatternGeometry& pattern,
const std::vector<int>& connectToNeighborsVi,
const VCGPolyMesh& tileInto,
std::vector<int>& tileIntoEdgesToTiledVi);
static std::shared_ptr<PatternGeometry> tilePattern(
std::vector<ConstPatternGeometry>& patterns,
const std::vector<int>& connectToNeighborsVi,
const VCGPolyMesh& tileInto,
const std::vector<int>& perSurfaceFacePatternIndices,
std::vector<size_t>& tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>>& perPatternIndexTiledPatternEdgeIndex);
std::unordered_set<VertexIndex> getInterfaceNodes(
const std::vector<size_t>& numberOfNodesPerSlot);
bool isInterfaceConnected(
const std::unordered_set<VertexIndex>& interfaceNodes);
void deleteDanglingVertices() override;
void deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu) override;
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer>& pu)
override;
};
#endif // FLATPATTERNGEOMETRY_HPP

3
vcgtrimesh.cpp
View File

@ -121,8 +121,7 @@ bool VCGTriMesh::save(const std::string plyFilename)
#ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor, const bool &shouldEnable) const
polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing(const std::optional<std::array<float, 3> > &desiredColor, const bool &shouldEnable) const
{
auto vertices = getVertices();
auto faces = getFaces();

38
vcgtrimesh.hpp
View File

@ -1,8 +1,8 @@
#ifndef VCGTRIMESH_HPP
#define VCGTRIMESH_HPP
#include "mesh.hpp"
#include <vcg/complex/complex.h>
#include <optional>
#include "mesh.hpp"
#ifdef POLYSCOPE_DEFINED
#include <polyscope/surface_mesh.h>
@ -22,36 +22,38 @@ class VCGTriMeshVertex : public vcg::Vertex<VCGTriMeshUsedTypes,
vcg::vertex::Color4b,
vcg::vertex::VFAdj,
vcg::vertex::Qualityd,
vcg::vertex::VEAdj>
{};
class VCGTriMeshFace
: public vcg::Face<VCGTriMeshUsedTypes, vcg::face::FFAdj, vcg::face::VFAdj,
vcg::face::VertexRef, vcg::face::BitFlags,
vcg::vertex::VEAdj> {};
class VCGTriMeshFace : public vcg::Face<VCGTriMeshUsedTypes,
vcg::face::FFAdj,
vcg::face::VFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::Normal3d> {};
class VCGTriMeshEdge : public vcg::Edge<VCGTriMeshUsedTypes, vcg::edge::VertexRef, vcg::edge::VEAdj>
{};
class VCGTriMeshEdge : public vcg::Edge<VCGTriMeshUsedTypes,
vcg::edge::VertexRef,
vcg::edge::VEAdj> {};
class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
std::vector<VCGTriMeshFace>,
std::vector<VCGTriMeshEdge>>,
public Mesh
{
public:
public Mesh {
public:
VCGTriMesh();
VCGTriMesh(const std::string &filename);
bool load(const std::filesystem::path &meshFilePath) override;
bool load(std::istringstream &offInputStream);
VCGTriMesh(const std::string& filename);
bool load(const std::filesystem::path& meshFilePath) override;
bool load(std::istringstream& offInputStream);
Eigen::MatrixX3d getVertices() const;
Eigen::MatrixX3i getFaces() const;
bool save(const std::string plyFilename);
template <typename MeshElement> size_t getIndex(const MeshElement &element) {
template <typename MeshElement>
size_t getIndex(const MeshElement& element) {
return vcg::tri::Index<VCGTriMesh>(*this, element);
}
#ifdef POLYSCOPE_DEFINED
polyscope::SurfaceMesh *registerForDrawing(
const std::optional<std::array<double, 3>> &desiredColor = std::nullopt,
const bool &shouldEnable = true) const;
polyscope::SurfaceMesh* registerForDrawing(
const std::optional<std::array<float, 3>>& desiredColor = std::nullopt,
const bool& shouldEnable = true) const;
#endif
Eigen::MatrixX2i getEdges() const;
void updateEigenEdgeAndVertices();