Refactoring

This commit is contained in:
iasonmanolas 2022-01-14 15:02:27 +02:00
parent 282e1609c8
commit 4650af984d
15 changed files with 708 additions and 279 deletions

View File

@ -889,25 +889,26 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const Node &edgeNode = pMesh->nodes[ev];
internalForcesContributionFromThisEdge[evi].first = edgeNode.vi;
const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
// refElemOtherVertex can be precomputed
const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
if (edgeNode.referenceElement != &e) {
internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
}
const size_t vi = edgeNode.vi;
// #pragma omp parallel for schedule(static) num_threads(6)
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
&& fixedVertices.at(edgeNode.vi).contains(dofi);
if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi};
// Axial force computation
const double e_k = element.length - element.initialLength;
const double e_kDeriv = computeDerivativeElementLength(e, dui);
const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
const VertexPointer &rev_j = edgeNode.referenceElement->cV(0);
const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1);
// refElemOtherVertex can be precomputed
const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j;
const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex];
if (edgeNode.referenceElement != &e) {
internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi;
}
const size_t vi = edgeNode.vi;
// #pragma omp parallel for schedule(static) num_threads(6)
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi]
&& fixedVertices.at(edgeNode.vi)
.contains(dofi);
if (!isDofConstrainedFor_ev) {
DifferentiateWithRespectTo dui{ev, dofi};
// Axial force computation
const double e_k = element.length - element.initialLength;
const double e_kDeriv = computeDerivativeElementLength(e, dui);
const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(axialForce_dofi)) {
@ -984,7 +985,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
internalForcesContributionFromThisEdge[evi].second[dofi]
= axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
+ torsionalForce_dofi;
}
}
if (edgeNode.referenceElement != &e) {
const bool isDofConstrainedFor_refElemOtherVertex
= isVertexConstrained[refElemOtherVertexNode.vi]
@ -1072,7 +1073,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
}
#endif
}
pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN();
// pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN();
@ -2736,54 +2736,52 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
//#endif
void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const
void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePath) const
{
std::filesystem::create_directories(folderPath);
nlohmann::json json;
json[jsonLabels.gamma] = gamma;
json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false";
json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false";
json[jsonLabels.beVerbose] = beVerbose ? "true" : "false";
json[jsonLabels.Dtini] = Dtini;
json[jsonLabels.xi] = xi;
json[GET_VARIABLE_NAME(shouldDraw)] = shouldDraw;
json[GET_VARIABLE_NAME(beVerbose)] = beVerbose;
json[GET_VARIABLE_NAME(shouldCreatePlots)] = shouldCreatePlots;
json[GET_VARIABLE_NAME(Dtini)] = Dtini;
json[GET_VARIABLE_NAME(xi)] = xi;
json[GET_VARIABLE_NAME(gamma)] = gamma;
json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold;
if (maxDRMIterations.has_value()) {
json[jsonLabels.maxDRMIterations] = maxDRMIterations.value();
json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
}
if (debugModeStep.has_value()) {
json[jsonLabels.debugModeStep] = debugModeStep.value();
json[GET_VARIABLE_NAME(debugModeStep)] = debugModeStep.value();
}
if (totalTranslationalKineticEnergyThreshold.has_value()) {
json[jsonLabels.totalTranslationalKineticEnergyThreshold]
json[GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold)]
= totalTranslationalKineticEnergyThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[jsonLabels.averageResidualForcesCriterionThreshold]
= averageResidualForcesCriterionThreshold.value();
}
if (averageResidualForcesCriterionThreshold.has_value()) {
json[jsonLabels.averageResidualForcesCriterionThreshold]
json[GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold)]
= averageResidualForcesCriterionThreshold.value();
}
if (linearGuessForceScaleFactor.has_value()) {
json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value();
json[GET_VARIABLE_NAME(linearGuessForceScaleFactor)] = linearGuessForceScaleFactor.value();
}
if (viscousDampingFactor.has_value()) {
json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
}
json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping;
// if (intermediateResultsSaveStep.has_value()) {
// json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
// }
if (saveIntermediateBestStates.has_value()) {
json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
? "true"
: "false";
}
// if (saveIntermediateBestStates.has_value()) {
// json[GET_VARIABLE_NAME(saveIntermediateBestStates)] = saveIntermediateBestStates.value()
// ? "true"
// : "false";
// }
const std::string jsonFilename = "drmSettings.json";
std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
std::filesystem::create_directories(jsonFilePath.parent_path());
std::ofstream jsonFile(jsonFilePath);
std::cout << "Saving DRM settings to:" << jsonFilePath << std::endl;
jsonFile << json;
jsonFile.close();
}
bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath)
@ -2807,41 +2805,40 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
std::ifstream ifs(jsonFilePath.string());
ifs >> json;
if (json.contains(jsonLabels.shouldDraw)) {
shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false;
if (json.contains(GET_VARIABLE_NAME(shouldDraw))) {
shouldDraw = json.at(GET_VARIABLE_NAME(shouldDraw)) == "true" ? true : false;
}
if (json.contains(jsonLabels.beVerbose)) {
beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false;
if (json.contains(GET_VARIABLE_NAME(beVerbose))) {
beVerbose = json.at(GET_VARIABLE_NAME(beVerbose)) == "true" ? true : false;
}
if (json.contains(jsonLabels.shouldCreatePlots)) {
shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false;
if (json.contains(GET_VARIABLE_NAME(shouldCreatePlots))) {
shouldCreatePlots = json.at(GET_VARIABLE_NAME(shouldCreatePlots)) == "true" ? true : false;
}
if (json.contains(jsonLabels.Dtini)) {
Dtini = json.at(jsonLabels.Dtini);
if (json.contains(GET_VARIABLE_NAME(Dtini))) {
Dtini = json.at(GET_VARIABLE_NAME(Dtini));
}
if (json.contains(jsonLabels.xi)) {
xi = json.at(jsonLabels.xi);
if (json.contains(GET_VARIABLE_NAME(xi))) {
xi = json.at(GET_VARIABLE_NAME(xi));
}
if (json.contains(jsonLabels.maxDRMIterations)) {
maxDRMIterations = json.at(jsonLabels.maxDRMIterations);
if (json.contains(GET_VARIABLE_NAME(maxDRMIterations))) {
maxDRMIterations = json.at(GET_VARIABLE_NAME(maxDRMIterations));
}
if (json.contains(jsonLabels.debugModeStep)) {
debugModeStep = json.at(jsonLabels.debugModeStep);
if (json.contains(GET_VARIABLE_NAME(debugModeStep))) {
debugModeStep = json.at(GET_VARIABLE_NAME(debugModeStep));
}
if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) {
if (json.contains(GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold))) {
totalTranslationalKineticEnergyThreshold = json.at(
jsonLabels.totalTranslationalKineticEnergyThreshold);
GET_VARIABLE_NAME(totalTranslationalKineticEnergyThreshold));
}
if (json.contains(jsonLabels.gamma)) {
gamma = json.at(jsonLabels.gamma);
if (json.contains(GET_VARIABLE_NAME(gamma))) {
gamma = json.at(GET_VARIABLE_NAME(gamma));
}
if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) {
if (json.contains(GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold))) {
averageResidualForcesCriterionThreshold = json.at(
jsonLabels.averageResidualForcesCriterionThreshold);
GET_VARIABLE_NAME(averageResidualForcesCriterionThreshold));
}
if (json.contains(jsonLabels.linearGuessForceScaleFactor)) {
linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor);
if (json.contains(GET_VARIABLE_NAME(linearGuessForceScaleFactor))) {
linearGuessForceScaleFactor = json.at(GET_VARIABLE_NAME(linearGuessForceScaleFactor));
}
// if (json.contains(GET_VARIABLE_NAME(intermediateResultsSaveStep))) {
@ -2854,6 +2851,7 @@ bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePat
? true
: false;
}
// std::cout << json.dump() << std::endl;
return true;
}

View File

@ -33,7 +33,6 @@ class DRMSimulationModel : public SimulationModel
public:
struct Settings
{
// bool isDebugMode{false};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
@ -46,8 +45,8 @@ public:
int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
double totalResidualForcesNormThreshold{1e-8};
double totalExternalForcesNormPercentageTermination{1e-5};
std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep;
std::optional<double> totalTranslationalKineticEnergyThreshold;
@ -55,11 +54,11 @@ public:
std::optional<double> linearGuessForceScaleFactor;
// std::optional<int> intermediateResultsSaveStep;
std::optional<bool> saveIntermediateBestStates;
Settings() {}
void save(const std::filesystem::path &folderPath) const;
bool load(const std::filesystem::path &filePath);
std::optional<double> viscousDampingFactor;
bool useKineticDamping{true};
Settings() {}
void save(const std::filesystem::path &jsonFilePath) const;
bool load(const std::filesystem::path &filePath);
struct JsonLabels
{
const std::string shouldDraw{"shouldDraw"};
@ -67,14 +66,16 @@ public:
const std::string shouldCreatePlots{"shouldCreatePlots"};
const std::string Dtini{"DtIni"};
const std::string xi{"xi"};
const std::string gamma{"gamma"};
const std::string totalResidualForcesNormThreshold;
const std::string maxDRMIterations{"maxIterations"};
const std::string debugModeStep{"debugModeStep"};
const std::string totalTranslationalKineticEnergyThreshold{
"totalTranslationaKineticEnergyThreshold"};
const std::string gamma{"gamma"};
const std::string averageResidualForcesCriterionThreshold{
"averageResidualForcesThreshold"};
const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"};
const std::string viscousDampingFactor{"viscousDampingFactor"};
};
static JsonLabels jsonLabels;
};

View File

@ -25,13 +25,13 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const {
return eigenEdgeNormals;
}
bool VCGEdgeMesh::save(const std::string &plyFilename)
bool VCGEdgeMesh::save(const std::filesystem::__cxx11::path &meshFilePath)
{
std::string filename = plyFilename;
std::string filename = meshFilePath;
if (filename.empty()) {
filename = std::filesystem::current_path().append(getLabel() + ".ply").string();
} else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) {
filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply").string();
} else if (std::filesystem::is_directory(std::filesystem::path(meshFilePath))) {
filename = std::filesystem::path(meshFilePath).append(getLabel() + ".ply").string();
}
assert(std::filesystem::path(filename).extension().string() == ".ply");
unsigned int mask = 0;
@ -197,11 +197,11 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
return true;
}
bool VCGEdgeMesh::load(const std::string &plyFilename)
bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath)
{
std::string usedPath = plyFilename;
if (std::filesystem::path(plyFilename).is_relative()) {
usedPath = std::filesystem::absolute(plyFilename).string();
std::string usedPath = meshFilePath;
if (std::filesystem::path(meshFilePath).is_relative()) {
usedPath = std::filesystem::absolute(meshFilePath).string();
}
assert(std::filesystem::exists(usedPath));
Clear();
@ -217,15 +217,15 @@ bool VCGEdgeMesh::load(const std::string &plyFilename)
getEdges(eigenEdges);
getVertices(eigenVertices);
vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
label = std::filesystem::path(plyFilename).stem().string();
label = std::filesystem::path(meshFilePath).stem().string();
const bool printDebugInfo = false;
if (printDebugInfo) {
std::cout << plyFilename << " was loaded successfuly." << std::endl;
std::cout << meshFilePath << " was loaded successfuly." << std::endl;
std::cout << "Mesh has " << EN() << " edges." << std::endl;
}
label = std::filesystem::path(plyFilename).stem().string();
label = std::filesystem::path(meshFilePath).stem().string();
return true;
}

View File

@ -61,8 +61,9 @@ public:
void set(const std::vector<double> &vertexPositions, const std::vector<int> &edges);
void removeDuplicateVertices();
void deleteDanglingVertices();
void deleteDanglingVertices(vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
virtual void deleteDanglingVertices();
virtual void deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu);
void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const;
@ -72,8 +73,8 @@ public:
// bool loadUsingNanoply(const std::string &plyFilename);
bool load(const std::string &plyFilename) override;
bool save(const std::string &plyFilename = std::string()) override;
bool load(const std::filesystem::path &meshFilePath) override;
bool save(const std::filesystem::path &meshFilePath = std::filesystem::path()) override;
bool createSpanGrid(const size_t squareGridDimension);
bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight);

View File

@ -11,7 +11,7 @@
namespace PolygonalRemeshing {
//std::shared_ptr<VCGPolyMesh>
std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
inline std::shared_ptr<VCGPolyMesh> remeshWithPolygons(VCGTriMesh &surface)
{
// const std::string tileIntoFilePath(
// "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply");

View File

@ -1,7 +1,7 @@
#ifndef MESH_HPP
#define MESH_HPP
//#include <Eigen/Core>
#include <filesystem>
#include <string>
class Mesh
@ -11,8 +11,8 @@ protected:
public:
virtual ~Mesh() = default;
virtual bool load(const std::string &filePath) { return false; }
virtual bool save(const std::string &filePath) { return false; }
virtual bool load(const std::filesystem::path &meshFilePath) { return false; }
virtual bool save(const std::filesystem::path &meshFilePath) { return false; }
std::string getLabel() const;
void setLabel(const std::string &newLabel);

View File

@ -5,6 +5,7 @@
#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>
@ -64,14 +65,25 @@ class VCGPolyMesh
public Mesh
{
public:
// virtual bool load(const std::string &filename) override {
// int mask;
// vcg::tri::io::Importer<VCGPolyMesh>::LoadMask(filename.c_str(), mask);
// int error = vcg::tri::io::Importer<VCGPolyMesh>::Open(
// *this, filename.c_str(), mask);
// if (error != 0) {
// return false;
// }
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;
@ -84,11 +96,6 @@ public:
// // 0) {
// // std::cout << "Could not load tri mesh" << std::endl;
// // }
// vcg::tri::UpdateTopology<VCGPolyMesh>::FaceFace(*this);
// vcg::tri::UpdateTopology<VCGPolyMesh>::VertexEdge(*this);
// // vcg::tri::UpdateTopology<VCGPolyMesh>::VertexFace(*this);
// vcg::tri::UpdateNormal<VCGPolyMesh>::PerVertexNormalized(*this);
// return true;
// }
Eigen::MatrixX3d getVertices() const
{
@ -113,6 +120,34 @@ public:
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())
{
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";

317
reducedmodelevaluator.cpp Normal file
View File

@ -0,0 +1,317 @@
#include "reducedmodelevaluator.hpp"
#include "hexagonremesher.hpp"
#include "trianglepatterngeometry.hpp"
#include <filesystem>
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
ReducedModelEvaluator::ReducedModelEvaluator()
{
}
std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResults)
{
// std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
// std::filesystem::path tileIntoSurfaceFilePath{
// "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
// assert(std::filesystem::exists(tileIntoSurfaceFilePath));
// pTileIntoSurface->load(tileIntoSurfaceFilePath);
//Set required file paths
const std::filesystem::path tileInto_triMesh_filename
= "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
"instantMeshes_plane_34.ply";
const std::filesystem::path reducedPatternFilePath
= "/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
const std::filesystem::path intermediateResultsDirectoryPath
= "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
// const std::filesystem::path drmSettingsFilePath
// = "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
// "defaultDRMSimulationSettings.json";
DRMSimulationModel::Settings drmSimulationSettings;
// drmSimulationSettings.load(drmSettingsFilePath);
drmSimulationSettings.linearGuessForceScaleFactor = 1;
drmSimulationSettings.debugModeStep = 1000;
drmSimulationSettings.beVerbose = true;
constexpr bool shouldRerunFullPatternSimulation = false;
const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
"22Hex_randomBending1",
"22Hex_randomBending2",
"22Hex_randomBending3",
"22Hex_randomBending4",
"22Hex_randomBending5",
"22Hex_randomBending6",
"22Hex_randomBending7",
"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_randomBending20",
"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.3N",
"22Hex_shear_2N",
"22Hex_shear_5N",
"22Hex_axial_10N",
"22Hex_axial_20N",
"22Hex_cylinder_0.05N",
"22Hex_cylinder_0.1N",
"22Hex_s_0.05N",
"22Hex_s_0.1N"};
//Load surface
std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
VCGTriMesh tileInto_triMesh;
const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
assert(surfaceLoadSuccessfull);
return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
}();
const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
(optimizationResults.baseTriangle.cP(1)
+ optimizationResults.baseTriangle.cP(
2))
/ 2);
pTileIntoSurface->moveToCenter();
const double scaleFactor = optimizedBaseTriangleHeight
/ pTileIntoSurface->getAverageFaceRadius();
vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
// tileIntoSurface.registerForDrawing();
// polyscope::show();
//Tile full pattern into surface
std::vector<PatternGeometry> fullPatterns(1);
fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
//// Base triangle pattern might contain dangling vertices.Remove those
fullPatterns[0].interfaceNodeIndex = 3;
fullPatterns[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,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledFullVi,
perPatternIndexTiledFullPatternEdgeIndices);
pTiledFullPattern->setLabel("Tiled_full_patterns");
// pTiledFullPattern->registerForDrawing();
//Tile reduced pattern into surface
PatternGeometry reducedPattern;
reducedPattern.load(reducedPatternFilePath);
reducedPattern.deleteDanglingVertices();
assert(reducedPattern.interfaceNodeIndex == 1);
std::vector<PatternGeometry> reducedPatterns(1);
reducedPatterns[0].copy(reducedPattern);
const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);
std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledReducedVi;
std::shared_ptr<PatternGeometry> pTiledReducedPattern
= PatternGeometry::tilePattern(reducedPatterns,
{0},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledReducedVi,
perPatternIndexTiledReducedPatternEdgeIndices);
pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
// pTiledReducedPattern->registerForDrawing();
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
fullToReducedViMap; //of only the common vertices
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
reducedToFullViMap; //of only the common vertices
for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
}
constructInverseMap(fullToReducedViMap, reducedToFullViMap);
std::vector<size_t> tilledFullPatternInterfaceVi;
tilledFullPatternInterfaceVi.clear();
tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
size_t viIndex = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
}
//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
pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
if (optimizationResults.fullPatternYoungsModulus == 0) {
std::cerr << "Full pattern's young modulus not found." << std::endl;
std::terminate();
}
pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
optimizationResults.fullPatternYoungsModulus);
pTiledFullPattern_simulationMesh->reset();
////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];
ReducedModelOptimization::Results::applyOptimizationResults_elements(
optimizationResults, pTiledReducedPattern_simulationMesh);
pTiledReducedPattern_simulationMesh->reset();
std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
const std::filesystem::path scenariosDirectoryPath
= "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
const std::filesystem::path tiledReducedPatternJobFilePath
= std::filesystem::path(scenariosDirectoryPath)
.append(jobLabel)
.append("SimulationJobs")
.append("Reduced")
.append(SimulationJob::jsonDefaultFileName);
//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);
// pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
// pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
// polyscope::show();
//Save reduced job
const std::filesystem::path tesellatedResultsFolderPath
= std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
const std::filesystem::path surfaceFolderPath = std::filesystem::path(
tesellatedResultsFolderPath)
.append(pTileIntoSurface->getLabel());
const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
.append(scenarioLabel);
const std::filesystem::path reducedJobDirectoryPath
= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
std::filesystem::create_directories(reducedJobDirectoryPath);
pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
//Run scenario
////Full
const std::string patternLabel = [&]() {
const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
const int numberOfOccurences = std::count_if(patternLabel.begin(),
patternLabel.end(),
[](char c) { return c == '#'; });
if (numberOfOccurences == 0) {
return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
+ optimizationResults.baseTriangleFullPattern.getLabel();
} else if (numberOfOccurences == 1) {
return optimizationResults.baseTriangleFullPattern.getLabel();
}
}();
const auto fullResultsFolderPath
= std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
std::filesystem::remove_all(fullResultsFolderPath);
}
const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
scenarioDirectoryPath)
.append(patternLabel)
.append("SimulationJob");
SimulationResults simulationResults_tiledFullPattern_drm;
if (std::filesystem::exists(fullResultsFolderPath)) {
//Load full pattern results
assert(std::filesystem::exists(fullPatternJobFolderPath));
simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
fullPatternJobFolderPath);
simulationResults_tiledFullPattern_drm.converged = true;
} else {
//Full
std::cout << "Executing:" << jobLabel << std::endl;
DRMSimulationModel drmSimulationModel;
simulationResults_tiledFullPattern_drm
= drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
}
std::filesystem::create_directories(fullResultsFolderPath);
simulationResults_tiledFullPattern_drm.save(
std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
if (!simulationResults_tiledFullPattern_drm.converged) {
std::cerr << "Full pattern simulation failed." << std::endl;
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_tiledReducedPattern
= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
// simulationResults_tiledReducedPattern.registerForDrawing();
// simulationResults_tiledFullPattern_drm.registerForDrawing();
// polyscope::show();
//measure distance
const double distance_fullDrmToReduced
= simulationResults_tiledFullPattern_drm
.computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
double distance_fullSumOfAllVerts = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
.displacements[fullToReducedPair.first]
.getTranslation()
.norm();
}
const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
/ distance_fullSumOfAllVerts;
distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
}
//#ifndef POLYSCOPE_DEFINED
// return distances_drm2reduced;
//#else
//report distance
// csvFile csv_results("", flse);
csvFile csv_results({}, false);
// csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
// .append("results.csv")
// .string(),
// false);
csv_results /*<< "Job Label"*/
<< "drm2Reduced"
<< "norm_drm2Reduced";
csv_results << endrow;
double sumOfNormalizedFull2Reduced = 0;
for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
<< distance_normalizedFullDrmToReduced;
csv_results << endrow;
sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
}
std::cout << "Average normalized error per scenario:"
<< sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
return distances_normalizedDrm2reduced;
//#endif
}

