suzanne.soy/FreeCAD
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

only use scale value for rotation-only solving, not the transformed clusters. This gives the more expected results

Browse Source
This commit is contained in:
Stefan Tröger 2013-12-22 12:11:01 +01:00
parent 66f9afbef3
commit 0a750b1bc5
3 changed files with 45 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp
View File

@ -21,6 +21,7 @@
#define DCM_KERNEL_IMP_H
#include "../kernel.hpp"
#include <Base/Console.h>
#include <boost/math/special_functions.hpp>
#include <Eigen/QR>
@ -62,9 +63,7 @@ void Dogleg<Kernel>::calculateStep(const Eigen::MatrixBase<Derived>& g, const Ei
const typename Kernel::Vector h_gn = jacobi.fullPivLu().solve(-residual);
const double eigen_error = (jacobi*h_gn + residual).norm();
#ifdef USE_LOGGING
BOOST_LOG(log)<< "Gauss Newton error: "<<eigen_error;
if(!boost::math::isfinite(h_gn.norm())) {
if(!boost::math::isfinite(h_gn.norm())) {
BOOST_LOG(log)<< "Unnormal gauss-newton detected: "<<h_gn.norm();
}
@ -81,9 +80,6 @@ void Dogleg<Kernel>::calculateStep(const Eigen::MatrixBase<Derived>& g, const Ei
// compute the dogleg step
if(h_gn.norm() <= delta) {
h_dl = h_gn;
#ifdef USE_LOGGING
BOOST_LOG(log)<< "Gauss Newton Step: "<<eigen_error;
#endif
}
else if((alpha*h_sd).norm() >= delta) {
//h_dl = alpha*h_sd;
@ -94,8 +90,6 @@ void Dogleg<Kernel>::calculateStep(const Eigen::MatrixBase<Derived>& g, const Ei
BOOST_LOG(log)<< "Unnormal dogleg descent detected: "<<h_dl.norm();
}
BOOST_LOG(log)<< "Steepest Descend step";
#endif
}
else {
@ -131,8 +125,6 @@ void Dogleg<Kernel>::calculateStep(const Eigen::MatrixBase<Derived>& g, const Ei
BOOST_LOG(log)<< "Unnormal dogleg beta detected: "<<beta;
}
BOOST_LOG(log)<< "Dogleg Step";
#endif
}
};
@ -421,21 +413,22 @@ void Kernel<Scalar, Nonlinear>::MappedEquationSystem::setAccess(AccessType t) {
if(t==complete) {
new(&Jacobi) MatrixMap(&m_jacobi(0,0),m_eqns,m_params,DynStride(m_eqns,1));
new(&Parameter) VectorMap(&m_parameter(0),m_params,DynStride(1,1));
new(&Residual) VectorMap(&m_residual(0), m_eqns, DynStride(1,1));
new(&Residual) VectorMap(&m_residual(0), m_eqns, DynStride(1,1));
}
else if(t==rotation) {
int num = m_param_trans_offset;
new(&Jacobi) MatrixMap(&m_jacobi(0,0),m_eqn_trans_offset,num,DynStride(m_eqns,1));
new(&Parameter) VectorMap(&m_parameter(0),num,DynStride(1,1));
new(&Residual) VectorMap(&m_residual(0), m_eqn_trans_offset, DynStride(1,1));
new(&Residual) VectorMap(&m_residual(0), m_eqn_trans_offset, DynStride(1,1));
}
else if(t==general) {
int num = m_params - m_param_trans_offset;
int eq_num = m_eqns - m_eqn_trans_offset;
new(&Jacobi) MatrixMap(&m_jacobi(m_eqn_trans_offset,m_param_trans_offset),eq_num,num,DynStride(m_eqns,1));
new(&Parameter) VectorMap(&m_parameter(m_param_trans_offset),num,DynStride(1,1));
new(&Residual) VectorMap(&m_residual(m_eqn_trans_offset), eq_num, DynStride(1,1));
new(&Residual) VectorMap(&m_residual(m_eqn_trans_offset), eq_num, DynStride(1,1));
}
m_access = t;
};
@ -443,11 +436,11 @@ template<typename Scalar, template<class> class Nonlinear>
bool Kernel<Scalar, Nonlinear>::MappedEquationSystem::hasAccessType(AccessType t) {
if(t==rotation)
return (m_param_rot_offset>0);
return (m_param_rot_offset>0 && m_eqn_rot_offset>0);
else if(t==general)
return (m_param_trans_offset<m_params);
return (m_param_trans_offset<m_params && m_eqn_trans_offset<m_eqns);
else
return (m_params>0);
return (m_params>0 && m_eqns>0);
};
template<typename Scalar, template<class> class Nonlinear>

45
src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp
View File

@ -107,12 +107,12 @@ SystemSolver<Sys>::Rescaler::Rescaler(boost::shared_ptr<Cluster> c, Mes& m) : cl
template<typename Sys>
void SystemSolver<Sys>::Rescaler::operator()() {
mes.Scaling = scaleClusters();
mes.Scaling = scaleClusters(calculateScale());
rescales++;
};
template<typename Sys>
typename SystemSolver<Sys>::Scalar SystemSolver<Sys>::Rescaler::scaleClusters() {
typename SystemSolver<Sys>::Scalar SystemSolver<Sys>::Rescaler::calculateScale() {
typedef typename Cluster::cluster_iterator citer;
std::pair<citer, citer> cit = cluster->clusters();
@ -134,6 +134,12 @@ typename SystemSolver<Sys>::Scalar SystemSolver<Sys>::Rescaler::scaleClusters()
sc = (s>sc) ? s : sc;
}
return sc;
}
template<typename Sys>
typename SystemSolver<Sys>::Scalar SystemSolver<Sys>::Rescaler::scaleClusters(Scalar sc) {
//if no scaling-value returned we can use 1
sc = (Kernel::isSame(sc,0, 1e-10)) ? 1. : sc;
@ -357,20 +363,33 @@ void SystemSolver<Sys>::solveCluster(boost::shared_ptr<Cluster> cluster, Sys& sy
//if(!has_cycle) {
#ifdef USE_LOGGING
BOOST_LOG(log)<< "non-cyclic system dedected: solve rotation only";
BOOST_LOG(log)<< "non-cyclic system dedected: solve rotation only";
#endif
//cool, lets do uncylic. first all rotational constraints with rotational parameters
mes.setAccess(rotation);
//cool, lets do uncylic. first all rotational constraints with rotational parameters
mes.setAccess(rotation);
//solve can be done without catching exceptions, because this only fails if the system is
//unsolvable. We need the rescaler here two as unscaled rotations may be unprecise if the
//parts are big (as ne the normals are always 1) and then distances may not be possible
//to solve correctly with only translations
Rescaler re(cluster, mes);
re();
//rotations need to be calculated in a scaled manner. thats because the normales used for
//rotation calculation are always 1, no matter how big the part is. This can lead to problems
//when for example two rotated faces have a precision error on the parallel normals but a distance
//at the outer edges is far bigger than the precision as the distance from normal origin to outer edge
//is bigger 1. that would lead to unsolvable translation-only systems.
//solve need to catch exceptions to reset the mes scaling on failure
Rescaler re(cluster, mes);
mes.Scaling = 1./(re.calculateScale()*SKALEFAKTOR);
try {
sys.kernel().solve(mes, re);
mes.Scaling = 1.;
}
catch(...) {
mes.Scaling = 1.;
throw;
}
//now let's see if we have to go on with the translations
if(mes.hasAccessType(general)) {
//now let's see if we have to go on with the translations
mes.setAccess(general);
mes.recalculate();
@ -392,6 +411,8 @@ void SystemSolver<Sys>::solveCluster(boost::shared_ptr<Cluster> cluster, Sys& sy
done = false;
}
}
};
//};
//not done already? try it the hard way!

3
src/Mod/Assembly/App/opendcm/module3d/solver.hpp
View File

@ -84,7 +84,8 @@ struct SystemSolver : public Job<Sys> {
void operator()();
Scalar scaleClusters();
Scalar calculateScale();
Scalar scaleClusters(Scalar sc);
void collectPseudoPoints(boost::shared_ptr<Cluster> parent,
LocalVertex cluster,
std::vector<typename Kernel::Vector3,