/*
 * \file mesh.cpp
 * \author Angus Leeming
 */

#include <boost/numeric/ublas/config.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <boost/random.hpp>
//#include <boost/nondet_random.hpp>

#include <iostream>
#include <vector>


namespace {

typedef std::size_t size_type;

int get_random(int n);
void perturb_mesh(std::vector<double> & x, double total_length);

} // namespace anon


int main(int argc, char * argv[])
{
	int const nsprings = 10;

	// Initial length of all springs is 1.
	// x is the position of the "nodes" at the end of each spring.
	std::vector<double> x(nsprings);
	for (std::vector<double>::size_type i = 0; i != nsprings; ++i)
		x[i] = double(i);

	perturb_mesh(x, x[nsprings-1]);

	std::cout << "Node[0] position: " << x[0] << '\n';;
	for (std::vector<double>::size_type i = 1; i != nsprings; ++i)
		std::cout << "Node[" << i << "] position: " << x[i]
			  << ", element length: " << (x[i] - x[i-1])
			  << '\n';
	
	return 0;
}


namespace {

void perturb_mesh(std::vector<double> & x,
		  double total_length)
{
	// Use a mass-spring analogy to perturb the mesh.
	//
	//                 k_{n-1}           k_{n}
	//   F_{n-1} <-----------------o-----------------> F_{n}
	//                \el_{n-1}    n    \el_{n}
	//
	// where F is the force applied,
	//       k is the spring stiffness and
	//       \el is the element lenght.
	//
	// At equilibrium  F_{n-1} - F_{n} = 0
	// for each node n, 1 <= n < nnodes-1
	// where F_{n} = k_{n} \times \el_{n}
	//
	// Together with  \el_{total} = \sum_{n}\el_{n}
	// this allows us to formulate the system
	//
	// |k_0  -k_1            ||\el_0|   |0          |
	// |      k_1  -k_2      ||\el_1| = |0          |
	// |            k_2  -k_3||\el_2|   |0          |
	// |  1     1     1     1||\el_3|   |\el_{total}|
	//
	// where |...| indicate a matrix or vector rather than a determinant,
	//
	// and then solve it to give the unknown element lengths, \el_{n}.

	using namespace boost::numeric;

	// There are one fewer springs than there are nodes.
	size_type system_size = x.size() - 1;
	ublas::scalar_vector<double> const zero(system_size, 0.0);

	// Assign the spring stiffnesses in the range 0.8 <= k <= 1.2
	ublas::vector<double> k(zero);
	ublas::vector<double>::iterator kit = k.begin();
	ublas::vector<double>::iterator const kend = k.end();
	for (; kit != kend; ++kit)
		*kit = 0.8 + (1.0E-04 * get_random(4000));

	// Populate the matrix A.
	ublas::sparse_matrix<double> A(system_size, system_size);
	for (size_type i = 0; i != A.size1() - 1; ++i) {
		A(i,i)   = k(i);
		A(i,i+1) = -k(i+1);
	}
	for (size_type i = 0; i != A.size2(); ++i)
		A(A.size1() - 1, i) = 1;

	// On return from lu_factorize,
	// A contains the LU factorization and pm the row pivots.
	ublas::permutation_matrix<size_type> pm(A.size1());
	lu_factorize(A, pm);

	// Initialise el so that it contains the 'b' of Ax = b
	ublas::vector<double> el(zero);
	el(el.size() - 1) = total_length;

        // On return from lu_substitute,
	// el contains the 'x' of Ax=b.
        lu_substitute(A, pm, el);

	// Ok, we have a collection of element lengths.
	// Use them to ascertain the position of each node.
	std::vector<double>::iterator xit = x.begin();
	std::vector<double>::iterator const xend = x.end();
	ublas::vector<double>::const_iterator elit = el.begin();

	x[0] = 0;
	double length = 0;
	++xit;
	for (; xit != xend; ++xit, ++elit) {
		length += *elit;
		*xit = length;
	}
}


// Change this type to get different random number generators.
// Note that boost::random_device is *much* slower than boost::mt19937
// (Run times are three times longer.)
typedef boost::mt19937 rng_t;
//typedef boost::random_device rng_t;

int get_random(int n)
{
	// Use "rng_t &" here because boost::random_device is noncopyable.
	typedef boost::variate_generator<rng_t &, boost::uniform_int<> >
		variate_generator;
	// Produces randomness out of thin air.
	static rng_t rng;
	// Distribution that maps to 1..6. Arbitrary choice...
	static boost::uniform_int<> mapping(1,6);
	// Glues randomness with mapping.
	static variate_generator die(rng, mapping);

	// Reset mapping to the real, desired range.
	die.distribution() = boost::uniform_int<>(0,n-1);
	// Roll the die and generate a random number.
	return die();
}

} // namespace anon