14
reducedmodelevaluator.hpp Normal file
View File

@ -0,0 +1,14 @@
#ifndef REDUCEDMODELEVALUATOR_HPP
#define REDUCEDMODELEVALUATOR_HPP
#include "reducedmodeloptimizer_structs.hpp"
class ReducedModelEvaluator
{
public:
ReducedModelEvaluator();
static std::vector<double> evaluateReducedModel(
ReducedModelOptimization::Results &optimizationResults);
};
#endif // REDUCEDMODELEVALUATOR_HPP

View File

@ -7,7 +7,7 @@
#include "unordered_map"
#include <string>
namespace ReducedPatternOptimization {
namespace ReducedModelOptimization {
enum BaseSimulationScenario {
Axial,
@ -107,21 +107,29 @@ inline int getParameterIndex(const std::string &s)
struct Settings
{
inline static std::string defaultFilename{"OptimizationSettings.json"};
std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy;
std::vector<std::vector<OptimizationParameterIndex>> optimizationStrategy = {
{E, A, I2, I3, J, R, Theta}};
// {A, I2, I3, J, R, Theta}};
std::optional<std::vector<double>>
optimizationVariablesGroupsWeights; //TODO:should be removed in the future if not a splitting strategy is used for the optimization
enum NormalizationStrategy { NonNormalized, Epsilon };
inline static std::vector<std::string> normalizationStrategyStrings{"NonNormalized", "Epsilon"};
NormalizationStrategy normalizationStrategy{Epsilon};
std::array<xRange, NumberOfOptimizationVariables> variablesRanges;
std::array<xRange, NumberOfOptimizationVariables> variablesRanges{xRange{"E", 0.001, 1000},
xRange{"A", 0.001, 1000},
xRange{"I2", 0.001, 1000},
xRange{"I3", 0.001, 1000},
xRange{"J", 0.001, 1000},
xRange{"R", 0.05, 0.95},
xRange{"Theta", -30, 30}};
int numberOfFunctionCalls{100000};
double solverAccuracy{1e-3};
double translationNormalizationEpsilon{1e-4};
double rotationNormalizationEpsilon{vcg::math::ToRad(3.0)};
double translationEpsilon{4e-3};
double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
struct ObjectiveWeights
{
double translational{1};
double rotational{1};
double translational{1.2};
double rotational{0.8};
bool operator==(const ObjectiveWeights &other) const
{
return this->translational == other.translational
@ -155,15 +163,15 @@ struct Settings
void setDefault()
{
ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
numberOfFunctionCalls = 100000;
// ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
// ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
// ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
// ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
// ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
// ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
// ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
// variablesRanges = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
// numberOfFunctionCalls = 100000;
enum OptimizationParameterComparisonScenarioIndex {
AllVar,
@ -179,12 +187,12 @@ struct Settings
NumberOfScenarios
};
const std::vector<
std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
optimizationParameters = [&]() {
std::vector<
std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>>
std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>>
optimizationParameters(NumberOfScenarios);
using namespace ReducedPatternOptimization;
using namespace ReducedModelOptimization;
optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
optimizationParameters[GeoYM] = {{R, Theta, E}};
optimizationParameters[MatGeo] = {{A, I2, I3, J, R, Theta}};
@ -198,43 +206,28 @@ struct Settings
return optimizationParameters;
}();
constexpr OptimizationParameterComparisonScenarioIndex scenario = MatGeo;
constexpr OptimizationParameterComparisonScenarioIndex scenario = AllVar;
optimizationStrategy = optimizationParameters[scenario];
if (scenario == YMGeo_Mat) {
optimizationVariablesGroupsWeights = {0.15, 0.85};
}
normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
//#ifdef POLYSCOPE_DEFINED
//#else
translationNormalizationEpsilon = 4e-3;
rotationNormalizationEpsilon = vcg::math::ToRad(3.0);
// normalizationStrategy = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
// translationNormalizationEpsilon = 0;
// rotationNormalizationEpsilon = vcg::math::ToRad(0.0);
// solverAccuracy = 1e-2;
// translationNormalizationParameter = 0;
// rotationNormalizationParameter = 0;
// translationNormalizationParameter = 1e-3;
// rotationNormalizationParameter = vcg::math::ToRad(3.0);
// translationNormalizationParameter = 5e-2;
// rotationNormalizationParameter = vcg::math::ToRad(3.0);
//#endif
solverAccuracy = 1e-2;
//TODO refactor that
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.5;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].rotational
= 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.5;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].rotational
= 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.2;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].rotational
= 2
- perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 1.95;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].rotational
= 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 1.2;
perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].rotational
= 2 - perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].translational = 1.174011;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].translational = 1.690460;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending].translational = 1.812089;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].translational = 0.127392;
// perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle].translational = 0.503043;
// for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
// baseScenario++) {
// perBaseScenarioObjectiveWeights[baseScenario].rotational
// = 2 - perBaseScenarioObjectiveWeights[baseScenario].translational;
// }
}
void save(const std::filesystem::path &saveToPath)
@ -242,10 +235,10 @@ struct Settings
assert(std::filesystem::is_directory(saveToPath));
nlohmann::json json;
json[JsonKeys::OptimizationStrategy] = optimizationStrategy;
json[GET_VARIABLE_NAME(optimizationStrategy)] = optimizationStrategy;
if (optimizationVariablesGroupsWeights.has_value()) {
json[JsonKeys::OptimizationStrategyGroupWeights] = optimizationVariablesGroupsWeights
.value();
json[GET_VARIABLE_NAME(ptimizationStrategyGroupWeights)]
= optimizationVariablesGroupsWeights.value();
}
//write x ranges
const std::array<std::tuple<std::string, double, double>, NumberOfOptimizationVariables>
@ -267,8 +260,8 @@ struct Settings
// = xRange.toString();
// }
json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls;
json[JsonKeys::SolverAccuracy] = solverAccuracy;
json[GET_VARIABLE_NAME(numberOfFunctionCalls)] = numberOfFunctionCalls;
json[GET_VARIABLE_NAME(solverAccuracy)] = solverAccuracy;
//Objective weights
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios> objectiveWeightsPairs;
std::transform(perBaseScenarioObjectiveWeights.begin(),
@ -279,11 +272,15 @@ struct Settings
objectiveWeights.rotational);
});
json[JsonKeys::ObjectiveWeights] = objectiveWeightsPairs;
json[GET_VARIABLE_NAME(translationEpsilon)] = translationEpsilon;
json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(
angularDistanceEpsilon);
std::filesystem::path jsonFilePath(
std::filesystem::path(saveToPath).append(defaultFilename));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
jsonFile.close();
}
bool load(const std::filesystem::path &jsonFilepath)
@ -298,14 +295,14 @@ struct Settings
std::ifstream ifs(jsonFilepath);
ifs >> json;
if (json.contains(JsonKeys::OptimizationStrategy)) {
if (json.contains(GET_VARIABLE_NAME(optimizationStrategy))) {
optimizationStrategy
= std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>(
(json.at(JsonKeys::OptimizationStrategy)));
= std::vector<std::vector<ReducedModelOptimization::OptimizationParameterIndex>>(
(json.at(GET_VARIABLE_NAME(optimizationStrategy))));
}
if (json.contains(JsonKeys::OptimizationStrategyGroupWeights)) {
if (json.contains(GET_VARIABLE_NAME(optimizationStrategyGroupWeights))) {
optimizationVariablesGroupsWeights = std::vector<double>(
json[JsonKeys::OptimizationStrategyGroupWeights]);
json[GET_VARIABLE_NAME(optimizationStrategyGroupWeights)]);
}
//read x ranges
if (json.contains(JsonKeys::OptimizationVariables)) {
@ -328,8 +325,12 @@ struct Settings
}
}
numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
solverAccuracy = json.at(JsonKeys::SolverAccuracy);
if (json.contains(GET_VARIABLE_NAME(numberOfFunctionCalls))) {
numberOfFunctionCalls = json.at(GET_VARIABLE_NAME(numberOfFunctionCalls));
}
if (json.contains(GET_VARIABLE_NAME(solverAccuracy))) {
solverAccuracy = json.at(GET_VARIABLE_NAME(solverAccuracy));
}
//Objective weights
if (json.contains(JsonKeys::ObjectiveWeights)) {
std::array<std::pair<double, double>, NumberOfBaseSimulationScenarios>
@ -342,6 +343,15 @@ struct Settings
objectiveWeightsPair.second};
});
}
if (json.contains(GET_VARIABLE_NAME(translationalNormalizationEpsilon))) {
translationEpsilon
= json[GET_VARIABLE_NAME(translationalNormalizationEpsilon)];
}
if (json.contains(GET_VARIABLE_NAME(angularDistanceEpsilon))) {
angularDistanceEpsilon = vcg::math::ToRad(
static_cast<double>(json[GET_VARIABLE_NAME(angularDistanceEpsilon)]));
}
// perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
@ -381,7 +391,8 @@ struct Settings
}
os << "Function Calls";
os << "Solution Accuracy";
os << "Normalization strategy";
os << "Normalization trans epsilon";
os << "Normalization rot epsilon(deg)";
os << JsonKeys::ObjectiveWeights;
os << "Optimization parameters";
// os << std::endl;
@ -397,8 +408,8 @@ struct Settings
}
os << numberOfFunctionCalls;
os << solverAccuracy;
os << normalizationStrategyStrings[normalizationStrategy] + "_"
+ std::to_string(translationNormalizationEpsilon);
os << std::to_string(translationEpsilon);
os << std::to_string(vcg::math::ToDeg(angularDistanceEpsilon));
std::string objectiveWeightsString;
objectiveWeightsString += "{";
for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios;
@ -436,8 +447,8 @@ struct Settings
return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls
&& settings1.variablesRanges == settings2.variablesRanges
&& settings1.solverAccuracy == settings2.solverAccuracy && haveTheSameObjectiveWeights
&& settings1.translationNormalizationEpsilon
== settings2.translationNormalizationEpsilon;
&& settings1.translationEpsilon
== settings2.translationEpsilon;
}
inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
@ -504,7 +515,7 @@ struct Settings
std::vector<std::pair<std::string, double>> optimalXNameValuePairs;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
double fullPatternYoungsModulus;
double fullPatternYoungsModulus{0};
PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
@ -646,7 +657,7 @@ struct Settings
settings.save(saveToPath);
#ifdef POLYSCOPE_DEFINED12
#ifdef POLYSCOPE_DEFINED
std::cout << "Saved optimization results to:" << saveToPath << std::endl;
#endif
}
@ -748,14 +759,16 @@ struct Settings
template<typename MeshType>
static void applyOptimizationResults_innerHexagon(
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
{
std::unordered_map<std::string, double>
optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
assert(optimalXVariables.contains("R") && optimalXVariables.contains("Theta"));
assert(
(optimalXVariables.contains("R") && optimalXVariables.contains("Theta"))
|| (optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")));
if (optimalXVariables.contains("HexSize")) {
applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"),
optimalXVariables.at("HexAngle"),
@ -807,8 +820,8 @@ struct Settings
}
static void applyOptimizationResults_elements(
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const std::shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const std::shared_ptr<SimulationMesh> &pReducedPattern_simulationMesh)
{
assert(CrossSectionType::name == RectangularBeamDimensions::name);
// Set reduced parameters fron the optimization results
@ -817,55 +830,37 @@ struct Settings
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
const std::string ALabel = "A";
if (optimalXVariables.contains(ALabel)) {
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
}
}
assert(optimalXVariables.contains(ALabel));
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";
if (optimalXVariables.contains(ymLabel)) {
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setMaterial(elementMaterial);
}
}
assert(optimalXVariables.contains(ymLabel));
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
const std::string JLabel = "J";
if (optimalXVariables.contains(JLabel)) {
const double J = optimalXVariables.at(JLabel);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.J = J;
}
}
assert(optimalXVariables.contains(JLabel));
const double J = optimalXVariables.at(JLabel);
const std::string I2Label = "I2";
if (optimalXVariables.contains(I2Label)) {
const double I2 = optimalXVariables.at(I2Label);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.I2 = I2;
}
}
assert(optimalXVariables.contains(I2Label));
const double I2 = optimalXVariables.at(I2Label);
const std::string I3Label = "I3";
if (optimalXVariables.contains(I3Label)) {
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.dimensions.inertia.I3 = I3;
}
assert(optimalXVariables.contains(I3Label));
const double I3 = optimalXVariables.at(I3Label);
for (int ei = 0; ei < pReducedPattern_simulationMesh->EN(); ei++) {
Element &e = pReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
e.setMaterial(elementMaterial);
e.dimensions.inertia.J = J;
e.dimensions.inertia.I2 = I2;
e.dimensions.inertia.I3 = I3;
}
pTiledReducedPattern_simulationMesh->reset();
}
#if POLYSCOPE_DEFINED

