/* * \file mesh.cpp * \author Angus Leeming */ #include #include #include #include #include #include #include #include //#include #include #include namespace { typedef std::size_t size_type; int get_random(int n); void perturb_mesh(std::vector & 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 x(nsprings); for (std::vector::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::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 & 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 const zero(system_size, 0.0); // Assign the spring stiffnesses in the range 0.8 <= k <= 1.2 ublas::vector k(zero); ublas::vector::iterator kit = k.begin(); ublas::vector::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 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 pm(A.size1()); lu_factorize(A, pm); // Initialise el so that it contains the 'b' of Ax = b ublas::vector 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::iterator xit = x.begin(); std::vector::iterator const xend = x.end(); ublas::vector::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 > 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