diff --git a/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp b/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp
index 9e9f8085d..cd2ae8d5a 100644
--- a/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp
+++ b/src/Mod/Assembly/App/opendcm/core/imp/kernel_imp.hpp
@@ -21,6 +21,7 @@
#define DCM_KERNEL_IMP_H
#include "../kernel.hpp"
+#include
#include
#include
@@ -62,9 +63,7 @@ void Dogleg::calculateStep(const Eigen::MatrixBase& 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: "<::calculateStep(const Eigen::MatrixBase& 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: "<= delta) {
//h_dl = alpha*h_sd;
@@ -94,8 +90,6 @@ void Dogleg::calculateStep(const Eigen::MatrixBase& g, const Ei
BOOST_LOG(log)<< "Unnormal dogleg descent detected: "<::calculateStep(const Eigen::MatrixBase& g, const Ei
BOOST_LOG(log)<< "Unnormal dogleg beta detected: "<::MappedEquationSystem::MappedEquationSystem(int params
m_parameter(params), m_residual(equations),
m_params(params), m_eqns(equations), Scaling(1.),
Jacobi(&m_jacobi(0,0),equations,params,DynStride(equations,1)),
- Parameter(&m_parameter(0),params,DynStride(1,1)),
+ Parameter(&m_parameter(0),params,DynStride(1,1)),
Residual(&m_residual(0),equations,DynStride(1,1)) {
m_param_rot_offset = 0;
@@ -421,21 +413,22 @@ void Kernel::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 class Nonlinear>
bool Kernel::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_offset0);
+ return (m_params>0 && m_eqns>0);
};
template class Nonlinear>
diff --git a/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp b/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp
index b4b843c43..77e6e0fc7 100644
--- a/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp
+++ b/src/Mod/Assembly/App/opendcm/module3d/imp/solver_imp.hpp
@@ -107,12 +107,12 @@ SystemSolver::Rescaler::Rescaler(boost::shared_ptr c, Mes& m) : cl
template
void SystemSolver::Rescaler::operator()() {
- mes.Scaling = scaleClusters();
+ mes.Scaling = scaleClusters(calculateScale());
rescales++;
};
template
-typename SystemSolver::Scalar SystemSolver::Rescaler::scaleClusters() {
+typename SystemSolver::Scalar SystemSolver::Rescaler::calculateScale() {
typedef typename Cluster::cluster_iterator citer;
std::pair cit = cluster->clusters();
@@ -134,6 +134,12 @@ typename SystemSolver::Scalar SystemSolver::Rescaler::scaleClusters()
sc = (s>sc) ? s : sc;
}
+ return sc;
+}
+
+template
+typename SystemSolver::Scalar SystemSolver::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::solveCluster(boost::shared_ptr 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::solveCluster(boost::shared_ptr cluster, Sys& sy
done = false;
}
}
+ };
+
//};
//not done already? try it the hard way!
diff --git a/src/Mod/Assembly/App/opendcm/module3d/solver.hpp b/src/Mod/Assembly/App/opendcm/module3d/solver.hpp
index 0379b7d14..45e750ac9 100644
--- a/src/Mod/Assembly/App/opendcm/module3d/solver.hpp
+++ b/src/Mod/Assembly/App/opendcm/module3d/solver.hpp
@@ -84,7 +84,8 @@ struct SystemSolver : public Job {
void operator()();
- Scalar scaleClusters();
+ Scalar calculateScale();
+ Scalar scaleClusters(Scalar sc);
void collectPseudoPoints(boost::shared_ptr parent,
LocalVertex cluster,
std::vector