Commit 0abb4580 authored by Pries, Jason's avatar Pries, Jason

Initial attempt at Torque post processing

Problems with discontinuous material boundaries have fields which suffer from jump discontinuities
This means the divergence of the Maxwell stress tensor is impulsive at material interfaces
Evaluating net force/torque using body forces requires both volumetric and surface integrals in these cases
This is a nontrivial excercise

In the meantime, the next step will be to add a special calculation method for cylindrical airgaps
parent b9c7a73e
......@@ -3,8 +3,7 @@ cmake_minimum_required(VERSION 3.2)
project(Oersted)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage -fopenmp")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 -fopenmp")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage -fopenmp")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
include_directories(./lib/)
......
//
// Created by p7k on 8/23/17.
//
#include "Element.h"
//
// Created by p7k on 8/23/17.
//
#include "Facet.h"
//
// Created by p7k on 8/23/17.
//
#include "Node.h"
std::ostream& operator<<(std::ostream& os, XY const &xy) {
os << xy.Value[0] << "," << xy.Value[1];
return os;
}
\ No newline at end of file
......@@ -28,6 +28,12 @@ public:
void y(double_t val) { Value[1] = val; };
friend double_t hypot(XY const &xy0, XY const &xy1) {
std::hypot(xy0.Value[0] - xy1.Value[0], xy0.Value[1] - xy1.Value[1]);
}
friend std::ostream& operator<<(std::ostream& os, XY const &xy);
protected:
double_t Value[2];
};
......
#include "Triangle.h"
std::fstream &operator<<(std::fstream &fs, Triangle<1> const& t) {
fs << t.node(0) << ',' << t.node(1) << ',' << t.node(2) << std::endl;
return fs;
}
std::fstream &operator<<(std::fstream &fs, Triangle<2> const& t) {
fs << t.node(0) << ',' << t.node(1) << ',' << t.node(5) << std::endl;
fs << t.node(2) << ',' << t.node(3) << ',' << t.node(1) << std::endl;
fs << t.node(4) << ',' << t.node(5) << ',' << t.node(3) << std::endl;
fs << t.node(1) << ',' << t.node(3) << ',' << t.node(5) << std::endl;
return fs;
}
......@@ -20,9 +20,9 @@ public:
Triangle(size_t id, std::array<size_t, NumNodes> const &n) : Facet{id}, Node{n} {};
void node(size_t const &i, size_t n) { Node[i] = n; };
void node(size_t i, size_t n) { Node[i] = n; };
size_t const &node(size_t const &i) const { return Node[i]; };
size_t const &node(size_t i) const { return Node[i]; };
template<size_t Q>
std::array<double_t, NumNodes> basis(size_t q) const;
......@@ -33,6 +33,15 @@ public:
template<size_t Q>
std::array<double_t, NumNodes> db(size_t q) const;
template<size_t Q>
std::array<double_t, NumNodes> daa(size_t q) const;
template<size_t Q>
std::array<double_t, NumNodes> dbb(size_t q) const;
template<size_t Q>
std::array<double_t, NumNodes> dab(size_t q) const;
template<size_t Q>
Oe::Matrix<double_t, 2, 2> jacobian(std::vector<XY> const &nodes, size_t q) const {
Oe::Matrix<double_t, 2, 2> value;
......@@ -106,9 +115,96 @@ public:
auto const &daq = da<Q>(q);
auto const &dbq = db<Q>(q);
for(size_t i = 0; i != NumNodes; ++i) {
df.dx(q, Node[i], ID) += J(0,0) * daq[i] + J(0,1) * dbq[i];
df.dy(q, Node[i], ID) += J(1,0) * daq[i] + J(1,1) * dbq[i];
for (size_t i = 0; i != NumNodes; ++i) {
df.dx(q, Node[i], ID) += J(0, 0) * daq[i] + J(0, 1) * dbq[i];
df.dy(q, Node[i], ID) += J(1, 0) * daq[i] + J(1, 1) * dbq[i];
}
}
}
template<size_t Q>
void second_derivative(SecondDerivativeMatrixGroup &ddf, std::vector<XY> const &nodes) const {
for (size_t q = 0; q != TriangleQuadrature<Q>::size; ++q) {
auto const &daq = da<Q>(q);
auto const &dbq = db<Q>(q);
auto const &daaq = daa<Q>(q);
auto const &dbbq = dbb<Q>(q);
auto const &dabq = dab<Q>(q);
// Calculate derivatives of mapped (x,y) coordinates w.r.t. master (a,b) coordinates
double_t dxda{0.0}, dxdb{0.0}, dyda{0.0}, dydb{0.0};
Oe::Vector3d d2x = Oe::Vector3d::Zero();
Oe::Vector3d d2y = Oe::Vector3d::Zero();
for (size_t i = 0; i != NumNodes; ++i) {
auto const &p = nodes[Node[i]];
dxda += daq[i] * p.x();
dxdb += dbq[i] * p.x();
dyda += daq[i] * p.y();
dydb += dbq[i] * p.y();
d2x(0) += daaq[i] * p.x();
d2x(1) += dbbq[i] * p.x();
d2x(2) += dabq[i] * p.x();
d2y(0) += daaq[i] * p.y();
d2y(1) += dbbq[i] * p.y();
d2y(2) += dabq[i] * p.y();
}
// First derivative Jacobian
Oe::Matrix2d J1;
J1(0, 0) = dxda;
J1(0, 1) = dyda;
J1(1, 0) = dxdb;
J1(1, 1) = dydb;
// J1 = J1.inverse();
// Second derivative Jacobian
Oe::Matrix3d J2;
J2(0, 0) = dxda * dxda;
J2(0, 1) = dyda * dyda;
J2(0, 2) = dxda * dyda * 2.0;
J2(1, 0) = dxdb * dxdb;
J2(1, 1) = dydb * dydb;
J2(1, 2) = dxdb * dydb * 2.0;
J2(2, 0) = dxda * dxdb;
J2(2, 1) = dyda * dydb;
J2(2, 2) = dxda * dydb + dxdb * dyda;
// J2 = J2.inverse();
//
Oe::Vector2d b1;
Oe::Vector2d v1;
Oe::Vector3d b2;
Oe::Vector3d v2;
for (size_t i = 0; i != NumNodes; ++i) {
// Calculate dL_i/dx and dL_i/dy
b1(0) = daq[i];
b1(1) = dbq[i];
v1 = J1.inverse() * b1;
// rhs = d^2 L_i / (da db) - dL_i/dx * dx/(da db) - dL_i/dy * dy/(da db)
b2(0) = daaq[i] - v1(0) * d2x(0) - v1(1) * d2y(0);
b2(1) = dbbq[i] - v1(0) * d2x(1) - v1(1) * d2y(1);
b2(2) = dabq[i] - v1(0) * d2x(2) - v1(1) * d2y(2);
v2 = J2.inverse() * b2;
ddf.dxx(q, Node[i], ID) = v2(0);
ddf.dyy(q, Node[i], ID) = v2(1);
ddf.dxy(q, Node[i], ID) = v2(2);
}
}
}
......@@ -139,6 +235,9 @@ protected:
std::array<size_t, NumNodes> Node;
};
std::fstream &operator<<(std::fstream &fs, Triangle<1> const& t);
std::fstream &operator<<(std::fstream &fs, Triangle<2> const& t);
template<>
template<size_t Q>
inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::basis(size_t q) const {
......@@ -196,6 +295,42 @@ inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::db(size_t q) con
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::daa(size_t q) const {
std::array<double_t, NumNodes> value;
value[0] = 0.0;
value[1] = 0.0;
value[2] = 0.0;
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::dbb(size_t q) const {
std::array<double_t, NumNodes> value;
value[0] = 0.0;
value[1] = 0.0;
value[2] = 0.0;
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::dab(size_t q) const {
std::array<double_t, NumNodes> value;
value[0] = 0.0;
value[1] = 0.0;
value[2] = 0.0;
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::da(size_t q) const {
......@@ -232,4 +367,58 @@ inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::db(size_t q) con
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::daa(size_t q) const {
std::array<double_t, NumNodes> value;
double_t const &a = TriangleQuadrature<Q>::a[q];
double_t const &b = TriangleQuadrature<Q>::b[q];
value[0] = +4.0;
value[1] = -8.0;
value[2] = +4.0;
value[3] = +0.0;
value[4] = +0.0;
value[5] = +0.0;
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::dbb(size_t q) const {
std::array<double_t, NumNodes> value;
double_t const &a = TriangleQuadrature<Q>::a[q];
double_t const &b = TriangleQuadrature<Q>::b[q];
value[0] = +4.0;
value[1] = +0.0;
value[2] = +0.0;
value[3] = +0.0;
value[4] = +4.0;
value[5] = -8.0;
return value;
}
template<>
template<size_t Q>
inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::dab(size_t q) const {
std::array<double_t, NumNodes> value;
double_t const &a = TriangleQuadrature<Q>::a[q];
double_t const &b = TriangleQuadrature<Q>::b[q];
value[0] = +4.0;
value[1] = -4.0;
value[2] = +0.0;
value[3] = +4.0;
value[4] = +0.0;
value[5] = -4.0;
return value;
}
#endif //OERSTED_TRIANGLE_H
\ No newline at end of file
......@@ -15,12 +15,14 @@ class MagneticMaterialInterface {
public:
virtual ~MagneticMaterialInterface(){};
virtual void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy) = 0;
virtual void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) = 0;
};
class LinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
LinearIsotropicMagneticMaterial(double relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};
LinearIsotropicMagneticMaterial(double_t relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};
void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) override {
for(size_t const &i : index) {
......@@ -32,6 +34,13 @@ public:
}
}
void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy) override {
for(size_t const &i : index) {
Fx(i) *= Reluctivity;
Fy(i) *= Reluctivity;
}
}
protected:
double Reluctivity;
};
......@@ -119,6 +128,37 @@ public:
}
}
void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy) override {
for(size_t const &i : index) {
const double_t &Bx = Fx(i);
const double_t &By = Fy(i);
double_t B = sqrt(Bx * Bx + By * By);
double_t M;
if (B > BSaturation) {
M = MSaturation;
} else {
size_t j = (size_t)(B / DeltaB);
const auto &s = MBSpline[j];
M = (s[0] * B + s[1]) * B + s[2];
}
if (B > 0.0) {
double_t Mx = Bx * (M / B);
double_t My = By * (M / B);
Fx(i) = Bx / mu_0 - Mx;
Fy(i) = By / mu_0 - My;
} else {
Fx(i) = 0.0;
Fy(i) = 0.0;
}
}
}
protected:
std::vector<std::array<double_t,3>> MBSpline;
double_t DeltaB;
......
......@@ -2,7 +2,7 @@
#define OERSTED_PHYSICSCONSTANTS_H
static constexpr double speed_of_light = 299792458.0;
static constexpr double mu_0 = 4 * M_PI * 1e-7;
static constexpr double mu_0 = 4.0 * M_PI * 1e-7;
static constexpr double epsilon_0 = 1.0 / (mu_0 * speed_of_light * speed_of_light);
#endif //OERSTED_PHYSICSCONSTANTS_H
......@@ -9,8 +9,7 @@ set(SOURCE_FILES
./src/SparseMatrixGroup.cpp ./src/SparseMatrixGroup.h
./src/DiagonalMatrixGroup.cpp ./src/DiagonalMatrixGroup.h
./src/DerivativeMatrixGroup.cpp ./src/DerivativeMatrixGroup.h
)
./src/SecondDerivativeMatrixGroup.cpp ./src/SecondDerivativeMatrixGroup.h)
add_library(matrix SHARED ${SOURCE_FILES})
......
#include "SecondDerivativeMatrixGroup.h"
#ifndef OERSTED_SECONDDERIVATIVEMATRIXGROUP_H
#define OERSTED_SECONDDERIVATIVEMATRIXGROUP_H
#include "MatrixGroup.h"
#include "SparseMatrixGroup.h"
class SecondDerivativeMatrixGroup {
public:
SecondDerivativeMatrixGroup() {};
SecondDerivativeMatrixGroup(size_t const qpts, size_t const nodes, size_t const nels, size_t const el_nnz) :
Dxx{qpts, nodes, nels, el_nnz}, Dyy{qpts, nodes, nels, el_nnz} , Dxy{qpts, nodes, nels, el_nnz} {};
auto const &dxx(size_t const q) const { return Dxx(q); };
auto const &dyy(size_t const q) const { return Dyy(q); };
auto const &dxy(size_t const q) const { return Dxy(q); };
auto const &dyx(size_t const q) const { return Dxy(q); };
auto &dxx(size_t const q, size_t const i, size_t const j) { return Dxx(q, i, j); };
auto &dyy(size_t const q, size_t const i, size_t const j) { return Dyy(q, i, j); };
auto &dxy(size_t const q, size_t const i, size_t const j) { return Dxy(q, i, j); };
auto &dyx(size_t const q, size_t const i, size_t const j) { return Dxy(q, i, j); };
void transform(Oe::SparseMatrix<double_t> &A) {
Dxx.transform(A);
Dyy.transform(A);
Dxy.transform(A);
}
size_t size() const { return Dxx.size(); };
protected:
SparseMatrixGroup Dxx;
SparseMatrixGroup Dyy;
SparseMatrixGroup Dxy;
};
#endif //OERSTED_SECONDDERIVATIVEMATRIXGROUP_H
#ifndef OERSTED_MESH_H
#define OERSTED_MESH_H
// TODO: Ensure output mesh is independent of memory location of sketch elements. Something (potentially sorting by shared_ptr address) could be causing discrepancies
#include <algorithm>
#include <fstream>
#include <numeric>
......@@ -49,15 +51,13 @@ public:
};
template<size_t P>
class MeshData{
class MeshData {
public:
std::vector<Point> Nodes;
// std::vector<std::vector<std::vector<size_t>>> Edges;
std::vector<Triangle<P>> Triangles;
std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> Boundaries;
std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t, size_t>>> Boundaries;
std::map<std::shared_ptr<const Contour>, std::vector<size_t>> Regions;
};
......@@ -239,7 +239,7 @@ public:
MeshData<P> mesh_data() const {
std::vector<NodeData> data;
std::vector<Point> points;
std::vector<Triangle<P>> triangles;
std::vector<Triangle<P>> tris;
std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> boundaries;
std::map<std::shared_ptr<const Contour>, std::vector<size_t>> regions;
......@@ -283,22 +283,22 @@ public:
}
std::sort(data.begin(), data.end());
triangles.resize(size_triangles());
tris.resize(size_triangles());
size_t i = 0;
size_t n = 0;
while (i < data.size()) {
points.push_back(data[i].Coordinates);
size_t k = data[i].EdgeID * P + data[i].NodeID;
triangles[data[i].ElementID].node(k,n);
size_t k = (data[i].EdgeID * P + data[i].NodeID) % Triangle<P>::NumNodes;
tris[data[i].ElementID].node(k,n);
boundaries[data[i].Constraint].emplace(data[i].Parameter,n);
size_t j = i + 1;
while (j < data.size() && data[i].Coordinates == data[j].Coordinates) {
size_t k = data[j].EdgeID * P + data[j].NodeID;
triangles[data[j].ElementID].node(k,n);
size_t k = (data[j].EdgeID * P + data[j].NodeID) % Triangle<P>::NumNodes;
tris[data[j].ElementID].node(k,n);
boundaries[data[j].Constraint].emplace(data[j].Parameter,n);
......@@ -309,11 +309,11 @@ public:
++n;
}
for(size_t i = 0; i != triangles.size(); ++i) {
triangles[i].id(i);
for(size_t i = 0; i != tris.size(); ++i) {
tris[i].id(i);
}
return MeshData<P>{points,triangles,boundaries,regions};
return MeshData<P>{points,tris,boundaries,regions};
}
std::vector<size_t> boundary_nodes(size_t i) const { return BoundaryConstraints[i]->nodes(*this); }
......
......@@ -21,6 +21,8 @@ public:
auto end() const { return Map.end(); }
VariableMap &operator[](size_t i) { return Map[i]; };
protected:
bool Antiperiodic{false};
......
......@@ -44,6 +44,21 @@ DiagonalMatrixGroup FiniteElementMesh<P, Q>::determinant() const {
return mat;
}
template<size_t P, size_t Q>
DiagonalMatrixGroup FiniteElementMesh<P, Q>::weights() const {
DiagonalMatrixGroup mat(TriangleQuadrature<Q>::size, Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].template determinant<Q>(mat, Nodes);
}
for (size_t q= 0; q != TriangleQuadrature<Q>::size; ++q) {
mat(q) *= TriangleQuadrature<Q>::w[q];
}
return mat;
}
template<size_t P, size_t Q>
SparseMatrixGroup FiniteElementMesh<P, Q>::basis() const {
SparseMatrixGroup mat(TriangleQuadrature<Q>::size, Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
......@@ -66,6 +81,18 @@ DerivativeMatrixGroup FiniteElementMesh<P, Q>::derivative() const {
return df;
}
template<size_t P, size_t Q>
SecondDerivativeMatrixGroup FiniteElementMesh<P, Q>::second_derivative() const {
// TODO: Template based on derivative order template<size_t P, size_t Q, size_t O> DerivativeMatrixGroup<O> FiniteELementMesh<P,Q>::derivative<O>() const {...};
SecondDerivativeMatrixGroup df(TriangleQuadrature<Q>::size, Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].template second_derivative<Q>(df, Nodes);
}
return df;
}
template<size_t P, size_t Q>
std::vector<std::vector<XY>> FiniteElementMesh<P, Q>::quadrature_points() const {
std::vector<std::vector<XY>> qp;
......@@ -127,11 +154,11 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_dis
double_t dy0 = Nodes[boundary_nodes[ii]].y() - Nodes[boundary_nodes[i]].y();
for (Triangle<P> &t : Triangles) {
for (size_t k = 0; k != 3; ++k) { // TODO: This essentially assumes Element order = 1
for (size_t k = 0; k != 3 * P; ++k) { // TODO: This might not work for Element order > 2 having internal nodes, node numbering should occur in a spiral
if (t.node(k) == boundary_nodes[i]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
size_t k0 = t.node(0 * P);
size_t k1 = t.node(1 * P);
size_t k2 = t.node(2 * P);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
......@@ -143,9 +170,9 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_dis
}
if (t.node(k) == boundary_nodes[ii]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
size_t k0 = t.node(0 * P);
size_t k1 = t.node(1 * P);
size_t k2 = t.node(2 * P);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[boundary_nodes[i]].x() + Nodes[boundary_nodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[boundary_nodes[i]].y() + Nodes[boundary_nodes[ii]].y()) / 2.0;
......@@ -253,8 +280,8 @@ void FiniteElementMesh<P, Q>::write_scalar(Oe::VectorXd const &v, std::string pa
// Write header
size_t loc = 3;
fs << loc << ',' << Triangles.size() << std