Tbb,refactoring
This commit is contained in:
iasonmanolas
parent
846289448b
commit
90a00c5415
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -222,6 +222,7 @@ void DRMSimulationModel::reset()
|
|||
mCurrentSimulationStep = 0;
|
||||
history.clear();
|
||||
constrainedVertices.clear();
|
||||
isRigidSupport.clear();
|
||||
pMesh.reset();
|
||||
plotYValues.clear();
|
||||
plotHandle.reset();
|
||||
|
|
@ -808,134 +809,111 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
|
|||
//#ifdef ENABLE_OPENMP
|
||||
//#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);
|
||||
|
|
@ -949,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)
|
||||
|
||||
|
|
@ -1247,17 +1261,7 @@ void DRMSimulationModel::updateNodalAccelerations()
|
|||
for (VertexType &v : pMesh->vert) {
|
||||
Node &node = pMesh->nodes[v];
|
||||
const VertexIndex vi = pMesh->getIndex(v);
|
||||
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
|
||||
if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) {
|
||||
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.translational;
|
||||
} else if (dofi == DoF::Nx) {
|
||||
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalJ;
|
||||
} else if (dofi == DoF::Ny) {
|
||||
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI3;
|
||||
} else if (dofi == DoF::Nr) {
|
||||
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2;
|
||||
}
|
||||
}
|
||||
node.acceleration = node.force.residual / node.mass_6d;
|
||||
#ifdef POLYSCOPE_DEFINED
|
||||
if (std::isnan(node.acceleration.norm())) {
|
||||
std::cout << "acceleration " << vi << ":" << node.acceleration.toString() << std::endl;
|
||||
|
|
@ -2578,3 +2582,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 {}
|
||||
|
|
|
|||
|
|
@ -25,28 +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:
|
||||
|
|
@ -158,9 +162,9 @@ private:
|
|||
|
||||
void updateElementalFrames();
|
||||
|
||||
VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const;
|
||||
VectorType computeDerivativeOfR(const EdgeType &e,
|
||||
const DifferentiateWithRespectTo &dui) const;
|
||||
|
||||
// bool isRigidSupport(const VertexType &v) const;
|
||||
|
||||
static double computeDerivativeOfNorm(const VectorType &x,
|
||||
const VectorType &derivativeOfX);
|
||||
|
|
|
|||
|
|
@ -148,6 +148,8 @@ struct Node {
|
|||
};
|
||||
|
||||
Mass mass;
|
||||
Vector6d mass_6d;
|
||||
Vector6d damping_6d;
|
||||
VertexIndex vi;
|
||||
CoordType initialLocation;
|
||||
CoordType initialNormal;
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -57,6 +57,15 @@ struct Vector6d : public std::array<double, 6> {
|
|||
return result;
|
||||
}
|
||||
|
||||
Vector6d operator/(const Vector6d &v) const
|
||||
{
|
||||
Vector6d result;
|
||||
for (size_t i = 0; i < 6; i++) {
|
||||
result[i] = this->operator[](i) / v[i];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
Vector6d operator+(const Vector6d &v) const {
|
||||
Vector6d result;
|
||||
for (size_t i = 0; i < 6; i++) {
|
||||
|
|
|
|||
Loading…
Reference in New Issue