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

refactoring. Saving of drm settings. TBB instead of openmp

This commit is contained in:
iasonmanolas 2021-07-21 18:00:36 +03:00
parent af4449a9dd
commit 3c879db90f
7 changed files with 3231 additions and 284 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
CMakeLists.txt
View File

@ -43,6 +43,17 @@ download_project(PROJ threed-beam-fea
)
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##TBB
download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR})
message(${TBB_BINARY_DIR})
###Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED)
@ -53,6 +64,7 @@ endif(MSVC)
#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp)
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ltbb")
add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
@ -64,9 +76,8 @@ target_include_directories(${PROJECT_NAME}
)
if(${USE_POLYSCOPE})
find_package(OpenMP REQUIRED)
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA OpenMP::OpenMP_CXX)
target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
else()
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA)
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread)
endif()
target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)

309
drmsimulationmodel.cpp
View File

@ -222,7 +222,7 @@ void DRMSimulationModel::reset()
mCurrentSimulationStep = 0;
history.clear();
constrainedVertices.clear();
rigidSupports.clear();
isRigidSupport.clear();
pMesh.reset();
plotYValues.clear();
plotHandle.reset();
@ -807,136 +807,113 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
// omp_lock_t writelock;
// omp_init_lock(&writelock);
//#ifdef ENABLE_OPENMP
//#pragma omp parallel for schedule(static) num_threads(8)
//#pragma omp parallel for schedule(static) num_threads(5)
//#endif
for (int ei = 0; ei < pMesh->EN(); ei++) {
const EdgeType &e = pMesh->edge[ei];
const SimulationMesh::VertexType &ev_j = *e.cV(0);
const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1);
const Element &element = pMesh->elements[e];
const Element::LocalFrame &ef = element.frame;
const VectorType t3CrossN_j = Cross(ef.t3, ev_j.cN());
const VectorType t3CrossN_jplus1 = Cross(ef.t3, ev_jplus1.cN());
const double theta1_j = ef.t1.dot(t3CrossN_j);
const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1);
const double theta2_j = ef.t2.dot(t3CrossN_j);
const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1);
const bool shouldBreak = mCurrentSimulationStep == 12970;
const double theta3_j = computeTheta3(e, ev_j);
const double theta3_jplus1 = computeTheta3(e, ev_jplus1);
// element.rotationalDisplacements_j.theta1 = theta1_j;
// element.rotationalDisplacements_jplus1.theta1 = theta1_jplus1;
// element.rotationalDisplacements_j.theta2 = theta2_j;
// element.rotationalDisplacements_jplus1.theta2 = theta2_jplus1;
// element.rotationalDisplacements_j.theta3 = theta3_j;
// element.rotationalDisplacements_jplus1.theta3 = theta3_jplus1;
std::vector<std::pair<int, Vector6d>> internalForcesContributionFromThisEdge(4,
{-1,
Vector6d()});
for (VertexIndex evi = 0; evi < 2; evi++) {
const SimulationMesh::VertexType &ev = *e.cV(evi);
const Node &edgeNode = pMesh->nodes[ev];
internalForcesContributionFromThisEdge[evi].first = edgeNode.vi;
std::for_each(
std::execution::par_unseq,
pMesh->edge.begin(),
pMesh->edge.end(),
[&](const EdgeType &e)
// for (int ei = 0; ei < pMesh->EN(); ei++)
{
const int ei = pMesh->getIndex(e);
// const EdgeType &e = pMesh->edge[ei];
const SimulationMesh::VertexType &ev_j = *e.cV(0);
const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1);
const Element &element = pMesh->elements[e];
const Element::LocalFrame &ef = element.frame;
const VectorType t3CrossN_j = Cross(ef.t3, ev_j.cN());
const VectorType t3CrossN_jplus1 = Cross(ef.t3, ev_jplus1.cN());
const double theta1_j = ef.t1.dot(t3CrossN_j);
const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1);
const double theta2_j = ef.t2.dot(t3CrossN_j);
const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1);
const bool shouldBreak = mCurrentSimulationStep == 12970;
const double theta3_j = computeTheta3(e, ev_j);
const double theta3_jplus1 = computeTheta3(e, ev_jplus1);
// element.rotationalDisplacements_j.theta1 = theta1_j;
// element.rotationalDisplacements_jplus1.theta1 = theta1_jplus1;
// element.rotationalDisplacements_j.theta2 = theta2_j;
// element.rotationalDisplacements_jplus1.theta2 = theta2_jplus1;
// element.rotationalDisplacements_j.theta3 = theta3_j;
// element.rotationalDisplacements_jplus1.theta3 = theta3_jplus1;
std::vector<std::pair<int, Vector6d>>
internalForcesContributionFromThisEdge(4, {-1, Vector6d()});
for (VertexIndex evi = 0; evi < 2; evi++) {
const SimulationMesh::VertexType &ev = *e.cV(evi);
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;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(axialForce_dofi)) {
std::cerr << "nan in axial" << evi << std::endl;
}
#endif
// Torsional force computation
const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi);
const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi);
const double theta1Diff = theta1_jplus1 - theta1_j;
const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv;
const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff
* element.rigidity.torsional;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(torsionalForce_dofi)) {
std::cerr << "nan in torsional" << evi << std::endl;
}
#endif
// First bending force computation
////theta2_j derivative
const double theta2_j_deriv = computeDerivativeTheta2(e, 0, evi, dofi);
////theta2_jplus1 derivative
const double theta2_jplus1_deriv = computeDerivativeTheta2(e, 1, evi, dofi);
////1st in bracket term
const double firstBendingForce_inBracketsTerm_0 = theta2_j_deriv * 2 * theta2_j;
////2nd in bracket term
const double firstBendingForce_inBracketsTerm_1 = theta2_jplus1_deriv
* theta2_j;
////3rd in bracket term
const double firstBendingForce_inBracketsTerm_2 = theta2_j_deriv
* theta2_jplus1;
////4th in bracket term
const double firstBendingForce_inBracketsTerm_3 = 2 * theta2_jplus1_deriv
* theta2_jplus1;
// 3rd term computation
const double firstBendingForce_inBracketsTerm
= firstBendingForce_inBracketsTerm_0 + firstBendingForce_inBracketsTerm_1
+ firstBendingForce_inBracketsTerm_2 + firstBendingForce_inBracketsTerm_3;
const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm
* element.rigidity.firstBending;
// Second bending force computation
////theta2_j derivative
const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui);
////theta2_jplus1 derivative
const double theta3_jplus1_deriv = computeDerivativeTheta3(e, ev_jplus1, dui);
////1st in bracket term
const double secondBendingForce_inBracketsTerm_0 = theta3_j_deriv * 2
* theta3_j;
////2nd in bracket term
const double secondBendingForce_inBracketsTerm_1 = theta3_jplus1_deriv
* theta3_j;
////3rd in bracket term
const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv
* theta3_jplus1;
////4th in bracket term
const double secondBendingForce_inBracketsTerm_3 = 2 * theta3_jplus1_deriv
* theta3_jplus1;
// 3rd term computation
const double secondBendingForce_inBracketsTerm
= secondBendingForce_inBracketsTerm_0 + secondBendingForce_inBracketsTerm_1
+ secondBendingForce_inBracketsTerm_2
+ secondBendingForce_inBracketsTerm_3;
double secondBendingForce_dofi = secondBendingForce_inBracketsTerm
* element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi].second[dofi]
= axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
+ torsionalForce_dofi;
}
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) {
const bool isDofConstrainedFor_refElemOtherVertex
= isVertexConstrained[refElemOtherVertexNode.vi]
&& fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi);
if (!isDofConstrainedFor_refElemOtherVertex) {
DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi};
////theta3_j derivative
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)) {
std::cerr << "nan in axial" << evi << std::endl;
}
#endif
// Torsional force computation
const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi);
const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi);
const double theta1Diff = theta1_jplus1 - theta1_j;
const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv;
const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff
* element.rigidity.torsional;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(torsionalForce_dofi)) {
std::cerr << "nan in torsional" << evi << std::endl;
}
#endif
// First bending force computation
////theta2_j derivative
const double theta2_j_deriv = computeDerivativeTheta2(e, 0, evi, dofi);
////theta2_jplus1 derivative
const double theta2_jplus1_deriv = computeDerivativeTheta2(e, 1, evi, dofi);
////1st in bracket term
const double firstBendingForce_inBracketsTerm_0 = theta2_j_deriv * 2
* theta2_j;
////2nd in bracket term
const double firstBendingForce_inBracketsTerm_1 = theta2_jplus1_deriv
* theta2_j;
////3rd in bracket term
const double firstBendingForce_inBracketsTerm_2 = theta2_j_deriv
* theta2_jplus1;
////4th in bracket term
const double firstBendingForce_inBracketsTerm_3 = 2 * theta2_jplus1_deriv
* theta2_jplus1;
// 3rd term computation
const double firstBendingForce_inBracketsTerm
= firstBendingForce_inBracketsTerm_0
+ firstBendingForce_inBracketsTerm_1
+ firstBendingForce_inBracketsTerm_2
+ firstBendingForce_inBracketsTerm_3;
const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm
* element.rigidity.firstBending;
// Second bending force computation
////theta2_j derivative
const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui);
////theta3_jplus1 derivative
////theta2_jplus1 derivative
const double theta3_jplus1_deriv = computeDerivativeTheta3(e,
ev_jplus1,
dui);
@ -950,25 +927,61 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv
* theta3_jplus1;
////4th in bracket term
const double secondBendingForce_inBracketsTerm_3 = theta3_jplus1_deriv * 2
const double secondBendingForce_inBracketsTerm_3 = 2 * theta3_jplus1_deriv
* theta3_jplus1;
// 4th term computation
// 3rd term computation
const double secondBendingForce_inBracketsTerm
= secondBendingForce_inBracketsTerm_0
+ secondBendingForce_inBracketsTerm_1
+ secondBendingForce_inBracketsTerm_2
+ secondBendingForce_inBracketsTerm_3;
const double secondBendingForce_dofi = secondBendingForce_inBracketsTerm
* element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi + 2].second[dofi]
= secondBendingForce_dofi;
double secondBendingForce_dofi = secondBendingForce_inBracketsTerm
* element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi].second[dofi]
= axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
+ torsionalForce_dofi;
}
if (edgeNode.referenceElement != &e) {
const bool isDofConstrainedFor_refElemOtherVertex
= isVertexConstrained[refElemOtherVertexNode.vi]
&& fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi);
if (!isDofConstrainedFor_refElemOtherVertex) {
DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi};
////theta3_j derivative
const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui);
////theta3_jplus1 derivative
const double theta3_jplus1_deriv = computeDerivativeTheta3(e,
ev_jplus1,
dui);
////1st in bracket term
const double secondBendingForce_inBracketsTerm_0 = theta3_j_deriv * 2
* theta3_j;
////2nd in bracket term
const double secondBendingForce_inBracketsTerm_1 = theta3_jplus1_deriv
* theta3_j;
////3rd in bracket term
const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv
* theta3_jplus1;
////4th in bracket term
const double secondBendingForce_inBracketsTerm_3 = theta3_jplus1_deriv
* 2 * theta3_jplus1;
// 4th term computation
const double secondBendingForce_inBracketsTerm
= secondBendingForce_inBracketsTerm_0
+ secondBendingForce_inBracketsTerm_1
+ secondBendingForce_inBracketsTerm_2
+ secondBendingForce_inBracketsTerm_3;
const double secondBendingForce_dofi = secondBendingForce_inBracketsTerm
* element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi + 2].second[dofi]
= secondBendingForce_dofi;
}
}
}
}
}
internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge;
}
internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge;
});
//#pragma omp parallel for schedule(static) num_threads(8)
@ -1251,7 +1264,7 @@ void DRMSimulationModel::updateNodalMasses()
pMesh->nodes[v].damping_6d[DoF::Nr] = 2
* std::sqrt(rotationalSumSk_I2
* pMesh->nodes[v].mass_6d[DoF::Nr]);
pMesh->nodes[v].damping_6d = pMesh->nodes[v].damping_6d * 1e-2;
pMesh->nodes[v].damping_6d = pMesh->nodes[v].damping_6d * 1e-3;
}
assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk
/ pMesh->nodes[v].mass.translational
@ -1296,8 +1309,8 @@ void DRMSimulationModel::updateNodalVelocities()
Node &node = pMesh->nodes[v];
if (mSettings.useViscousDamping) {
const Vector6d massOverDt = node.mass_6d / Dt;
const Vector6d visciousDampingFactor(viscuousDampingConstant / 2);
// const Vector6d &visciousDampingFactor = node.damping_6d;
// const Vector6d visciousDampingFactor(viscuousDampingConstant / 2);
const Vector6d &visciousDampingFactor = node.damping_6d;
const Vector6d denominator = massOverDt + visciousDampingFactor / 2;
const Vector6d firstTermNominator = massOverDt - visciousDampingFactor / 2;
const Vector6d firstTerm = node.velocity * firstTermNominator / denominator;
@ -2609,3 +2622,19 @@ mesh->currentTotalPotentialEnergykN*/
return results;
}
bool DRMSimulationModel::Settings::save(const filesystem::__cxx11::path &folderPath) const
{
bool returnValue = true;
std::filesystem::create_directories(folderPath);
nlohmann::json json;
json[jsonLabels.meshFilename]
= std::filesystem::relative(std::filesystem::path(meshFilename),
std::filesystem::path(jsonFilename).parent_path())
.string();
const std::string jsonFilename = "drmSettings.json";
std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename));
jsonFile << json;
}
bool DRMSimulationModel::Settings::load(const filesystem::__cxx11::path &filePath) const {}