View File

@ -113,25 +113,27 @@ class SimulationJob {
} jsonLabels;
public:
std::shared_ptr<SimulationMesh> pMesh;
std::string label{"empty_job"};
std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
SimulationJob(
const std::shared_ptr<SimulationMesh> &m, const std::string &label,
const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
const std::unordered_map<VertexIndex, Vector6d> &ef = {},
const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
: pMesh(m), label(label), constrainedVertices(cv),
nodalExternalForces(ef), nodalForcedDisplacements(fd) {}
SimulationJob() {}
SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
inline static std::string jsonDefaultFileName = "SimulationJob.json";
std::shared_ptr<SimulationMesh> pMesh;
std::string label{"empty_job"};
std::unordered_map<VertexIndex, std::unordered_set<int>> constrainedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
SimulationJob(const std::shared_ptr<SimulationMesh> &m,
const std::string &label,
const std::unordered_map<VertexIndex, std::unordered_set<int>> &cv,
const std::unordered_map<VertexIndex, Vector6d> &ef = {},
const std::unordered_map<VertexIndex, Eigen::Vector3d> &fd = {})
: pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef),
nodalForcedDisplacements(fd)
{}
SimulationJob() {}
SimulationJob(const std::string &jsonFilename) { load(jsonFilename); }
bool isEmpty()
{
return constrainedVertices.empty() && nodalExternalForces.empty()
&& nodalForcedDisplacements.empty() && pMesh == nullptr;
bool isEmpty()
{
return constrainedVertices.empty() && nodalExternalForces.empty()
&& nodalForcedDisplacements.empty() && pMesh == nullptr;
}
void remap(const std::unordered_map<size_t, size_t> &sourceToDestinationViMap,
SimulationJob &destination_simulationJob) const
@ -489,6 +491,7 @@ struct SimulationResults
/*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments
* vector contains the same info using euler angles
*/
inline const static std::string defaultJsonFilename{"SimulationResults.json"};
bool converged{false};
std::shared_ptr<SimulationJob> pJob;
SimulationHistory history;
@ -553,6 +556,16 @@ struct SimulationResults
std::filesystem::create_directories(resultsFolderPath);
Eigen::write_binary(std::filesystem::path(resultsFolderPath).append(filename).string(), m);
nlohmann::json json;
json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
std::filesystem::path jsonFilePath(
std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
jsonFile.close();
saveDeformedModel(resultsFolderPath.string());
}
@ -578,6 +591,28 @@ struct SimulationResults
this->pJob = pJob;
load(loadFromPath);
}
static double computeDistance(
const SimulationResults &resultsA,
const SimulationResults &resultsB,
const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap)
{
double distance = 0;
for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
const double vertexToVertexDistance
= (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
- resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
.norm();
distance += vertexToVertexDistance;
}
return distance;
}
double computeDistance(const SimulationResults &other,
const std::unordered_map<VertexIndex, VertexIndex> &thisToOtherViMap) const
{
return computeDistance(*this, other, thisToOtherViMap);
}
#ifdef POLYSCOPE_DEFINED
void unregister() const
{
@ -732,7 +767,7 @@ struct SimulationResults
}
#endif
private:
void load(const std::filesystem::path &loadFromPath)
bool load(const std::filesystem::path &loadFromPath)
{
converged = true; //assuming it has converged
assert(pJob != nullptr);
@ -756,6 +791,19 @@ private:
* Eigen::AngleAxisd(displacements[vi][5],
Eigen::Vector3d::UnitZ());
}
const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
.append(defaultJsonFilename);
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Simulation results could not be loaded because filepath does "
"not exist:"
<< jsonFilepath << std::endl;
return false;
}
std::ifstream ifs(jsonFilepath);
nlohmann::json json;
ifs >> json;
internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
}
};

