taucssolver.cpp

finite element mesh 参数化有限元网格划分

//=============================================================================
//   Simplified version of:
//   Example code for the full-day course
//
//   M. Botsch, M. Pauly, L. Kobbelt, P. Alliez, B. Levy,
//   "Geometric Modeling Based on Polygonal Meshes"
//   held at SIGGRAPH 2007, San Diego, and Eurographics 2008, Crete.
//
//   Copyright (C) 2007 by  Computer Graphics Laboratory, ETH Zurich, 
//                      and Computer Graphics Group,      RWTH Aachen
//
//                                                                            
//-----------------------------------------------------------------------------
//                                                                            
//                                License                                     
//                                                                            
//   This program is free software; you can redistribute it and/or
//   modify it under the terms of the GNU General Public License
//   as published by the Free Software Foundation; either version 2
//   of the License, or (at your option) any later version.
//   
//   This program is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//   GNU General Public License for more details.
//   
//   You should have received a copy of the GNU General Public License
//   along with this program; if not, write to the Free Software
//   Foundation, Inc., 51 Franklin Street, Fifth Floor, 
//   Boston, MA  02110-1301, USA.
//                                                                            
//=============================================================================
//=============================================================================
//
//  CLASS TaucsSolver - IMPLEMENTATION
//
//=============================================================================


//== INCLUDES =================================================================
#include "stdafx.h"
#include "TaucsSolver.h"
#include <iostream>


//== IMPLEMENTATION ==========================================================


//-----------------------------------------------------------------------------


TaucsSolver::
TaucsSolver()
: F(NULL), n_rows(0) 
{
}


//-----------------------------------------------------------------------------


TaucsSolver::
~TaucsSolver()
{
	delete_matrices();
}


//-----------------------------------------------------------------------------


void
TaucsSolver::
delete_matrices()
{
	if (F) {
		taucs_linsolve(NULL, &F, 0, NULL, NULL, NULL, NULL); 
		F = NULL;
	}
}


//-----------------------------------------------------------------------------


void
TaucsSolver::
begin_row()
{
	if (colptr.empty() || colptr.back() != (int)values.size())
	{
		colptr.push_back(values.size());
		n_rows = colptr.size()-1;
	}
}


void
TaucsSolver::
add_value(int _i, double _val)
{
	if (_i <= n_rows)
	{
		values.push_back(_val);
		rowind.push_back(_i);
	}
}


void
TaucsSolver::
end_row()
{
	if (colptr.empty() || colptr.back() != (int)values.size())
	{
		colptr.push_back(values.size());
		n_rows = colptr.size()-1;
	}
}


//-----------------------------------------------------------------------------


bool
TaucsSolver::
factorize()
{
	// delete old matrices
	delete_matrices();


	// setup ccs matrix
	A.n        = colptr.size()-1;
	A.m        = colptr.size()-1;
	A.flags    = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
	A.colptr   = &colptr[0];
	A.rowind   = &rowind[0];
	A.values.d = &values[0];

	// Factor matrix
	char* options[] = {"taucs.factor.LLT=true", NULL};
	int rc = taucs_linsolve(&A, &F, 0, NULL, NULL, options, NULL);
	if(rc != TAUCS_SUCCESS)
      return false;	

	return true;
}

//-----------------------------------------------------------------------------

bool
TaucsSolver::
solve(std::vector<double>& _b, std::vector<double>& _x)
{
	const unsigned int N = A.n;

	if (N != _b.size() || N != _x.size())
	{
		std::cerr << "TaucsSolver: matrix size doesn't match vector size\n";
		return false;
	}

	char* options[] = {"taucs.factor=false", NULL};
    int rc = taucs_linsolve(&A, &F, 1, &_x[0], &_b[0], options, NULL);
	if(rc != TAUCS_SUCCESS)
      return false;

	return true;
}