2671
drmsimulationmodel.cpp.orig Executable file
View File

File diff suppressed because it is too large Load Diff

24
drmsimulationmodel.hpp
View File

@ -25,30 +25,32 @@ class DRMSimulationModel
public:
struct Settings
{
bool isDebugMode{false};
int debugModeStep{100000};
// bool isDebugMode{false};
std::optional<int> debugModeStep{0};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
int drawingStep{1};
double totalTranslationalKineticEnergyThreshold{1e-8};
double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
double residualForcesMovingAverageNormThreshold{1e-8};
// int drawingStep{0};
// double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
// double residualForcesMovingAverageNormThreshold{1e-8};
double Dtini{0.1};
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{50};
int desiredGradualExternalLoadsSteps{1};
std::optional<double> shouldUseTranslationalKineticEnergyThreshold;
// int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
bool useAverage{false};
double averageResidualForcesCriterionThreshold{1e-5};
std::optional<double> averageResidualForcesCriterionThreshold{1e-5};
Settings() {}
bool save(const std::filesystem::path &folderPath) const;
bool load(const std::filesystem::path &filePath) const;
bool useViscousDamping{false};
bool useKineticDamping{true};
struct JsonLabels
{};
};
private:

254
drmsimulationmodel.hpp.orig Executable file
View File