View File

@ -86,39 +86,39 @@ std::vector<vcg::Point3d> PatternGeometry::getVertices() const {
return verts;
}
PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern) {
PatternGeometry PatternGeometry::createTile(PatternGeometry &pattern)
{
const size_t fanSize = PatternGeometry().getFanSize();
PatternGeometry fan(createFan(pattern));
PatternGeometry tile(fan);
const size_t fanSize = PatternGeometry().getFanSize();
PatternGeometry fan(createFan(pattern));
PatternGeometry tile(fan);
if (fanSize % 2 == 1) {
vcg::Matrix44d R;
auto rotationAxis = vcg::Point3d(0, 0, 1);
R.SetRotateDeg(180, rotationAxis);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
}
vcg::Matrix44d T;
const double centerAngle = 2 * M_PI / fanSize;
const double triangleHeight = std::sin((M_PI - centerAngle) / 2)
* PatternGeometry().triangleEdgeSize;
T.SetTranslate(0, -2 * triangleHeight, 0);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
if (fanSize % 2 == 1) {
vcg::Matrix44d R;
auto rotationAxis = vcg::Point3d(0, 0, 1);
R.SetRotateDeg(180, rotationAxis);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, R);
}
vcg::Matrix44d T;
const double centerAngle = 2 * M_PI / fanSize;
const double triangleHeight =
std::sin((M_PI - centerAngle) / 2) * PatternGeometry().triangleEdgeSize;
T.SetTranslate(0, -2 * triangleHeight, 0);
vcg::tri::UpdatePosition<PatternGeometry>::Matrix(fan, T);
PatternGeometry fanOfFan = createFan(fan);
vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
PatternGeometry fanOfFan = createFan(fan);
vcg::tri::Append<PatternGeometry, PatternGeometry>::Mesh(tile, fanOfFan);
vcg::tri::Clean<PatternGeometry>::MergeCloseVertex(tile, 0.0000005);
vcg::tri::Allocator<PatternGeometry>::CompactEveryVector(tile);
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(tile);
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(tile);
// for (size_t vi = 0; vi < pattern.vn; vi++) {
// tile.vert[vi].C() = vcg::Color4b::Blue;
// }
// for (size_t vi = 0; vi < pattern.vn; vi++) {
// tile.vert[vi].C() = vcg::Color4b::Blue;
// }
tile.setLabel("tilled_" + pattern.getLabel());
tile.updateEigenEdgeAndVertices();
return tile;
tile.setLabel("tilled_" + pattern.getLabel());
tile.updateEigenEdgeAndVertices();
return tile;
}
PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) {
@ -145,9 +145,9 @@ PatternGeometry::PatternGeometry(PatternGeometry &other) {
vcg::tri::UpdateTopology<PatternGeometry>::EdgeEdge(*this);
}
bool PatternGeometry::load(const std::string &plyFilename)
bool PatternGeometry::load(const std::filesystem::__cxx11::path &meshFilePath)
{
if (!VCGEdgeMesh::load(plyFilename)) {
if (!VCGEdgeMesh::load(meshFilePath)) {
return false;
}
vcg::tri::UpdateTopology<PatternGeometry>::VertexEdge(*this);
@ -777,6 +777,22 @@ double PatternGeometry::computeBaseTriangleHeight() const
return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP());
}
void PatternGeometry::deleteDanglingVertices()
{
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> pu;
VCGEdgeMesh::deleteDanglingVertices(pu);
if (!pu.remap.empty()) {
interfaceNodeIndex
= pu.remap[interfaceNodeIndex]; //TODO:Could this be automatically be determined?
}
}
void PatternGeometry::deleteDanglingVertices(
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VertexPointer> &pu)
{
VCGEdgeMesh::deleteDanglingVertices(pu);
}
vcg::Triangle3<double> PatternGeometry::getBaseTriangle() const
{
return baseTriangle;
@ -814,6 +830,7 @@ PatternGeometry::PatternGeometry(
updateEigenEdgeAndVertices();
}
//TODO: refactor such that it takes as an input a single pattern and tiles into the desired faces without requiring the each face will contain a single pattern
std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
std::vector<ConstPatternGeometry> &patterns,
const std::vector<int> &connectToNeighborsVi,

View File

@ -46,7 +46,7 @@ private:
* vcglib interface.
* */
PatternGeometry(PatternGeometry &other);
bool load(const std::string &plyFilename) override;
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);
@ -124,6 +124,9 @@ private:
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;
};
#endif // FLATPATTERNGEOMETRY_HPP

