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

Changed default normal. Identical cantilever results with panettaD

This commit is contained in:
iasonmanolas 2022-07-12 13:00:11 +03:00
parent 9d9f7dbacf
commit e70b4e9863
3 changed files with 255 additions and 170 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
CMakeLists.txt
View File

@ -48,10 +48,10 @@ if(${USE_POLYSCOPE})
add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR})
add_compile_definitions(POLYSCOPE_DEFINED)
target_sources(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/deps/imgui/imgui/misc/cpp/imgui_stdlib.cpp)
if(${USE_POLYSCOPE})
message("Using polyscope here")
message("Using polyscope..")
target_include_directories(${PROJECT_NAME} PUBLIC ${POLYSCOPE_SOURCE_DIR}/include)
target_link_libraries(${PROJECT_NAME} PUBLIC polyscope)
endif()
endif()
##vcglib devel branch
@ -62,9 +62,6 @@ download_project(PROJ vcglib_devel
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${vcglib_devel_SOURCE_DIR} ${vcglib_devel_BINARY_DIR})
#add_subdirectory("/home/iason/Coding/Libraries/vcglib" ${CMAKE_CURRENT_BINARY_DIR}/vcglib)
#target_include_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/vcglib")
target_sources(${PROJECT_NAME} PUBLIC ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
##matplot++ lib
@ -90,27 +87,11 @@ download_project(PROJ threed-beam-fea
)
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##TBB
set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb)
download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${TBB_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
link_directories(${TBB_BINARY_DIR})
###Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED)
find_package(Eigen3 3.4 REQUIRED)
if(MSVC)
add_compile_definitions(_HAS_STD_BYTE=0)
endif(MSVC)
#link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs)
#set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
if(${MYSOURCES_STATIC_LINK})
message("STATIC LINK MY SOURCES:" ${MYSOURCES_STATIC_LINK})
endif()
@ -149,33 +130,56 @@ download_project(PROJ ENSMALLEN
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
# add_subdirectory(${ENSMALLEN_SOURCE_DIR} ${ENSMALLEN_BINARY_DIR})
# target_link_libraries(${PROJECT_NAME} INTERFACE ensmallen)
target_include_directories(${PROJECT_NAME}
PUBLIC ${ENSMALLEN_SOURCE_DIR}/include)
add_compile_definitions(USE_ENSMALLEN)
#endif()
#Chrono
#set(Chrono_DIR "/home/iason/Coding/Libraries/chrono-7.0.3/build/cmake")
set(Chrono_DIR "/home/iason/Coding/Libraries/chrono-7.0.3/build/cmake")
LIST(APPEND CMAKE_PREFIX_PATH "${CMAKE_INSTALL_PREFIX}/../Chrono/lib")
#set(Chrono_DIR "/home/iason/Coding/build/external dependencies/CHRONO-src/cmake")
find_package(Chrono #[[COMPONENTS FEA]]
CONFIG)
if (NOT Chrono_FOUND)
message("Could not find Chrono or one of its required modules")
return()
#LIST(APPEND CMAKE_PREFIX_PATH "${CMAKE_INSTALL_PREFIX}/../Chrono/lib")
#find_package(Chrono #[[COMPONENTS FEA]]
# CONFIG)
#if (NOT Chrono_FOUND)
# message("Could not find Chrono or one of its required modules")
# return()
#endif()
#set_target_properties(${PROJECT_NAME} PROPERTIES
# COMPILE_FLAGS "${CHRONO_CXX_FLAGS} ${EXTRA_COMPILE_FLAGS}"
# COMPILE_DEFINITIONS "CHRONO_DATA_DIR=\"${CHRONO_DATA_DIR}\""
# LINK_FLAGS "${CHRONO_LINKER_FLAGS}")
#target_link_libraries(${PROJECT_NAME} PRIVATE ${CHRONO_LIBRARIES})
target_include_directories(${PROJECT_NAME} PRIVATE ${CHRONO_INCLUDE_DIRS} #[["/home/iason/Coding/Libraries/chrono-7.0.3/src"]] #[["/usr/include/irrlicht"]] )
#target_link_directories(${PROJECT_NAME} PUBLIC "/home/iason/Coding/Libraries/chrono-7.0.3/build/Debug/lib/")
target_link_libraries(${PROJECT_NAME} PRIVATE "/home/iason/Coding/Libraries/chrono-7.0.3/build/Debug/lib/libChronoEngine.so" #[[ChronoEngine_irrlicht]])
#add_subdirectory("/home/iason/Coding/Projects/Approximating shapes with flat patterns/Elastic rods/ElasticRods")
#set(MESHFEM_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/meshfem)
#download_project(PROJ MESHFEM
# GIT_REPOSITORY https://github.com/MeshFEM/MeshFEM.git
# GIT_TAG 3847ffaeb6e9c4bbe15f0f81f5c773c7c481c059
# BINARY_DIR ${MESHFEM_BINARY_DIR}
# PREFIX ${EXTERNAL_DEPS_DIR}
# ${UPDATE_DISCONNECTED_IF_AVAILABLE}
#)
#add_subdirectory(${MESHFEM_SOURCE_DIR} ${MESHFEM_BINARY_DIR})
#target_link_libraries(${PROJECT_NAME} PUBLIC MeshFEM)
##target_include_directories(${PROJECT_NAME}
##PUBLIC ${MESHFEM_SOURCE_DIR}/include)
##add_compile_definitions(USE_ENSMALLEN)
##TBB
if(NOT TARGET tbb)
set(TBB_BINARY_DIR ${CMAKE_CURRENT_BINARY_DIR}/tbb)
download_project(PROJ TBB
GIT_REPOSITORY https://github.com/wjakob/tbb.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
BINARY_DIR ${TBB_BINARY_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF)
add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR})
endif()
set_target_properties(${PROJECT_NAME} PROPERTIES
COMPILE_FLAGS "${CHRONO_CXX_FLAGS} ${EXTRA_COMPILE_FLAGS}"
COMPILE_DEFINITIONS "CHRONO_DATA_DIR=\"${CHRONO_DATA_DIR}\""
LINK_FLAGS "${CHRONO_LINKER_FLAGS}")
target_link_libraries(${PROJECT_NAME} PRIVATE ${CHRONO_LIBRARIES})
target_include_directories(${PROJECT_NAME} PRIVATE ${CHRONO_INCLUDE_DIRS})
link_directories(${TBB_BINARY_DIR})
#find_package(OpenMP REQUIRED)
#target_link_libraries(${PROJECT_NAME} OpenMP::OpenMP_CXX)
#target_include_directories(${PROJECT_NAME} OpenMP::OpenMP_CXX)
#message(STATUS "OpenMP_C_INCLUDE_DIRS: ${OpenMP_CXX_INCLUDE_DIRS}")

238
der_leimer.cpp
View File

@ -2,6 +2,7 @@
#include <unordered_set>
#define MAT(mat, m, row, col) mat[m * col + row]
//#define EXPLICIT_BC
DER_leimer::DER_leimer() {}
@ -255,11 +256,12 @@ SimulationResults DER_leimer::constructSimulationResults(const std::shared_ptr<S
// pJob->pMesh->registerForDrawing();
// polyscope::show();
for (int vi = 0; vi < config->numVertices(); vi++) {
std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
<< " " << config->vertices(vi, 2) << std::endl;
// std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
// << " " << config->vertices(vi, 2) << std::endl;
for (int dofi = 0; dofi < 3; dofi++) {
const double initialPos_dofi = pJob->pMesh->vert[vi].cP()[dofi];
const double displacement_dofi = config->vertices(vi, dofi) - initialPos_dofi;
const double deformedPos_dofi = config->vertices(vi, dofi);
const double displacement_dofi = deformedPos_dofi - initialPos_dofi;
simulationResults.displacements[vi][dofi] = displacement_dofi;
}
for (int dofi = 3; dofi < 6; dofi++) {
@ -353,7 +355,7 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptr<Simulation
// std::vector<Anchor> anchors{a0};
const std::shared_ptr<VCGEdgeMesh> pMesh = pJob->pMesh;
const int numRods = pMesh->EN();
const int numVertsPerRod = 2;
constexpr int numVertsPerRod = 2;
std::vector<int> vInd;
vInd.reserve(pMesh->EN());
for (int ei = 0; ei < pMesh->EN(); ei++) {
@ -391,12 +393,12 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptr<Simulation
a.bary = edgeVi;
a.position = ep->cP(edgeVi).ToEigenVector<Eigen::Vector3d>();
a.normal = ep->cV(edgeVi)->cN().ToEigenVector<Eigen::Vector3d>();
a.k_pos = 1e4;
a.k_pos = 1e6;
if (bc.second == fixedBC) {
a.k_dir = 0;
} else {
assert(bc.second == rigidBC);
a.k_dir = 1e4;
a.k_dir = 1e6;
}
anchors.push_back(a);
}
@ -569,7 +571,7 @@ bool DER_leimer::simulateOneStepGrid(RodConfig *config)
iter++;
if (iter > 50 /*&& forceResidual < 1e-5*/) {
if (iter > 100 /*|| solverResidual < 1e-10*/ || forceResidual < 1e-5) {
// for (int vi = 0; vi < config->numVertices(); vi++) {
// std::cout << vi << " pos:" << config->vertices(vi, 0) << " " << config->vertices(vi, 1)
// << " " << config->vertices(vi, 2) << std::endl;
@ -594,11 +596,11 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
std::vector<Eigen::Triplet<int>> *Hu_index,
std::vector<double> *Hu_value)
{
const bool computeHessian = [&]() {
if (H_index && H_value && Ju && Hu_index && Hu_value)
return true;
return false;
}();
// const bool computeHessian = [&]() {
// if (H_index && H_value && Ju && Hu_index && Hu_value)
// return true;
// return false;
// }();
const bool hasExternalForces = externalForces.rows() != 0;
@ -624,8 +626,10 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
nterms += config.rods[i]->numSegments(); // stretching
nterms += 2 * njoints; // bending
nterms += njoints; // twisting
#ifndef EXPLICIT_BC
if (params.anchorPoints) //??
nterms += 3 * config.rods[i]->numVertices();
#endif
// if (hasExternalForces) {
// nterms += config.rods[i]->numVertices();
// }
@ -635,12 +639,14 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
}
ndofs += config.numVertices() * 3; // DOF for vertices
//if (params.anchorPoints)
// nterms += 3 * config.numVertices();
if (params.anchorPoints)
nterms += 3 * config.numVertices();
//if (params.gravityEnabled)
// nterms += config.numVertices(); // floor force
#ifndef EXPLICIT_BC
nterms += 4 * config.numAnchors();
#endif
// Barycentric Coordinates are handled as DOF's since they can change when rods are sliding
int baryoffset = ndofs;
@ -669,8 +675,8 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
ndofs += 3 * config.numIntersections(); // euler angles as DOFs for intersections
r.resize(nterms);
r.setConstant(std::numeric_limits<double>::infinity());
// r.setZero();
// r.setConstant(std::numeric_limits<double>::infinity());
r.setZero();
gravityPE = 0;
gravityForces.resize(ndofs);
@ -864,16 +870,25 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
// const Eigen::Vector3d debug_forceLeimer = jacobianEntry * segr;
// assert(debug_forceLeimer == debug_forcePanetta_j);
for (int j = 0; j < 3; j++) {
// Gradients for dEs/dx_i and dEs/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m1)) {
// Gradients for dEs/dx_i and dEs/dx_(i+1)
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(
Eigen::Triplet<double>(roffset + si_rod, 3 * m1 + j, jacobianEntry[j]));
// }
// if (!pJob->constrainedVertices.contains(m2)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m2)) {
#endif
J.push_back(
Eigen::Triplet<double>(roffset + si_rod, 3 * m2 + j, -jacobianEntry[j]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
Eigen::Matrix3d Hlen = (1 / len) * Eigen::Matrix<double, 3, 3>::Identity()
- (1 / (len * len * len)) * (v1 - v2)
@ -945,6 +960,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
uistiffoffset + 4 * nsegsTotal + dofoffset + si_rod,
df_dt * (len - restlen)));
}
#endif
}
}
roffset += rod.numSegments();
@ -1067,18 +1083,28 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
// const Eigen::Vector3d debug_force1Leimer_xip1 = dEb1_dxip1 * debug_vertr1;
for (int j = 0; j < 3; j++) {
// Gradients for dEb1/dx_(i-1), dEb1/dx_i and dEb1/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m0)) {
J.push_back(
Eigen::Triplet<double>(roffset + 2 * idx, 3 * m0 + j, dEb1_dxim1[j]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + 2 * idx,
3 * m0 + j,
dEb1_dxim1[j]));
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(
Eigen::Triplet<double>(roffset + 2 * idx, 3 * m1 + j, dEb1_dxi[j]));
// }
// if (!pJob->constrainedVertices.contains(m2)) {
J.push_back(
Eigen::Triplet<double>(roffset + 2 * idx, 3 * m2 + j, dEb1_dxip1[j]));
// }
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m2)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + 2 * idx,
3 * m2 + j,
dEb1_dxip1[j]));
#ifdef EXPLICIT_BC
}
#endif
}
// Gradients for dk2/dx_(i-1) and dk2/dx_i
Eigen::Vector3d dk21 = 1.0 / (v1 - v0).norm()
@ -1100,20 +1126,29 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
// const Eigen::Vector3d debug_force2Leimer_xip1 = dEb2_dxip1 * debug_vertr2;
for (int j = 0; j < 3; j++) {
// Gradients for dEb2/dx_(i-1), dEb2/dx_i and dEb2/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m0)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
3 * m0 + j,
dEb2_dxim1[j]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
J.push_back(
Eigen::Triplet<double>(roffset + 2 * idx + 1, 3 * m1 + j, dEb2_dxi[j]));
// }
// if (!pJob->constrainedVertices.contains(m2)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
3 * m1 + j,
dEb2_dxi[j]));
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m2)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + 2 * idx + 1,
3 * m2 + j,
dEb2_dxip1[j]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
// Material curvature gradients dd/dtheta1 and dd/theta2
@ -1146,6 +1181,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
0.5 * sqrt(factor2 /* * width * width * width*/)
* Dd12.dot(kb)));
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
double chi = 1 + t01.dot(t12);
double len01 = (v1 - v0).norm();
@ -2014,6 +2050,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
int stop = 0;
}
}
#endif
}
} else {
r[roffset + 2 * idx] = 0;
@ -2083,19 +2120,27 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
const Eigen::Vector3d debug_leimerTwistForce_pos2 = jacob_theta2
* r[roffset + idx];
for (int j = 0; j < 3; j++) {
// if (!pJob->constrainedVertices.contains(m0)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(
Eigen::Triplet<double>(roffset + idx, 3 * m0 + j, -jacob_theta1[j]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + idx,
3 * m1 + j,
jacob_theta1[j] - jacob_theta2[j]));
// }
// if (!pJob->constrainedVertices.contains(m2)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m2)) {
#endif
J.push_back(
Eigen::Triplet<double>(roffset + idx, 3 * m2 + j, jacob_theta2[j]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
// Gradients dEt/dtheta1 and dEt/dtheta2
const double debug_leimerTwistForce_theta1 = sqrt(debug_factor)
@ -2111,6 +2156,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
dofoffset + thetaoffset + si_rod,
sqrt(debug_factor)));
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
double chi = 1 + t01.dot(t12);
double len01 = (v1 - v0).norm();
@ -2374,6 +2420,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
int stop = 0;
}
}
#endif
}
} else
r[roffset + idx] = 0;
@ -2384,6 +2431,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
dofoffset += rod.numSegments();
}
#ifndef EXPLICIT_BC
// mesh sliding energy Eslide
// allows ribbons to slide tangentially on mesh surface they try to approximate
// see paper Section 7.1
@ -2566,6 +2614,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
}
}
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
// design variables - Ui anchor positions
for (int j = 0; j < 3; j++) {
@ -2637,6 +2686,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
int stop = 0;
}
}
#endif
}
// direction constraint for anchor
@ -2678,6 +2728,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
rod1offset + a.seg,
sqrt(factor) * dtheta1.dot(d1dtheta)));
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
double len = (vSeg1 - vSeg0).norm();
Eigen::Matrix3d dt_de1 = ((1 / len) * Eigen::Matrix<double, 3, 3>::Identity()
@ -2971,6 +3022,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
}
int stop = 0;
}
#endif
}
idx += 4;
@ -2978,6 +3030,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
roffset += idx;
}
#endif
// compute energy terms for rod intersections
int incCount = 0;
@ -3119,16 +3172,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
* (-k1 * ttilde + st01.cross(d2tilde));
for (int k = 0; k < 3; k++) {
// Gradients for dEb1/dx_(i-1), dEb1/dx_i and dEb1/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m0)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount,
m0 * 3 + k,
-sqrt(factor1) * dk11[k]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount,
m1 * 3 + k,
sqrt(factor1) * dk11[k]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
J.push_back(Eigen::Triplet<double>(roffset + incCount,
euleroffset + i * 3 + 0,
@ -3145,16 +3204,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
* (-k2 * ttilde + st01.cross(d1tilde));
for (int k = 0; k < 3; k++) {
// Gradients for dEb2/dx_(i-1), dEb2/dx_i and dEb2/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m0)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
m0 * 3 + k,
-sqrt(factor2) * dk21[k]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
m1 * 3 + k,
sqrt(factor2) * dk21[k]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
J.push_back(Eigen::Triplet<double>(roffset + incCount + 1,
euleroffset + i * 3 + 0,
@ -3183,16 +3248,22 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
Eigen::Vector3d dpsi_dx = 0.5 * kb / (cv1 - cv0).norm();
for (int k = 0; k < 3; k++) {
// Gradients dEt/dx_(i-1), dEt/dx_i and dEt/dx_(i+1)
// if (!pJob->constrainedVertices.contains(m0)) {
#ifdef EXPLICIT_BC
if (!pJob->constrainedVertices.contains(m0)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
m0 * 3 + k,
-sqrt(factor3) * dpsi_dx[k]));
// }
// if (!pJob->constrainedVertices.contains(m1)) {
#ifdef EXPLICIT_BC
}
if (!pJob->constrainedVertices.contains(m1)) {
#endif
J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
m1 * 3 + k,
sqrt(factor3) * dpsi_dx[k]));
// }
#ifdef EXPLICIT_BC
}
#endif
}
// Gradients dEt/dtheta1 and dEt/dtheta2
J.push_back(Eigen::Triplet<double>(roffset + incCount + 2,
@ -3283,6 +3354,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
}
}
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
double chi = 1 + ct01.dot(st01);
double len01 = (cv1 - cv0).norm();
@ -4078,6 +4150,7 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
H_value->push_back(-sqrt(factor3*width) * dpsi22);
*/
}
#endif
}
incCount += 3;
}
@ -4089,12 +4162,31 @@ void DER_leimer::rAndJGrid(const RodConfig &config,
exit(-1);
if (Jr) {
// std::vector<Eigen::Triplet<double>> J_nonConstrained;
// std::unordered_set<int> constrainedDof_global;
// for (const auto &constrainedViDof : pJob->constrainedVertices) {
// const int vi = constrainedViDof.first;
// std::unordered_set<int> &constrainedDof_local = constrainedViDof.second;
// }
// for (auto &triplet : J) {
// }
Jr->setFromTriplets(J.begin(), J.end());
}
#ifdef COMPUTE_HESSIAN
if (computeHessian) {
Ju->resize(nterms, udofs);
Ju->setFromTriplets(Jui.begin(), Jui.end());
}
#endif
}
Eigen::Vector3d DER_leimer::getPerpendicularVector(const Eigen::Vector3d &t) const
{
Eigen::Vector3d candidate1 = Eigen::Vector3d(1, 0, 0).cross(t),
candidate2 = Eigen::Vector3d(0, 1, 0).cross(t);
return (candidate1.norm() > candidate2.norm()) ? candidate1.normalized()
: candidate2.normalized();
}
DER_leimer::RodConfig *DER_leimer::readRodGrid(
@ -4203,29 +4295,8 @@ DER_leimer::RodConfig *DER_leimer::readRodGrid(
const Eigen::Vector3d v1 = rodState.centerlinePositions.row(
(si_rod + 1) % numVertices_rod); //the % numVertices_rod is not needed
const Eigen::Vector3d e = (v1 - v0).normalized();
Eigen::Vector3d segmentDirector = Eigen::Vector3d(0, 0, 1);
/*= ((pMesh->edge[rodIndex].cV(0)->cN() + pMesh->edge[rodIndex].cV(1)->cN()) / 2)
.ToEigenVector<Eigen::Vector3d>()
.normalized();*/
// = (e.cross(Eigen::Vector3d(0, 0, 1).normalized())).normalized();
// switch (rodIndex) {
// case 0:
// segmentDirector = Eigen::Vector3d(1, 0, 0);
// break;
// case 1:
// segmentDirector = Eigen::Vector3d(0, -1, 0);
// break;
// case 2:
// segmentDirector = Eigen::Vector3d(-1, 0, 0);
// segmentDirector = Eigen::Vector3d(-1, 0, 0);
// break;
// default:
// assert(false);
// }
rodState.directors.row(si_rod) = segmentDirector;
// rodState.directors.row(si_rod) = getPerpendicularVector(e);
rodState.directors.row(si_rod) = Eigen::Vector3d(0, 0, 1);
// rodState.directors.row(si_rod) = binorms[si_rod];
// for (int k = 0; k < 3; k++) {
@ -4745,3 +4816,10 @@ int DER_leimer::RodConfig::getNumDoF() const
ndofs += 3 * numIntersections(); // euler angles as DOFs for intersections
return ndofs;
}
void DER_leimer::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
{
std::cout << "This function is currently not implemented" << std::endl;
assert(false);
std::terminate();
}

3
der_leimer.hpp
View File

@ -139,6 +139,8 @@ public:
int getNumDoF() const;
};
void setStructure(const std::shared_ptr<SimulationMesh> &pMesh) override;
private:
bool simulateOneStepGrid(RodConfig *config);
@ -181,6 +183,7 @@ private:
SimulationResults constructSimulationResults(const std::shared_ptr<SimulationJob> &pJob,
RodConfig *config) const;
std::shared_ptr<SimulationJob> pJob;
Eigen::Vector3d getPerpendicularVector(const Eigen::Vector3d &t) const;
};
Eigen::Matrix3d eulerRotation(double e0, double e1, double e2, int deriv = -1, int deriv2 = -1);