Commit fef6d895 authored by Lionel Untereiner's avatar Lionel Untereiner
Browse files

Merge cgogn.u-strasbg.fr:~kraemer/CGoGN

parents fc45ade8 601d62f3
......@@ -17,6 +17,7 @@ find_package(GLEW REQUIRED)
find_package(Qt4 REQUIRED)
find_package(QGLViewer REQUIRED)
find_package(PythonLibs 2.7 REQUIRED)
find_package(SuiteSparse REQUIRED)
SET( QT_USE_QTOPENGL TRUE )
SET( QT_USE_QTXML TRUE )
......@@ -66,6 +67,7 @@ SET (EXT_INCLUDES
${QT_INCLUDE_DIR}
${QGLVIEWER_INCLUDE_DIR}
${PYTHON_INCLUDE_DIRS}
${SUITESPARSE_INCLUDE_DIRS}
)
# define libs for external libs
......@@ -81,6 +83,7 @@ SET (EXT_LIBS
${QT_LIBRARIES}
${QGLVIEWER_LIBRARIES}
${PYTHON_LIBRARIES}
${SUITESPARSE_LIBRARIES}
)
......
......@@ -6,12 +6,17 @@
#include "ui_surfaceDeformation.h"
#include "Utils/drawer.h"
#include "Container/fakeAttribute.h"
#include "OpenNL/linear_solver.h"
#include "Algo/LinearSolving/basic.h"
using namespace CGoGN;
using namespace SCHNApps;
// namespace CGoGN { namespace Utils { class Drawer; } }
enum SelectionMode
{
......@@ -20,6 +25,8 @@ enum SelectionMode
};
typedef NoNameIOAttribute<Eigen::Matrix3f> Eigen_Matrix3f ;
struct PerMapParameterSet
{
PerMapParameterSet() {}
......@@ -32,6 +39,19 @@ struct PerMapParameterSet
SelectionMode verticesSelectionMode;
std::vector<unsigned int> locked_vertices;
std::vector<unsigned int> handle_vertices;
VertexAttribute<PFP2::VEC3> positionInit;
VertexAttribute<PFP2::VEC3> vertexNormal;
EdgeAttribute<PFP2::REAL> edgeAngle;
EdgeAttribute<PFP2::REAL> edgeWeight;
VertexAttribute<PFP2::REAL> vertexArea;
VertexAttribute<PFP2::VEC3> diffCoord;
VertexAttribute<Eigen_Matrix3f> vertexRotationMatrix;
VertexAttribute<PFP2::VEC3> rotatedDiffCoord;
VertexAttribute<unsigned int> vIndex;
unsigned int nb_vertices;
LinearSolver<PFP2::REAL>* solver;
};
struct ParameterSet
......@@ -116,11 +136,14 @@ public slots:
void cb_selectLockedVertices(bool b);
void cb_selectHandleVertices(bool b);
void matchDiffCoord(View* view, MapHandlerGen* map);
void asRigidAsPossible(View* view, MapHandlerGen* map);
private:
SurfaceDeformationDockTab* m_dockTab;
QHash<View*, ParameterSet*> h_viewParams;
Utils::Drawer* m_drawer;
// Utils::Drawer* m_drawer;
bool b_refreshingUI;
......
......@@ -3,6 +3,11 @@
#include "Algo/Selection/raySelector.h"
#include "Algo/Selection/collector.h"
#include "Algo/Geometry/normal.h"
#include "Algo/Geometry/laplacian.h"
//#include "Utils/drawer.h"
#include <QKeyEvent>
#include <QMouseEvent>
......@@ -10,10 +15,7 @@
PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
verticesSelectionMode(LOCKED)
{
GenericMap* map = mh->getGenericMap();
lockingMarker = new CellMarker<VERTEX>(*map);
handleMarker = new CellMarker<VERTEX>(*map);
PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
AttributeContainer& cont = map->getAttributeContainer<VERTEX>();
......@@ -42,6 +44,34 @@ PerMapParameterSet::PerMapParameterSet(MapHandlerGen* mh) :
}
}
}
lockingMarker = new CellMarker<VERTEX>(*map);
handleMarker = new CellMarker<VERTEX>(*map);
positionInit = mh->addAttribute<PFP2::VEC3, VERTEX>("positionInit") ;
map->copyAttribute(positionInit, positionAttribute) ;
vertexNormal = mh->addAttribute<PFP2::VEC3, VERTEX>("vertexNormal");
edgeAngle = mh->addAttribute<PFP2::REAL, EDGE>("edgeAngle");
edgeWeight = mh->addAttribute<PFP2::REAL, EDGE>("edgeWeight");
vertexArea = mh->addAttribute<PFP2::REAL, VERTEX>("vertexArea");
diffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("diffCoord");
vertexRotationMatrix = mh->addAttribute<Eigen_Matrix3f, VERTEX>("vertexRotationMatrix") ;
rotatedDiffCoord = mh->addAttribute<PFP2::VEC3, VERTEX>("rotatedDiffCoord") ;
Algo::Surface::Geometry::computeNormalVertices<PFP2>(*map, positionAttribute, vertexNormal) ;
Algo::Surface::Geometry::computeAnglesBetweenNormalsOnEdges<PFP2>(*map, positionAttribute, edgeAngle) ;
Algo::Surface::Geometry::computeCotanWeightEdges<PFP2>(*map, positionAttribute, edgeWeight) ;
Algo::Surface::Geometry::computeVoronoiAreaVertices<PFP2>(*map, positionAttribute, vertexArea) ;
Algo::Surface::Geometry::computeLaplacianCotanVertices<PFP2>(*map, edgeWeight, vertexArea, positionAttribute, diffCoord) ;
for(unsigned int i = vertexRotationMatrix.begin(); i != vertexRotationMatrix.end() ; vertexRotationMatrix.next(i))
vertexRotationMatrix[i] = Eigen::Matrix3f::Identity();
vIndex = mh->addAttribute<unsigned int, VERTEX>("vIndex");
nb_vertices = static_cast<AttribMap*>(map)->computeIndexCells<VERTEX>(vIndex);
solver = new LinearSolver<PFP2::REAL>(nb_vertices);
}
PerMapParameterSet::~PerMapParameterSet()
......@@ -65,25 +95,22 @@ bool SurfaceDeformationPlugin::enable()
connect(m_window, SIGNAL(viewAndPluginUnlinked(View*, Plugin*)), this, SLOT(viewUnlinked(View*, Plugin*)));
connect(m_window, SIGNAL(currentViewChanged(View*)), this, SLOT(currentViewChanged(View*)));
Utils::ShaderColorPerVertex* s = new Utils::ShaderColorPerVertex();
registerShader(s);
m_drawer = new Utils::Drawer();
registerShader(m_drawer->getShader());
// m_drawer = new Utils::Drawer();
// registerShader(m_drawer->getShader());
return true;
}
void SurfaceDeformationPlugin::disable()
{
delete m_drawer;
// delete m_drawer;
}
void SurfaceDeformationPlugin::redraw(View* view)
{
if(selecting)
{
glDisable(GL_LIGHTING) ;
/* glDisable(GL_LIGHTING) ;
m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
m_drawer->lineWidth(2.0f) ;
m_drawer->begin(GL_LINES) ;
......@@ -102,7 +129,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
m_drawer->vertex(selectionCenter + selectionRadius * PFP2::VEC3(0,0,-1)) ;
m_drawer->end() ;
m_drawer->endList() ;
}
*/ }
ParameterSet* params = h_viewParams[view];
MapHandlerGen* mh = params->selectedMap;
......@@ -112,7 +139,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
if(!perMap->locked_vertices.empty() || !perMap->handle_vertices.empty())
{
glDisable(GL_LIGHTING) ;
/* glDisable(GL_LIGHTING) ;
m_drawer->newList(GL_COMPILE_AND_EXECUTE) ;
m_drawer->pointSize(4.0f) ;
m_drawer->begin(GL_POINTS) ;
......@@ -127,7 +154,7 @@ void SurfaceDeformationPlugin::redraw(View* view)
m_drawer->vertex(perMap->positionAttribute[perMap->handle_vertices[i]]) ;
m_drawer->end() ;
m_drawer->endList() ;
}
*/ }
}
}
......@@ -176,35 +203,26 @@ void SurfaceDeformationPlugin::mousePress(View* view, QMouseEvent* event)
neigh.collectAll(d) ;
const std::vector<Dart>& insideV = neigh.getInsideVertices() ;
if(perMap->verticesSelectionMode == LOCKED)
for(unsigned int i = 0; i < insideV.size(); ++i)
{
for(unsigned int i = 0; i < insideV.size(); ++i)
unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
if (!perMap->lockingMarker->isMarked(v))
{
unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
if (!perMap->lockingMarker->isMarked(v))
{
perMap->locked_vertices.push_back(v) ;
perMap->lockingMarker->mark(v);
}
perMap->locked_vertices.push_back(v) ;
perMap->lockingMarker->mark(v);
}
}
else
{
for(unsigned int i = 0; i < insideV.size(); ++i)
if(perMap->verticesSelectionMode == HANDLE)
{
unsigned int v = map->getEmbedding<VERTEX>(insideV[i]) ;
if(!perMap->handleMarker->isMarked(v))
{
perMap->handle_vertices.push_back(v) ;
perMap->handleMarker->mark(v);
}
if(!perMap->lockingMarker->isMarked(v))
{
perMap->locked_vertices.push_back(v) ;
perMap->lockingMarker->mark(v) ;
}
}
}
LinearSolving::resetSolver(perMap->solver, false) ;
view->updateGL() ;
}
}
......@@ -243,7 +261,8 @@ void SurfaceDeformationPlugin::mouseRelease(View* view, QMouseEvent* event)
void SurfaceDeformationPlugin::mouseMove(View* view, QMouseEvent* event)
{
ParameterSet* params = h_viewParams[view];
PerMapParameterSet* perMap = params->getCurrentMapParameterSet();
MapHandlerGen* map = params->selectedMap;
PerMapParameterSet* perMap = params->perMap[map->getName()];
if(dragging)
{
......@@ -257,9 +276,9 @@ void SurfaceDeformationPlugin::mouseMove(View* view, QMouseEvent* event)
dragPrevious = q ;
// matchDiffCoord() ;
// matchDiffCoord(view, map);
// for(unsigned int i = 0; i < 2; ++i)
// asRigidAsPossible();
asRigidAsPossible(view, map);
params->selectedMap->updateVBO(perMap->positionAttribute);
......@@ -470,6 +489,150 @@ void SurfaceDeformationPlugin::cb_selectHandleVertices(bool b)
}
}
void SurfaceDeformationPlugin::matchDiffCoord(View* view, MapHandlerGen* mh)
{
PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
PerMapParameterSet* perMap = h_viewParams[view]->perMap[mh->getName()];
LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
for(int coord = 0; coord < 3; ++coord)
{
LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
LinearSolving::startMatrix(perMap->solver) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
// LinearSolving::addRowsRHS_Laplacian_Topo<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->diffCoord, coord);
LinearSolving::endMatrix(perMap->solver) ;
LinearSolving::solve(perMap->solver) ;
LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
LinearSolving::resetSolver(perMap->solver, true) ;
}
}
void SurfaceDeformationPlugin::asRigidAsPossible(View* view, MapHandlerGen* mh)
{
PFP2::MAP* map = static_cast<MapHandler<PFP2>*>(mh)->getMap();
PerMapParameterSet* perMap = h_viewParams[view]->perMap[mh->getName()];
CellMarkerNoUnmark<VERTEX> m(*map) ;
for(Dart d = map->begin(); d != map->end(); map->next(d))
{
if(!m.isMarked(d))
{
m.mark(d) ;
Eigen::Matrix3f cov = Eigen::Matrix3f::Zero() ;
PFP2::VEC3 p = perMap->positionAttribute[d] ;
PFP2::VEC3 pInit = perMap->positionInit[d] ;
PFP2::REAL area = perMap->vertexArea[d] ;
Dart it = d ;
do
{
Dart neigh = map->phi1(it) ;
PFP2::VEC3 v = perMap->positionAttribute[neigh] - p ;
PFP2::VEC3 vv = perMap->positionInit[neigh] - pInit ;
for(unsigned int i = 0; i < 3; ++i)
for(unsigned int j = 0; j < 3; ++j)
cov(i,j) += v[i] * vv[j] * perMap->edgeWeight[it] / area ;
Dart dboundary = map->phi_1(it) ;
if(map->phi2(dboundary) == dboundary)
{
v = perMap->positionAttribute[dboundary] - p ;
vv = perMap->positionInit[dboundary] - p ;
for(unsigned int i = 0; i < 3; ++i)
for(unsigned int j = 0; j < 3; ++j)
cov(i,j) += v[i] * vv[j] * perMap->edgeWeight[dboundary] / area ;
}
it = map->alpha1(it) ;
} while(it != d) ;
Eigen::JacobiSVD<Eigen::Matrix3f> svd(cov, Eigen::ComputeFullU | Eigen::ComputeFullV) ;
Eigen::Matrix3f R = svd.matrixU() * svd.matrixV().transpose() ;
if(R.determinant() < 0)
{
Eigen::Matrix3f U = svd.matrixU() ;
for(unsigned int i = 0; i < 3; ++i)
U(i,2) *= -1 ;
R = U * svd.matrixV().transpose() ;
}
perMap->vertexRotationMatrix[d] = R ;
}
}
for(Dart d = map->begin(); d != map->end(); map->next(d))
{
if(m.isMarked(d))
{
m.unmark(d) ;
unsigned int degree = 0 ;
Eigen::Matrix3f r = Eigen::Matrix3f::Zero() ;
Dart it = d ;
do
{
r += perMap->vertexRotationMatrix[map->phi1(it)] ;
++degree ;
Dart dboundary = map->phi_1(it) ;
if(map->phi2(dboundary) == dboundary)
{
r += perMap->vertexRotationMatrix[dboundary] ;
++degree ;
}
it = map->alpha1(it) ;
} while(it != d) ;
r += perMap->vertexRotationMatrix[d] ;
r /= degree + 1 ;
PFP2::VEC3& dc = perMap->diffCoord[d] ;
Eigen::Vector3f rdc(dc[0], dc[1], dc[2]) ;
rdc = r * rdc ;
perMap->rotatedDiffCoord[d] = PFP2::VEC3(rdc[0], rdc[1], rdc[2]) ;
// Eigen::Vector3f rdc = Eigen::Vector3f::Zero() ;
// Dart it = d ;
// PFP::REAL vArea = perMap->vertexArea[d] ;
// Eigen::Matrix3f& rotM = perMap->vertexRotationMatrix[d] ;
// PFP::REAL val = 0 ;
// do
// {
// Dart ddn = map->phi1(it) ;
// PFP::REAL w = perMap->edgeWeight[it] / vArea ;
// PFP::VEC3 vv = (perMap->positionInit[ddn] - perMap->positionInit[it]) * w ;
// val += w ;
// Eigen::Matrix3f r = 0.5f * (perMap->vertexRotationMatrix[ddn] + rotM) ;
// Eigen::Vector3f vvr = r * Eigen::Vector3f(vv[0], vv[1], vv[2]) ;
// rdc += vvr ;
// Dart dboundary = map->phi_1(it) ;
// if(map->phi2(dboundary) == dboundary)
// {
// w = perMap->edgeWeight[dboundary] / vArea ;
// vv = (perMap->positionInit[dboundary] - perMap->positionInit[it]) * w ;
// val += w ;
// r = 0.5f * (perMap->vertexRotationMatrix[dboundary] + rotM) ;
// vvr = r * Eigen::Vector3f(vv[0], vv[1], vv[2]) ;
// rdc += vvr ;
// }
// it = map->alpha1(it) ;
// } while(it != d) ;
// rdc /= val ;
// perMap->rotatedDiffCoord[d] = PFP::VEC3(rdc[0], rdc[1], rdc[2]) ;
}
}
LinearSolving::initSolver(perMap->solver, perMap->nb_vertices, true, true) ;
for(int coord = 0; coord < 3; ++coord)
{
LinearSolving::setupVariables<PFP2>(*map, perMap->solver, perMap->vIndex, *perMap->lockingMarker, perMap->positionAttribute, coord) ;
LinearSolving::startMatrix(perMap->solver) ;
LinearSolving::addRowsRHS_Laplacian_Cotan<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->rotatedDiffCoord, coord) ;
LinearSolving::endMatrix(perMap->solver) ;
LinearSolving::solve(perMap->solver) ;
LinearSolving::getResult<PFP2>(*map, perMap->solver, perMap->vIndex, perMap->positionAttribute, coord) ;
LinearSolving::resetSolver(perMap->solver, true) ;
}
}
void SurfaceDeformationDockTab::refreshUI(ParameterSet* params)
......
......@@ -5,15 +5,15 @@ differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
surfaceDeformationPlugin = schnapps.loadPlugin("SurfaceDeformation");
obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/midRes/egea_remesh_9k.off");
#obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
v = schnapps.getView("view_0");
#v = schnapps.getView("view_0");
schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
#schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
#schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
#schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
schnapps.linkViewAndMap(v.getName(), obj.getName());
#schnapps.linkViewAndMap(v.getName(), obj.getName());
#differentialPropertiesPlugin.computeNormal(obj.getName());
#differentialPropertiesPlugin.computeCurvature(obj.getName());
#include <QSplashScreen>
#include "window.h"
#include <QFileInfo>
#include "PythonQt/PythonQt.h"
#include "PythonQt/gui/PythonQtScriptingConsole.h"
int main(int argc, char* argv[])
{
QApplication app(argc, argv);
......
......@@ -32,17 +32,12 @@
#define __CHOLESKY_SOLVER__
#include <cassert>
#include "OpenNL/sparse_matrix.h"
#include "CHOLMOD/cholmod.h"
#include <cholmod.h>
template< class MATRIX, class VECTOR >
template <typename CoeffType>
class Solver_CHOLESKY
{
public:
typedef MATRIX Matrix ;
typedef VECTOR Vector ;
typedef typename Vector::CoeffType CoeffType ;
Solver_CHOLESKY() {
N = 0 ;
factorized_ = false ;
......@@ -58,14 +53,14 @@ public:
bool factorized() { return factorized_ ; }
void factorize(const MATRIX &A_in) {
assert(A_in.n() == A_in.m()) ;
assert(A_in.has_symmetric_storage()) ;
void factorize(const Eigen::SparseMatrix<CoeffType>& A_in) {
assert(A_in.rows() == A_in.cols()) ;
// assert(A_in.has_symmetric_storage()) ;
N = A_in.n() ;
int NNZ = A_in.nnz() ;
N = A_in.rows() ;
// Translate sparse matrix into cholmod format
int NNZ = A_in.nonZeros() ;
cholmod_sparse* A = cholmod_allocate_sparse(N, N, NNZ, false, true, -1, CHOLMOD_REAL, &c);
int* colptr = static_cast<int*>(A->p) ;
......@@ -75,12 +70,12 @@ public:
// Convert local Matrix into CHOLMOD Matrix
int count = 0 ;
for(int j = 0; j < N; j++) {
const typename SparseMatrix<CoeffType>::Column& Cj = A_in.column(j) ;
colptr[j] = count ;
for(unsigned int ii = 0; ii < Cj.nb_coeffs(); ii++) {
a[count] = Cj.coeff(ii).a ;
rowind[count] = Cj.coeff(ii).index ;
count++ ;
colptr[j] = count;
for(typename Eigen::SparseMatrix<CoeffType>::InnerIterator it(A_in, j); it; ++it)
{
a[count] = it.value();
rowind[count] = it.row();
++count;
}
}
colptr[N] = NNZ ;
......@@ -94,10 +89,10 @@ public:
cholmod_free_sparse(&A, &c) ;
}
void solve(const VECTOR& b_in, VECTOR& x_out) {
void solve(const Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& b_in, Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>& x_out) {
assert(factorized_) ;
assert(L->n == b_in.dimension()) ;
assert(L->n == x_out.dimension()) ;
assert(L->n == b_in.rows()) ;
assert(L->n == x_out.rows()) ;
// Translate right-hand side into cholmod format
cholmod_dense* b = cholmod_allocate_dense(N, 1, N, CHOLMOD_REAL, &c) ;
......@@ -114,8 +109,8 @@ public:
x_out[i] = cx[i] ;
// Cleanup
cholmod_free_dense(&x, &c) ;
cholmod_free_dense(&b, &c) ;
cholmod_free_dense(&x, &c) ;
}
void reset() {
......
......@@ -37,6 +37,7 @@
#include <cstdlib>
#include <Eigen/Sparse>
#include "cholesky_solver.h"
template <class CoeffType>
......@@ -80,7 +81,7 @@ public:
delete A_ ;
delete x_ ;
delete b_ ;
// delete direct_solver_ ;
delete direct_solver_ ;
}
// __________________ Parameters ________________________
......@@ -188,9 +189,7 @@ private:
Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* x_ ;
Eigen::Matrix<CoeffType, Eigen::Dynamic, 1>* b_ ;
Eigen::ConjugateGradient<Eigen::SparseMatrix<CoeffType> >* symmetric_solver_;
Eigen::BiCGSTAB<Eigen::SparseMatrix<CoeffType> >* nonsymmetric_solver_;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<CoeffType> >* direct_solver_;
Solver_CHOLESKY<CoeffType>* direct_solver_ ;
} ;
#include "OpenNL/linear_solver.hpp"
......
......@@ -33,9 +33,7 @@ LinearSolver<CoeffType>::LinearSolver(unsigned int nb_variables) {
A_ = NULL ;
x_ = NULL ;
b_ = NULL ;
symmetric_solver_ = NULL;
nonsymmetric_solver_ = NULL;
direct_solver_ = NULL;
direct_solver_ = new Solver_CHOLESKY<CoeffType>() ;
}
template <typename CoeffType>
......@@ -138,7 +136,7 @@ void LinearSolver<CoeffType>::end_row() {
if(!matrix_already_set_) {
for(unsigned int i = 0; i < nf; i++) {
for(unsigned int j = 0; j < nf; j++) {
A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j];
A_->coeffRef(if_[i], if_[j]) += af_[i] * af_[j] ;
}
}
}
......@@ -156,11 +154,11 @@ void LinearSolver<CoeffType>::end_row() {
// Construct the matrix coefficients of the current row
if(!matrix_already_set_) {
for(unsigned int i = 0; i < nf; i++) {
A_->coeffRef(current_row_, if_[i]) += af_[i];
A_->coeffRef(current_row_, if_[i]) += af_[i] ;
}
}