View File

@ -7,8 +7,8 @@
#include <wrap/io_trimesh/import.h>
//#include <wrap/nanoply/include/nanoplyWrapper.hpp>
bool VCGTriMesh::load(const std::string &filename) {
assert(std::filesystem::exists(filename));
bool VCGTriMesh::load(const std::filesystem::__cxx11::path &meshFilePath) {
assert(std::filesystem::exists(meshFilePath));
unsigned int mask = 0;
mask |= vcg::tri::io::Mask::IOM_VERTCOORD;
mask |= vcg::tri::io::Mask::IOM_VERTNORMAL;
@ -20,7 +20,7 @@ bool VCGTriMesh::load(const std::string &filename) {
// std::filesystem::absolute(filename).string().c_str(), *this, mask)
// != 0) {
if (vcg::tri::io::Importer<VCGTriMesh>::Open(*this,
std::filesystem::absolute(filename).string().c_str())
std::filesystem::absolute(meshFilePath).string().c_str())
!= 0) {
std::cout << "Could not load tri mesh" << std::endl;
return false;
@ -32,7 +32,7 @@ bool VCGTriMesh::load(const std::string &filename) {
vcg::tri::UpdateTopology<VCGTriMesh>::VertexEdge(*this);
vcg::tri::UpdateNormal<VCGTriMesh>::PerVertexNormalized(*this);
label = std::filesystem::path(filename).stem().string();
label = std::filesystem::path(meshFilePath).stem().string();
return true;
}

View File

@ -39,7 +39,7 @@ class VCGTriMesh : public vcg::tri::TriMesh<std::vector<VCGTriMeshVertex>,
public:
VCGTriMesh();
VCGTriMesh(const std::string &filename);
bool load(const std::string &filename) override;
bool load(const std::filesystem::path &meshFilePath) override;
Eigen::MatrixX3d getVertices() const;
Eigen::MatrixX3i getFaces() const;
bool save(const std::string plyFilename);