@ -0,0 +1,254 @@
#ifndef BEAMFORMFINDER_HPP
#define BEAMFORMFINDER_HPP
#include "simulationmesh.hpp"
#include "matplot/matplot.h"
#include "simulation_structs.hpp"
#include <Eigen/Dense>
#include <filesystem>
#include <unordered_set>
struct SimulationJob;
enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
using DoFType = int;
using EdgeType = VCGEdgeMesh::EdgeType;
using VertexType = VCGEdgeMesh::VertexType;
struct DifferentiateWithRespectTo {
const VertexType &v;
const DoFType &dofi;
};
class DRMSimulationModel
{
public:
struct Settings
{
bool isDebugMode{false};
int debugModeStep{100000};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
int drawingStep{1};
double totalTranslationalKineticEnergyThreshold{1e-8};
double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
double residualForcesMovingAverageNormThreshold{1e-8};
double Dtini{0.1};
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{50};
int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
bool useAverage{false};
double averageResidualForcesCriterionThreshold{1e-5};
Settings() {}
bool useViscousDamping{false};
bool useKineticDamping{true};
};
private:
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
bool shouldTemporarilyDampForces{false};
bool shouldTemporarilyAmplifyForces{true};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
int externalLoadStep{1};
<<<<<<< HEAD
const double viscuousDampingConstant{100};
=======
>>>>>>> master
std::vector<bool> isVertexConstrained;
std::vector<bool> isRigidSupport;
double minTotalResidualForcesNorm{std::numeric_limits<double>::max()};
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> pMesh;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> constrainedVertices;
SimulationHistory history;
// Eigen::Tensor<double, 4> theta3Derivatives;
// std::unordered_map<MyKeyType, double, key_hash> theta3Derivatives;
bool shouldApplyInitialDistortion = false;
void reset();
void updateNodalInternalForces(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateNodalExternalForces(
const std::unordered_map<VertexIndex, Vector6d> &nodalForces,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void updateResidualForces();
void updateRotationalDisplacements();
void updateElementalLengths();
void updateNodalMasses();
void updateNodalAccelerations();
void updateNodalVelocities();
void updateNodalDisplacements();
void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
Vector6d computeElementTorsionalForce(const Element &element,
const Vector6d &displacementDifference,
const std::unordered_set<DoFType> &constrainedDof);
// BeamFormFinder::Vector6d computeElementFirstBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
// BeamFormFinder::Vector6d computeElementSecondBendingForce(
// const Element &element, const Vector6d &displacementDifference,
// const std::unordered_set<gsl::index> &constrainedDof);
void updateKineticEnergy();
void resetVelocities();
SimulationResults computeResults(const std::shared_ptr<SimulationJob> &pJob);
void updateNodePosition(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
void applyDisplacements(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
#ifdef POLYSCOPE_DEFINED
void draw(const string &screenshotsFolder= {});
#endif
void
updateNodalInternalForce(Vector6d &nodalInternalForce,
const Vector6d &elementInternalForce,
const std::unordered_set<DoFType> &nodalFixedDof);
Vector6d computeElementInternalForce(
const Element &elem, const Node &n0, const Node &n1,
const std::unordered_set<DoFType> &n0ConstrainedDof,
const std::unordered_set<DoFType> &n1ConstrainedDof);
Vector6d computeElementAxialForce(const ::EdgeType &e) const;
VectorType computeDisplacementDifferenceDerivative(
const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
double
computeDerivativeElementLength(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT1(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType
computeDerivativeOfNormal(const VertexType &v,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT3(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
VectorType computeDerivativeT2(const EdgeType &e,
const DifferentiateWithRespectTo &dui) const;
double computeDerivativeTheta2(const EdgeType &e, const VertexIndex &evi,
const VertexIndex &dwrt_evi,
const DoFType &dwrt_dofi) const;
void updateElementalFrames();
VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
// bool isRigidSupport(const VertexType &v) const;
static double computeDerivativeOfNorm(const VectorType &x,
const VectorType &derivativeOfX);
static VectorType computeDerivativeOfCrossProduct(
const VectorType &a, const VectorType &derivativeOfA, const VectorType &b,
const VectorType &derivativeOfB);
double computeTheta3(const EdgeType &e, const VertexType &v);
double computeDerivativeTheta3(const EdgeType &e, const VertexType &v,
const DifferentiateWithRespectTo &dui) const;
double computeTotalPotentialEnergy();
void computeRigidSupports();
void updateNormalDerivatives();
void updateT1Derivatives();
void updateT2Derivatives();
void updateT3Derivatives();
void updateRDerivatives();
double computeDerivativeTheta1(const EdgeType &e, const VertexIndex &evi,
const VertexIndex &dwrt_evi,
const DoFType &dwrt_dofi) const;
// void updatePositionsOnTheFly(
// const std::unordered_map<VertexIndex,
// std::unordered_set<gsl::index>>
// &fixedVertices);
void updateResidualForcesOnTheFly(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updatePositionsOnTheFly(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void updateNodeNormal(
VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices);
void applyForcedNormals(
const std::unordered_map<VertexIndex, VectorType> nodalForcedRotations);
void printCurrentState() const;
void printDebugInfo() const;
double computeTotalInternalPotentialEnergy();
void applySolutionGuess(const SimulationResults &solutionGuess,
const std::shared_ptr<SimulationJob> &pJob);
void updateNodeNr(VertexType &v);
public:
DRMSimulationModel();
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings(),
const SimulationResults &solutionGuess = SimulationResults());
static void runUnitTests();
};
template <typename PointType> PointType Cross(PointType p1, PointType p2) {
return p1 ^ p2;
}
inline size_t currentStep{0};
inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei,
const DoFType &dofi) {
const size_t numberOfVertices = 10;
const VertexIndex middleNodeIndex = numberOfVertices / 2;
// return vi == middleNodeIndex && dofi == 1;
return dofi == 1 && ((vi == 1 && ei == 0) || (vi == 9 && ei == 9));
}
inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei) {
const size_t numberOfVertices = 10;
const VertexIndex middleNodeIndex = numberOfVertices / 2;
return (vi == middleNodeIndex);
// return (vi == 0 || vi == numberOfVertices - 1) && currentStep == 1;
return (vi == 1 && ei == 0) || (vi == 9 && ei == 9);
}
#endif // BEAMFORMFINDER_HPP

149
topologyenumerator.cpp
View File

@ -191,7 +191,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector<size_t> &reduced
patternGeometryAllEdges.getVertices(),
intersectingEdges,
validEdges);
// statistics.print(setupString, perEdgeResultPath);
statistics.print(setupString, perEdgeResultPath);
statistics.reset();
}
} else {
std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges << " edges."
@ -372,6 +373,20 @@ std::vector<vcg::Point2i> TopologyEnumerator::getValidEdges(
return validEdges;
}
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");
if (saveTilledPattern) {
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,
@ -418,6 +433,8 @@ void TopologyEnumerator::computeValidPatterns(
// std::string previousPatternBinaryRepresentation(validEdges.size(),'0');
size_t patternIndex = 0;
bool validPatternsExist = false;
const bool exportTilledPattern = false;
const bool saveCompressedFormat = false;
do {
patternIndex++;
const std::string patternName = std::to_string(patternIndex);
@ -438,6 +455,7 @@ void TopologyEnumerator::computeValidPatterns(
PatternGeometry patternGeometry;
patternGeometry.add(allVertices, patternEdges);
patternGeometry.setLabel(patternName);
// Check if pattern contains intersecting edges
const bool patternContainsIntersectingEdges
@ -448,21 +466,13 @@ void TopologyEnumerator::computeValidPatterns(
statistics.numberOfPatternsWithIntersectingEdges++;
if (debugIsOn) {
if (savePlyFiles) {
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry);
auto intersectingPatternsPath = std::filesystem::path(resultsPath)
.append("Intersecting");
std::filesystem::create_directory(intersectingPatternsPath);
patternGeometry.save(
std::filesystem::path(intersectingPatternsPath).append(patternName).string()
+ ".ply");
tiledPatternGeometry.save(std::filesystem::path(intersectingPatternsPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
exportPattern(std::filesystem::path(resultsPath).append("Intersecting"),
patternGeometry,
exportTilledPattern);
}
} else {
continue; // should be uncommented in order to improve performance
}
continue; // should be uncommented in order to improve performance
}
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
@ -470,18 +480,10 @@ void TopologyEnumerator::computeValidPatterns(
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
if (debugIsOn /*|| savePlyFiles*/) {
if (savePlyFiles) {
auto danglingEdgesPath = std::filesystem::path(resultsPath)
.append("ExceedingAngleThreshold");
std::filesystem::create_directory(danglingEdgesPath);
patternGeometry.save(
std::filesystem::path(danglingEdgesPath).append(patternName).string()
+ ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
exportPattern(std::filesystem::path(resultsPath)
.append("ExceedingAngleThreshold"),
patternGeometry,
exportTilledPattern);
}
} else {
continue;
@ -493,18 +495,9 @@ void TopologyEnumerator::computeValidPatterns(
if (tiledPatternHasNodeWithValenceGreaterThanDesired) {
if (debugIsOn) {
if (savePlyFiles) {
auto danglingEdgesPath = std::filesystem::path(resultsPath)
.append("HighValencePatterns");
std::filesystem::create_directory(danglingEdgesPath);
patternGeometry.save(
std::filesystem::path(danglingEdgesPath).append(patternName).string()
+ ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
auto highValencePath = std::filesystem::path(resultsPath)
.append("HighValencePatterns");
exportPattern(highValencePath, patternGeometry, exportTilledPattern);
}
} else {
continue;
@ -522,28 +515,19 @@ void TopologyEnumerator::computeValidPatterns(
// Check dangling edges with vcg method
// const bool vcg_tiledPatternHasDangling =
// tiledPatternGeometry.hasUntiledDanglingEdges();
if (tiledPatternHasDanglingEdges && !hasFloatingComponents /*&& !hasArticulationPoints*/) {
if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents && !hasArticulationPoints*/) {
statistics.numberOfPatternsWithADanglingEdgeOrNode++;
if (debugIsOn) {
if (savePlyFiles) {
auto danglingEdgesPath = std::filesystem::path(resultsPath).append("Dangling");
std::filesystem::create_directory(danglingEdgesPath);
patternGeometry.save(
std::filesystem::path(danglingEdgesPath).append(patternName).string()
+ ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
exportPattern(danglingEdgesPath, patternGeometry, exportTilledPattern);
}
} else {
continue;
}
}
if (hasFloatingComponents && !hasArticulationPoints && !tiledPatternHasDanglingEdges) {
if (hasFloatingComponents /*&& !hasArticulationPoints && !tiledPatternHasDanglingEdges */) {
statistics.numberOfPatternsWithMoreThanASingleCC++;
if (debugIsOn) {
if (savePlyFiles) {
@ -599,34 +583,25 @@ void TopologyEnumerator::computeValidPatterns(
assert(colorsRegistered == eCC[0].first);
tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
if (exportTilledPattern) {
tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
}
}
} else {
continue;
}
}
if (hasArticulationPoints && !hasFloatingComponents && !tiledPatternHasDanglingEdges) {
if (hasArticulationPoints /*&& !hasFloatingComponents && !tiledPatternHasDanglingEdges */) {
statistics.numberOfPatternsWithArticulationPoints++;
if (debugIsOn) {
if (savePlyFiles) {
auto articulationPointsPath = std::filesystem::path(resultsPath)
.append("ArticulationPoints");
std::filesystem::create_directory(articulationPointsPath);
patternGeometry.save(
std::filesystem::path(articulationPointsPath).append(patternName).string()
+ ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
tiledPatternGeometry.save(std::filesystem::path(articulationPointsPath)
.append(patternName + "_tiled")
.string()
+ ".ply");
// std::cout << "Pattern:" << patternName << std::endl;
exportPattern(articulationPointsPath, patternGeometry, exportTilledPattern);
}
} else {
continue;
@ -649,39 +624,27 @@ void TopologyEnumerator::computeValidPatterns(
statistics.numberOfValidPatterns++;
validPatternsExist = true;
if (savePlyFiles) {
// if (numberOfDesiredEdges == 4) {
// std::cout << "Saving:"
// << std::filesystem::path(validPatternsPath)
// .append(patternName)
// .string() +
// ".ply"
// << std::endl;
// }
patternGeometry.save(
std::filesystem::path(validPatternsPath).append(patternName).string() + ".ply");
PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(
patternGeometry); // the marked nodes of hasDanglingEdges are
tiledPatternGeometry.save(
std::filesystem::path(validPatternsPath).append(patternName + "_tiled").string()
+ ".ply");
exportPattern(validPatternsPath, patternGeometry, exportTilledPattern);
}
PatternIO::Pattern pattern;
pattern.edges = patternEdges;
pattern.name = patternIndex;
patternSet.patterns.emplace_back(pattern);
// Save valid patterns
// if (patternIndex% patternSetBufferSize == 0) {
if (statistics.numberOfValidPatterns % patternSetBufferSize == 0) {
PatternIO::save(compressedPatternsFilePath, patternSet);
patternSet.patterns.clear();
patternSet.patterns.reserve(patternSetBufferSize);
if (saveCompressedFormat) {
PatternIO::Pattern pattern;
pattern.edges = patternEdges;
pattern.name = patternIndex;
patternSet.patterns.emplace_back(pattern);
// Save valid patterns
// if (patternIndex% patternSetBufferSize == 0) {
if (statistics.numberOfValidPatterns % patternSetBufferSize == 0) {
PatternIO::save(compressedPatternsFilePath, patternSet);
patternSet.patterns.clear();
patternSet.patterns.reserve(patternSetBufferSize);
}
}
}
// assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges);
} while (std::next_permutation(patternBinaryRepresentation.begin(),
patternBinaryRepresentation.end()));
if (!patternSet.patterns.empty()) {
if (!patternSet.patterns.empty() && saveCompressedFormat) {
PatternIO::save(compressedPatternsFilePath, patternSet);
}

91
topologyenumerator.hpp
View File

@ -1,12 +1,12 @@
#ifndef TOPOLOGYENUMERATOR_HPP
#define TOPOLOGYENUMERATOR_HPP
#include "nlohmann/json.hpp"
#include "patternIO.hpp"
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <filesystem>
#include <fstream>
#include <iostream>
//#include <nlohmann/json.hpp>
#include <string>
#include <vector>
@ -63,28 +63,26 @@ class TopologyEnumerator {
size_t numberOfPatternsWithADanglingEdgeOrNode{0};
size_t numberOfPatternsWithArticulationPoints{0};
size_t numberOfValidPatterns{0};
// nlohmann::json convertToJson() const {
// nlohmann::json json;
// json["numPossibleEdges"] = numberOfPossibleEdges;
// json["numCoincideEdges"] = numberOfCoincideEdges;
// json["numDuplicateEdges"] = numberOfDuplicateEdges;
// json["numValidEdges"] = numberOfValidEdges;
// json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs;
// json["numPatterns"] = numberOfPatterns;
// // json["numIntersectingEdgesOverAllPatterns"] =
// // numberOfIntersectingEdgesOverAllPatterns;
// json["numPatternsWithIntersectingEdges"] =
// numberOfPatternsWithIntersectingEdges;
// json["numPatternsWithNotASingleCC"] =
// numberOfPatternsWithMoreThanASingleCC;
// json["numPatternsWithDangling"] =
// numberOfPatternsWithADanglingEdgeOrNode;
// json["numPatternsWithArticulationPoints"] =
// numberOfPatternsWithArticulationPoints;
// json["numValidPatterns"] = numberOfValidPatterns;
nlohmann::json convertToJson() const
{
nlohmann::json json;
json["numPossibleEdges"] = numberOfPossibleEdges;
json["numCoincideEdges"] = numberOfCoincideEdges;
json["numDuplicateEdges"] = numberOfDuplicateEdges;
json["numValidEdges"] = numberOfValidEdges;
json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs;
json["numPatterns"] = numberOfPatterns;
// json["numIntersectingEdgesOverAllPatterns"] =
// numberOfIntersectingEdgesOverAllPatterns;
json["numPatternsWithIntersectingEdges"] = numberOfPatternsWithIntersectingEdges;
json["numPatternsWithNotASingleCC"] = numberOfPatternsWithMoreThanASingleCC;
json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode;
json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints;
json["numValidPatterns"] = numberOfValidPatterns;
return json;
}
// return json;
// }
void print(const std::string &setupString,
const std::filesystem::path &directoryPath) const {
std::cout << "The setup " << setupString << std::endl;
@ -112,20 +110,37 @@ class TopologyEnumerator {
<< " patterns found with a dangling node or edge" << std::endl;
std::cout << numberOfPatternsWithArticulationPoints
<< " patterns found with an articulation point" << std::endl;
std::cout << numberOfValidPatterns << " valid patterns were found"
<< std::endl;
// if (!directoryPath.empty()) {
// auto json = convertToJson();
std::cout << numberOfValidPatterns << " valid patterns were found" << std::endl;
if (!directoryPath.empty()) {
auto json = convertToJson();
// std::ofstream file;
// file.open(std::filesystem::path(directoryPath)
// .append("statistics.csv")
// .string());
// file << "setup," << setupString << "\n";
// for (const auto &el : json.items()) {
// file << el.key() << "," << el.value() << "\n";
// }
// }
std::ofstream file;
file.open(std::filesystem::path(directoryPath).append("statistics.csv").string());
file << "setup," << setupString << "\n";
for (const auto &el : json.items()) {
file << el.key() << ",";
}
file << "\n";
for (const auto &el : json.items()) {
file << el.value() << ",";
}
file << "\n";
}
}
void reset()
{
numberOfPossibleEdges = 0;
numberOfCoincideEdges = 0;
numberOfDuplicateEdges = 0;
numberOfValidEdges = 0;
numberOfIntersectingEdgePairs = 0;
numberOfPatterns = 0;
numberOfPatternsWithIntersectingEdges = 0;
numberOfPatternsWithMoreThanASingleCC = 0;
numberOfPatternsWithADanglingEdgeOrNode = 0;
numberOfPatternsWithArticulationPoints = 0;
numberOfValidPatterns = 0;
}
};
@ -173,8 +188,10 @@ private:
const size_t &numberOfDesiredEdges,
const std::filesystem::path &resultsPath,
const std::vector<vcg::Point3d> &vertices,
const std::unordered_map<size_t, std::unordered_set<size_t>>
&intersectingEdges,
const std::unordered_map<size_t, std::unordered_set<size_t>> &intersectingEdges,
const std::vector<vcg::Point2i> &validEdges);
void exportPattern(const std::filesystem::path &saveToPath,
PatternGeometry &patternGeometry,
const bool saveTilledPattern) const;
};
#endif // TOPOLOGYENUMERATOR_HPP