Page MenuHomec4science
Files F92359970

test_shape_matrix.cc
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 19, 16:41

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);
Eigen::Vector2d naturalCoordinates(Eigen::Vector2i, Eigen::Vector2d);
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 + 1) * dL, j * dL;
nd_pos.row(2) << (i + 1) * dL, (j + 1) * dL;
nd_pos.row(3) << i * dL, (j + 1) * 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 + 1, j),
glob_idx(i + 1, j + 1), glob_idx(i, j + 1)};
}
Eigen::Vector2d Mesh::naturalCoordinates(Eigen::Vector2i Ielem,
Eigen::Vector2d X) {
if (Ielem[0] < 0 or Ielem[1] < 0 or Ielem[0] > int(n_el) or
Ielem[1] > int(n_el)) {
throw std::runtime_error("element out of grid");
}
Real dL = L / n_el;
Eigen::Array2d c = X / dL - Ielem.cast<Real>();
c = c * 2. - 1.;
return c;
}
Eigen::Vector4d Mesh::getShape(Eigen::Vector2d X) {
auto I = getElement(X);
auto c = naturalCoordinates(I, X);
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, search_element) {
auto mesh = Mesh(1., 2);
auto test = [&mesh](Eigen::Vector2d pos, Eigen::Vector2d res) {
auto c = mesh.getElement(pos);
EXPECT_NEAR(res[0], c[0], 1e-10);
EXPECT_NEAR(res[1], c[1], 1e-10);
};
test({.25, .25}, {0, 0});
test({0, 0}, {0, 0});
test({.5, .5}, {1, 1});
test({.25, .5}, {0, 1});
test({.5, .25}, {1, 0});
}
TEST(GRID, natural_coordinates) {
auto mesh = Mesh(1., 2);
auto test = [&mesh](Eigen::Vector2i elem, Eigen::Vector2d pos,
Eigen::Vector2d res) {
auto c = mesh.naturalCoordinates(elem, pos);
EXPECT_NEAR(res[0], c[0], 1e-10);
EXPECT_NEAR(res[1], c[1], 1e-10);
};
test({0, 0}, {.5, .5}, {1, 1});
test({0, 0}, {.25, .25}, {0, 0});
test({1, 1}, {.5, .5}, {-1, -1});
test({0, 1}, {.5, .5}, {1, -1});
test({1, 0}, {.5, .5}, {-1, 1});
test({1, 1}, {.75, .75}, {0, 0});
}
TEST(GRID, creation) {
auto mesh = Mesh(1., 1);
Eigen::Vector4d shp = mesh.getShape({0.5, 0.5});
std::for_each(shp.begin(), shp.end(),
[](Real &v) { EXPECT_NEAR(v, 1. / 4., 1e-14); });
shp = mesh.getShape({0.25, 0.25});
EXPECT_NEAR(shp[0], 9. / 16., 1e-14);
EXPECT_NEAR(shp[1], 3. / 16., 1e-14);
EXPECT_NEAR(shp[2], 1. / 16., 1e-14);
EXPECT_NEAR(shp[3], 3. / 16., 1e-14);
shp = mesh.getShape({0.75, 0.75});
EXPECT_NEAR(shp[0], 1. / 16., 1e-14);
EXPECT_NEAR(shp[1], 3. / 16., 1e-14);
EXPECT_NEAR(shp[2], 9. / 16., 1e-14);
EXPECT_NEAR(shp[3], 3. / 16., 1e-14);
shp = mesh.getShape({0.75, 0.25});
EXPECT_NEAR(shp[0], 3. / 16., 1e-14);
EXPECT_NEAR(shp[1], 9. / 16., 1e-14);
EXPECT_NEAR(shp[2], 3. / 16., 1e-14);
EXPECT_NEAR(shp[3], 1. / 16., 1e-14);
shp = mesh.getShape({0.25, 0.75});
EXPECT_NEAR(shp[0], 3. / 16., 1e-14);
EXPECT_NEAR(shp[1], 1. / 16., 1e-14);
EXPECT_NEAR(shp[2], 3. / 16., 1e-14);
EXPECT_NEAR(shp[3], 9. / 16., 1e-14);
}
ShapeMatrix fill_shp_matrix(Mesh &mesh, Eigen::ArrayX2d &Xs) {
UInt n_points = Xs.rows();
UInt n_nodes = mesh.n_nodes;
ShapeMatrix shp(n_nodes, n_points);
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;
}
auto prepare_mesh_points_shape(int nb_elements, int nb_points) {
auto mesh = Mesh(1., nb_elements);
Eigen::ArrayX2d Xs = Eigen::ArrayX2d::Random(nb_points, 2);
Xs = (Xs + 1) / 2;
ShapeMatrix shp = fill_shp_matrix(mesh, Xs);
return std::make_tuple(mesh, Xs, shp);
}
TEST(SHAPE_MATRIX, filling) {
int n_elem = 10, n_points = 12;
auto &&[mesh, Xs, shp] = prepare_mesh_points_shape(n_elem, n_points);
EXPECT_EQ(shp.cols(), n_points);
EXPECT_EQ(shp.rows(), (n_elem + 1) * (n_elem + 1));
}
TEST(SHAPE_MATRIX, interpolation_constant_field) {
auto &&[mesh, Xs, shp] = prepare_mesh_points_shape(99, 1000);
Eigen::ArrayX2d constant_node_field(mesh.n_nodes, 2);
constant_node_field.col(0) = 1.;
constant_node_field.col(1) = 2.;
Eigen::ArrayX2d interpolated = shp.transpose() * constant_node_field.matrix();
for (auto &&v : interpolated.rowwise()) {
EXPECT_NEAR(v[0], 1., 1e-14);
EXPECT_NEAR(v[1], 2., 1e-14);
}
}
TEST(SHAPE_MATRIX, interpolate_position_field) {
auto &&[mesh, Xs, shp] = prepare_mesh_points_shape(100, 1000);
Eigen::ArrayX2d interpolated = shp.transpose() * mesh.nodes.matrix();
for (auto &&[X, v] : zip(Xs.rowwise(), interpolated.rowwise())) {
EXPECT_NEAR(v[0], X[0], 1e-10);
EXPECT_NEAR(v[1], X[1], 1e-10);
}
}
TEST(SHAPE_MATRIX, interpolate_gradient_field) {
auto &&[mesh, Xs, shp] = prepare_mesh_points_shape(100, 1000);
auto field = [](auto &X) {
Real a = 1, b = 2, c = 3, d = 4;
return Vector<2>{a * X[0] + b, c * X[1] + d};
};
Eigen::ArrayX2d gradient_node_field(mesh.n_nodes, 2);
for (auto &&[X, F] :
zip(mesh.nodes.rowwise(), gradient_node_field.rowwise())) {
F = field(X);
}
Eigen::ArrayX2d interpolated = shp.transpose() * gradient_node_field.matrix();
for (auto &&[X, v] : zip(Xs.rowwise(), interpolated.rowwise())) {
auto res = field(X);
EXPECT_NEAR(v[0], res[0], 1e-10);
EXPECT_NEAR(v[1], res[1], 1e-10);
}
}
TEST(SHAPE_MATRIX, interpolate_gradient_field_libmultiscale) {
auto &&[mesh, Xs, shp] = prepare_mesh_points_shape(100, 1000);
auto field = [](auto &X) {
Real a = 1, b = 2, c = 3, d = 4;
return Vector<2>{a * X[0] + b, c * X[1] + d};
};
ContainerArray<Real> gradient_node_field(mesh.n_nodes, 2);
for (auto &&[X, F] :
zip(mesh.nodes.rowwise(), gradient_node_field.array().rowwise())) {
F = field(X);
}
ContainerArray<Real> interpolated;
interpolated = shp.transpose() * gradient_node_field.matrix();
ContainerArray<Real> interpolated2 =
shp.transpose() * gradient_node_field.matrix();
for (auto &&[X, v] : zip(Xs.rowwise(), interpolated.array().rowwise())) {
auto res = field(X);
EXPECT_NEAR(v[0], res[0], 1e-10);
EXPECT_NEAR(v[1], res[1], 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;
}
}
// 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);
// }

Event Timeline

c4science · Help