Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/test/unit/test_shape_matrix.cc
4a52f2a5948cmaster

Last Change
Raw File

test_shape_matrix.cc
View Options

#include "shape_matrix.hh"
#include <aka_iterators.hh>
#include <gtest/gtest.h>
/* -------------------------------------------------------------------------- */
using namespace libmultiscale;
/* -------------------------------------------------------------------------- */
struct Mesh {
Mesh(Real L, UInt n_el);
void print();
Eigen::Vector4d getShape(Eigen::Vector2d);
Eigen::Vector2i getElement(Eigen::Vector2d);
Eigen::Vector4i getNodeIndexes(Eigen::Vector2i);
int glob_idx(int i, int j) { return i * n_nd + j; };
Real L;
UInt n_el;
UInt n_nd;
UInt n_elements;
UInt n_nodes;
Eigen::ArrayXXd elements;
Eigen::ArrayXXd nodes;
};
Mesh::Mesh(Real L, UInt n_el)
: L(L), n_el(n_el), n_nd(n_el + 1), n_elements(n_el * n_el),
n_nodes((n_el + 1) * (n_el + 1)), elements(n_elements, 4),
nodes(n_nodes, 2) {
for (UInt i = 0; i < n_el; ++i) {
for (UInt j = 0; j < n_el; ++j) {
UInt el = i * n_el + j;
auto nd_indexes = getNodeIndexes({i, j});
Real dL = L / n_el;
Eigen::ArrayXXd nd_pos(4, 2);
nd_pos.row(0) << i * dL, j * dL;
nd_pos.row(1) << i * dL, (j + 1) * dL;
nd_pos.row(2) << (i + 1) * dL, (j + 1) * dL;
nd_pos.row(3) << (i + 1) * dL, j * dL;
for (UInt n = 0; n < 4; ++n) {
elements(el, n) = nd_indexes[n];
nodes.row(nd_indexes[n]) = nd_pos.row(n);
}
}
}
}
void Mesh::print() {
Eigen::IOFormat NumpyFmt(Eigen::FullPrecision, 0, ", ", ",\n", "[", "]", "[",
"]");
std::cout << this->elements.format(NumpyFmt) << std::endl;
std::cout << this->nodes.format(NumpyFmt) << std::endl;
}
Eigen::Vector2i Mesh::getElement(Eigen::Vector2d X) {
return (n_el * X / L).cast<int>();
}
Eigen::Vector4i Mesh::getNodeIndexes(Eigen::Vector2i I) {
auto i = I[0];
auto j = I[1];
return Eigen::Vector4i{glob_idx(i, j), glob_idx(i, j + 1),
glob_idx(i + 1, j + 1), glob_idx(i + 1, j)};
}
Eigen::Vector4d Mesh::getShape(Eigen::Vector2d X) {
auto I = getElement(X);
// std::cout << "I " << I << "****" << std::endl;
if (I[0] < 0 or I[1] < 0 or I[0] > int(n_el) or I[1] > int(n_el)) {
throw std::runtime_error("point out of domain");
}
Eigen::Vector2d c = X / L - I.cast<Real>() * L / n_el;
c = c * 2. - Eigen::Vector2d::Constant(1.);
// std::cout << "c " << c << "****" << std::endl;
Eigen::Vector4d N;
N(0) = 1. / 4. * (1. - c(0)) * (1. - c(1)); /// N1(q_0)
N(1) = 1. / 4. * (1. + c(0)) * (1. - c(1)); /// N2(q_0)
N(2) = 1. / 4. * (1. + c(0)) * (1. + c(1)); /// N3(q_0)
N(3) = 1. / 4. * (1. - c(0)) * (1. + c(1)); /// N4(q_0)
return N;
}
////////////////////////////////////////////////////////////////
TEST(GRID, creation) {
auto mesh = Mesh(1., 10);
std::cout << mesh.getShape({0.5, 0.5}) << std::endl;
std::cout << mesh.getShape({0.501, 0.501}) << std::endl;
std::cout << mesh.getShape({0.49999, 0.49999}) << std::endl;
}
ShapeMatrix fill_shp_matrix(Mesh &mesh, Eigen::ArrayX2d &Xs) {
UInt n_points = Xs.rows();
UInt n_nodes = mesh.n_nodes;
std::cout << n_nodes << " " << n_points << std::endl;
ShapeMatrix shp(n_nodes, n_points);
std::cout << shp.rows() << " " << shp.cols() << std::endl;
for (int i = 0; i < Xs.rows(); ++i) {
Eigen::Vector2d X = Xs.row(i);
auto el = mesh.getElement(X);
auto node_indexes = mesh.getNodeIndexes(el);
auto &&shapes = mesh.getShape(X);
for (int j = 0; j < 4; ++j) {
auto &&I = node_indexes[j];
shp.insert(I, i) = shapes(j);
}
}
return shp;
}
TEST(SHAPE_MATRIX, filling) {
auto mesh = Mesh(1., 10);
Eigen::ArrayX2d Xs = Eigen::ArrayX2d::Random(12, 2);
Xs = (Xs + 1) / 2;
auto shp = fill_shp_matrix(mesh, Xs);
EXPECT_EQ(shp.cols(), 12);
EXPECT_EQ(shp.rows(), 121);
}
template <typename ViewType, typename Derived> struct make_view_impl {
make_view_impl(Eigen::DenseBase<Derived> &dense_obj) : dense_obj(dense_obj){};
struct iterator {
iterator(typename Derived::Scalar *ptr);
bool operator!=(iterator &other);
void operator++();
Eigen::Map<ViewType> operator*();
};
iterator begin() { return iterator(dense_obj.data()); };
iterator end();
Eigen::DenseBase<Derived> &dense_obj;
};
template <typename ViewType, typename Derived>
auto make_view(Eigen::DenseBase<Derived> &dense_obj) {
return make_view_impl<ViewType, Derived>(dense_obj);
}
TEST(SHAPE_MATRIX, interpolation_constant_field) {
auto mesh = Mesh(1., 99);
Eigen::ArrayX2d Xs = Eigen::ArrayX2d::Random(100, 2);
Xs = (Xs + 1) / 2;
ShapeMatrix shp = fill_shp_matrix(mesh, Xs);
EXPECT_EQ(shp.cols(), 100);
EXPECT_EQ(shp.rows(), 100 * 100);
Eigen::ArrayX2d constant_node_field(mesh.n_nodes, 2);
constant_node_field.col(0) = 1.;
constant_node_field.col(1) = 2.;
std::cout << shp.rows() << " " << shp.cols() << std::endl;
Eigen::ArrayX2d interpolated = shp.transpose() * constant_node_field.matrix();
for (auto &&v : interpolated.rowwise()) {
EXPECT_NEAR(v[0], 1., 1e-10);
EXPECT_NEAR(v[1], 2., 1e-10);
}
}
TEST(AKA_ITERATORS, zip) {
std::vector v1{1, 2, 3, 4};
std::vector v2{5, 6, 7, 8};
UInt cpt = 0;
for (auto &&[a, b] : zip(v1, v2)) {
EXPECT_EQ(a, v1[cpt]);
EXPECT_EQ(b, v2[cpt]);
++cpt;
}
}
c4science · Help