Commit e7e06b95 authored by JasonPries's avatar JasonPries
Browse files

Partial commit of work done at home

parent 0a553c0d
#include "Triangle.h"
\ No newline at end of file
#ifndef OERSTED_TRIANGLE_H
#define OERSTED_TRIANGLE_H
#include "Eigen"
#include "Eigen/Sparse"
#include "QuadratureRule.h"
#include "Node.h"
// TODO: Curved elements
template<size_t P>
class Triangle {
public:
static constexpr size_t N{(P + 1) * (P + 2) / 2};
Triangle() : Node{} {};
size_t const &node(size_t const &i) { return Node[i]; };
template<typename... Args>
Triangle(Args &&... args) : Node{std::forward<Args>(args)...} {};
template<size_t D>
Eigen::Matrix<double, 2, 2> jacobian(std::vector<XY> const &nodes) const; // TODO: Accept Eigen::Matrix<double,2,2> as input
template<size_t Q>
void basis(std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, TriangleQuadratureRule<Q>::size> &sp_arr, size_t &row);
template<size_t Q>
std::array<std::array<XY, N>, TriangleQuadratureRule<Q>::size> derivative(std::vector<XY> const &nodes) const; // TODO: Rewrite ala ::basis(...);
protected:
size_t Node[N];
};
template<>
template<>
Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &nodes) const {
Eigen::Matrix<double, 2, 2> value;
value(0, 0) = 1.0;
value(0, 1) = 0.0;
value(1, 0) = 0.0;
value(1, 1) = 1.0;
return value;
}
template<>
template<>
Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
Eigen::Matrix<double, 2, 2> value;
XY const &p0 = nodes[Node[0]];
XY const &p1 = nodes[Node[1]];
XY const &p2 = nodes[Node[2]];
double xx = p0.x() - p2.x();
double xy = p0.y() - p2.y();
double yx = p1.x() - p2.x();
double yy = p1.y() - p2.y();
double det = xx * yy - xy * yx;
value(0, 0) = yy / det;
value(0, 1) = -xy / det;
value(1, 0) = -yx / det;
value(1, 1) = xx / det;
return value;
};
template<>
template<size_t Q>
void Triangle<1>::basis(std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, TriangleQuadratureRule<Q>::size> &sp_arr, size_t &row) {
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
sp_arr[i].coeffRef(row, Node[0]) += TriangleQuadratureRule<Q>::a[i];
sp_arr[i].coeffRef(row, Node[1]) += TriangleQuadratureRule<Q>::b[i];
sp_arr[i].coeffRef(row, Node[2]) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
}
}
template<>
template<size_t Q>
std::array<std::array<XY, Triangle<1>::N>, TriangleQuadratureRule<Q>::size> Triangle<1>::derivative(std::vector<XY> const &nodes) const {
std::array<std::array<XY, Triangle<1>::N>, TriangleQuadratureRule<Q>::size> value{};
auto J = jacobian<1>(nodes);
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
value[i][0].x(J(0, 0));
value[i][0].y(J(1, 0));
value[i][1].x(J(0, 1));
value[i][1].y(J(1, 1));
value[i][2].x(-J(0, 0) - J(0, 1));
value[i][2].y(-J(1, 0) - J(1, 1));
}
return value;
};
#endif //OERSTED_TRIANGLE_H
\ No newline at end of file
......@@ -37,7 +37,8 @@ set(SOURCE_FILES
./src/Star.h ./src/Star.cpp
./src/Constellation.h ./src/Constellation.cpp
./src/Contour.h ./src/Contour.cpp)
./src/Contour.h ./src/Contour.cpp
)
add_library(sketch SHARED ${SOURCE_FILES})
......
......@@ -23,6 +23,13 @@ set(SOURCE_FILES
Mesh/util.h
Mesh/util.cpp
Physics/test_Physics.hpp
Physics/test_FiniteElementMesh.cpp
#Physics/test_Discretization.cpp
#Physics/test_Curl.cpp
#Physics/test_Physics.cpp
UseCases/test_UseCases.hpp
UseCases/test_Stator.cpp
UseCases/test_Rotor.cpp)
......
#include "test_Physics.hpp"
inline constexpr size_t operator "" _zu(unsigned long long int i) {
return (size_t) i;
}
class Master_Triangle_1 : public ::testing::Test {
public:
virtual void SetUp() {
node.push_back(XY{0.0, 0.0});
node.push_back(XY{1.0, 0.0});
node.push_back(XY{0.0, 1.0});
node.push_back(XY{-1.0, 0.0});
node.push_back(XY{0.0, -1.0});
tri.push_back(Triangle<1>(1_zu, 2_zu, 0_zu)); // Master
tri.push_back(Triangle<1>(2_zu, 3_zu, 0_zu)); // 90 degree
tri.push_back(Triangle<1>(3_zu, 4_zu, 0_zu)); // 180 degree
tri.push_back(Triangle<1>(4_zu, 1_zu, 0_zu)); // 270 degree
tri.push_back(Triangle<1>(3_zu, 2_zu, 0_zu)); // y-axis reflection
tri.push_back(Triangle<1>(1_zu, 4_zu, 0_zu)); // x-axis reflection
tri.push_back(Triangle<1>(4_zu, 3_zu, 0_zu)); // x=y reflection
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
};
TEST_F(Master_Triangle_1, jacobian_1) {
auto J = tri[0].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);
J = tri[1].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
J = tri[2].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);
J = tri[3].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), 1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
J = tri[4].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);
J = tri[5].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);
J = tri[6].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
}
TEST_F(Master_Triangle_1, basis_1) {
std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, 1> sp_arr;
sp_arr[0].resize(tri.size(), node.size());
sp_arr[0].reserve(Eigen::VectorXi::Constant(tri.size(),6));
for (size_t i = 0;i != tri.size();++i) {
tri[i].basis<1>(sp_arr, i);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(0)), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(1)), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(2)), 1.0 / 3.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_1) {
auto d = tri[0].derivative<1>(node);
EXPECT_NEAR(d[0][0].x(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][0].y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[0][1].x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[0][1].y(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][2].x(), -1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][2].y(), -1.0, FLT_EPSILON);
}
TEST_F(Master_Triangle_1, basis_2) {
std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, 3> sp_arr;
for (size_t i = 0;i!=3;++i) {
sp_arr[i].resize(tri.size(), node.size());
sp_arr[i].reserve(Eigen::VectorXi::Constant(tri.size(),6));
}
for (size_t i = 0; i != tri.size(); ++i) {
tri[i].basis<2>(sp_arr, i);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(0)), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(1)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(2)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(0)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(1)), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(2)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(0)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(1)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(2)), 2.0 / 3.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_2) {
auto d = tri[0].derivative<2>(node);
for (size_t i = 0; i != 3; ++i) {
EXPECT_NEAR(d[i][0].x(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][0].y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[i][1].x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[i][1].y(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][2].x(), -1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][2].y(), -1.0, FLT_EPSILON);
}
}
\ No newline at end of file
#ifndef OERSTED_TEST_PHYSICS_HPP
#define OERSTED_TEST_PHYSICS_HPP
#include <cstdio>
#include "gtest.h"
#include "Physics.hpp"
#endif //OERSTED_TEST_PHYSICS_HPP
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment