MatrixDiagonal.cpp
e7dc2669df68neighbour

MatrixDiagonal.cpp
View Options


			#include <catch2/catch.hpp>
			#include <xtensor/xrandom.hpp>
			#include <xtensor/xmath.hpp>
			#include <GooseFEM/GooseFEM.h>

			#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12));

			TEST_CASE("GooseFEM::MatrixDiagonal", "MatrixDiagonal.h")
			{

			SECTION("dot")
			{
			GooseFEM::Mesh::Quad4::Regular mesh(2, 2);

			xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nnode() * mesh.ndim()});
			xt::xtensor<double, 1> b = xt::random::rand<double>({mesh.nnode() * mesh.ndim()});
			xt::xtensor<double, 1> c = a * b;

			GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
			A.set(a);
			xt::xtensor<double, 1> C = A.Dot(b);

			REQUIRE(C.size() == c.size());
			REQUIRE(xt::allclose(C, c));
			}

			SECTION("solve")
			{
			GooseFEM::Mesh::Quad4::Regular mesh(2, 2);

			xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nnode() * mesh.ndim()});
			xt::xtensor<double, 1> b = xt::random::rand<double>({mesh.nnode() * mesh.ndim()});
			xt::xtensor<double, 1> c = a * b;

			GooseFEM::MatrixDiagonal A(mesh.conn(), mesh.dofs());
			A.set(a);
			xt::xtensor<double, 1> C = A.Dot(b);
			xt::xtensor<double, 1> B = A.Solve(C);

			REQUIRE(B.size() == b.size());
			REQUIRE(xt::allclose(B, b));
			}
			}