Commit 46ec5042 authored by JasonPries's avatar JasonPries
Browse files

Fix quadrature compilation issues with gcc 6.3.1-1

Definition of static constexpr double[] in TriangleQuadrature structs were behaving differently between different versions of gcc. The latest gcc 6.3.1-1 didn't like the fact that the ODR was being violated, whereas some previous versions didn't care.

The definitions need to be visible in the Physics headers, so the quadrature rules were moved to a separate library and the definitions were moved to a cpp file.
parent e661fd97
......@@ -18,6 +18,7 @@ add_subdirectory(./lib/)
include_directories(./src/)
include_directories(./src/Sketch/include/)
include_directories(./src/Mesh/include/)
include_directories(./src/Quadrature/include/)
include_directories(./src/Physics/include/)
add_subdirectory(./src/)
......
include_directories(./Sketch/)
include_directories(./Mesh/)
include_directories(./Quadrature/)
include_directories(./Physics/)
add_subdirectory(./Sketch/)
add_subdirectory(./Mesh/)
add_subdirectory(./Quadrature/)
add_subdirectory(./Physics/)
\ No newline at end of file
......@@ -5,6 +5,10 @@ set(SOURCE_FILES
#./src/PhysicsCommon.h
./src/Boundary.h ./src/Boundary.cpp
./src/BoundaryCondition.h ./src/BoundaryCondition.cpp
./src/Forcing.h ./src/Forcing.cpp
./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
......@@ -19,8 +23,6 @@ set(SOURCE_FILES
#./src/Operator.h ./src/Operator.cpp
./src/QuadratureRule.h ./src/QuadratureRule.cpp
./src/Region.h ./src/Region.cpp
./src/Triangle.h ./src/Triangle.cpp
......@@ -28,4 +30,4 @@ set(SOURCE_FILES
add_library(physics SHARED ${SOURCE_FILES})
#target_link_libraries(sketch ${Boost_LIBRARIES})
\ No newline at end of file
target_link_libraries(physics ${Boost_LIBRARIES})
\ No newline at end of file
......@@ -14,7 +14,7 @@
#include "../src/Node.h"
#include "../src/Physics.h"
#include "../src/PhysicsConstants.h"
#include "../src/QuadratureRule.h"
#include "Quadrature/src/TriangleQuadrature.h"
#include "../src/Region.h"
#include "../src/Triangle.h"
......
#ifndef OERSTED_BOUNDARY_H
#define OERSTED_BOUNDARY_H
#include <vector>
#include <cstddef>
template<size_t Dimension>
class Boundary {
};
......
......@@ -7,8 +7,9 @@
#include <Eigen/Sparse>
#include "Quadrature.hpp"
#include "MatrixGroup.h"
#include "QuadratureRule.h"
#include "FiniteElementMesh.h"
// TODO: DiscreteBoundaryCondition versus ContinuousBoundaryCondition? Somehow need to apply to top level model and have that be translated to the discrete model
......
#include "FiniteElementMesh.h"
#include "FiniteElementMesh.h"
\ No newline at end of file
......@@ -4,6 +4,7 @@
#include <memory>
#include <vector>
#include "Boundary.h"
#include "MatrixGroup.h"
#include "Node.h"
#include "Triangle.h"
......@@ -50,8 +51,8 @@ public:
auto const &boundary(size_t i) const { return Boundaries[i]; };
template<size_t Q>
DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> determinant() const {
DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Triangles.size());
DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
DiagonalMatrixGroup<TriangleQuadrature<Q>::size> mat(Triangles.size());
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].determinant<Q>(mat, Nodes);
......@@ -61,8 +62,8 @@ public:
}
template<size_t Q>
SparseMatrixGroup<TriangleQuadratureRule<Q>::size> basis() const {
SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
SparseMatrixGroup<TriangleQuadrature<Q>::size> basis() const {
SparseMatrixGroup<TriangleQuadrature<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].basis<Q>(mat, Nodes);
......@@ -72,8 +73,8 @@ public:
}
template<size_t Q>
DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> derivative() const {
DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> derivative() const {
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
for (size_t i = 0; i!= Triangles.size();++i) {
Triangles[i].derivative<Q>(df, Nodes);
......
......@@ -5,9 +5,11 @@
#include <vector>
#include "Eigen"
#include "Quadrature.hpp"
#include "Region.h"
#include "MatrixGroup.h"
#include "QuadratureRule.h"
#include "FiniteElementMesh.h"
class Forcing {
......@@ -26,8 +28,8 @@ class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Forcing {
public:
HomogeneousForcing(std::function<double(double)> f,
std::vector<std::shared_ptr<Region<2>>> regions,
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &basis,
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &weights)
SparseMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &basis,
DiagonalMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &weights)
: Function{f}, Regions{regions}, Basis(basis), Weights(weights) {
// TODO: May need to adjust implementation once quadrilaterals are added
......@@ -42,7 +44,7 @@ public:
Eigen::VectorXd operator()(double const t) const override {
Eigen::VectorXd output = (Basis[0] * (Weights(0).matrix().asDiagonal() * Indicator)); // TODO: ?SparseVector?
for (size_t i = 1; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
for (size_t i = 1; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
output += (Basis[i] * (Weights(i).matrix().asDiagonal() * Indicator));
}
......@@ -57,8 +59,8 @@ protected:
std::vector<std::shared_ptr<Region<2>>> Regions;
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Basis;
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Weights;
SparseMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &Basis;
DiagonalMatrixGroup<TriangleQuadrature<QuadratureOrder>::size> const &Weights;
Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
};
......
#include "MaterialProperties.h"
#include "MaterialProperties.h"
\ No newline at end of file
......@@ -50,7 +50,7 @@ protected:
std::shared_ptr<MagneticMaterialInterface> MagneticProperties;
};
MaterialProperties Air() {
inline MaterialProperties Air() {
return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
}
......
......@@ -3,6 +3,8 @@
// TODO: Const
#include <cstddef>
template<size_t N>
class SmallVector {
private:
......
......@@ -12,7 +12,6 @@
#include "Forcing.h"
#include "Node.h"
#include "PhysicsConstants.h"
#include "QuadratureRule.h"
#include "Triangle.h"
// TODO: These classes should be generic (e.g. CurlCurl, ScalarLaplacian, VectorLaplacian)
......@@ -54,7 +53,7 @@ protected:
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class PhysicsData {
public:
static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
static constexpr size_t QuadraturePoints = TriangleQuadrature<QuadratureOrder>::size;
/* TODO: template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder> class PhysicsData : public PhysicsInterface {...};
virtual Eigen::VectorXd init_residual() = 0;
......@@ -118,7 +117,7 @@ public:
r = -f;
J.setZero();
for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
// Caclulate flux density
Fx = Derivatives.dy(i).transpose() * v;
Fy = -Derivatives.dx(i).transpose() * v;
......@@ -146,8 +145,8 @@ public:
Basis = Domain.template basis<QuadratureOrder>();
Derivatives = Domain.template derivative<QuadratureOrder>();
for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
ElementWeights(i) *= TriangleQuadratureRule<QuadratureOrder>::w[i];
for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
ElementWeights(i) *= TriangleQuadrature<QuadratureOrder>::w[i];
}
};
......
#include "QuadratureRule.h"
\ No newline at end of file
......@@ -6,9 +6,10 @@
#include "Eigen"
#include "Eigen/Sparse"
#include "Quadrature.hpp"
#include "MatrixGroup.h"
#include "Node.h"
#include "QuadratureRule.h"
// TODO: Curved elements
......@@ -46,13 +47,13 @@ public:
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 determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const;
void determinant(DiagonalMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const;
template<size_t Q>
void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const;
void basis(SparseMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const;
template<size_t Q>
void derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes) const;
void derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size> &df, std::vector<XY> const &nodes) const;
protected:
size_t Node[NumNodes];
......@@ -61,7 +62,7 @@ protected:
/*
template<>
template<>
Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &nodes) const {
inline 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;
......@@ -75,7 +76,7 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &node
template<>
template<>
Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
inline 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]];
......@@ -98,8 +99,8 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &node
template<>
template<size_t Q>
void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const {
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const {
for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
XY const &p0 = nodes[Node[0]];
XY const &p1 = nodes[Node[1]];
XY const &p2 = nodes[Node[2]];
......@@ -112,22 +113,23 @@ void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::siz
}
}
template<>
template<size_t Q>
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const {
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
mat(i, Node[0], ID) += TriangleQuadratureRule<Q>::a[i];
mat(i, Node[1], ID) += TriangleQuadratureRule<Q>::b[i];
mat(i, Node[2], ID) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const {
for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
mat(i, Node[0], ID) += TriangleQuadrature<Q>::a[i];
mat(i, Node[1], ID) += TriangleQuadrature<Q>::b[i];
mat(i, Node[2], ID) += 1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i];
}
}
template<>
template<size_t Q>
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes) const {
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size> &df, std::vector<XY> const &nodes) const {
auto J = jacobian<1>(nodes);
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
df.dx(i, Node[0], ID) += J(0, 0);
df.dy(i, Node[0], ID) += J(1, 0);
......
project(Oersted_Quadrature)
set(SOURCE_FILES
./include/Quadrature.hpp
./src/TriangleQuadrature.h ./src/TriangleQuadrature.cpp
)
add_library(quadrature SHARED ${SOURCE_FILES})
\ No newline at end of file
#ifndef OERSTED_QUADRATURE_H
#define OERSTED_QUADRATURE_H
#include "../src/TriangleQuadrature.h"
#endif //OERSTED_QUADRATURE_H
#include "TriangleQuadrature.h"
constexpr double TriangleQuadrature<1>::a[];
constexpr double TriangleQuadrature<1>::b[];
constexpr double TriangleQuadrature<1>::w[];
constexpr double TriangleQuadrature<2>::a[];
constexpr double TriangleQuadrature<2>::b[];
constexpr double TriangleQuadrature<2>::w[];
\ No newline at end of file
#ifndef OERSTED_QUADRATURERULE_H
#define OERSTED_QUADRATURERULE_H
#ifndef OERSTED_TRIANGLEQUADRATURE_H
#define OERSTED_TRIANGLEQUADRATURE_H
#include <cstddef>
template<size_t Q>
struct TriangleQuadratureRule {
struct TriangleQuadrature {
};
template<>
struct TriangleQuadratureRule<1> {
struct TriangleQuadrature<1> {
static constexpr size_t size{1};
static constexpr double a[1]{1.0 / 3.0};
static constexpr double b[1]{1.0 / 3.0};
static constexpr double w[1]{1.0 / 2.0};
};
constexpr double TriangleQuadratureRule<1>::a[];
constexpr double TriangleQuadratureRule<1>::b[];
constexpr double TriangleQuadratureRule<1>::w[];
template<>
struct TriangleQuadratureRule<2> {
struct TriangleQuadrature<2> {
static constexpr size_t size{3};
static constexpr double a[3]{2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0};
static constexpr double b[3]{1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0};
static constexpr double w[3]{1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0};
};
constexpr double TriangleQuadratureRule<2>::a[];
constexpr double TriangleQuadratureRule<2>::b[];
constexpr double TriangleQuadratureRule<2>::w[];
#endif //OERSTED_QUADRATURERULE_H
#endif //OERSTED_TRIANGLEQUADRATURE_H
......@@ -36,4 +36,4 @@ set(SOURCE_FILES
add_executable(run_tests ${SOURCE_FILES})
target_link_libraries(run_tests gtest gtest_main sketch mesh)
\ No newline at end of file
target_link_libraries(run_tests gtest gtest_main sketch mesh physics quadrature)
\ No newline at end of file
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