Skip to content
Snippets Groups Projects
Commit c1d1f5c3 authored by Porrmann, Maik's avatar Porrmann, Maik
Browse files

Unify VecMatrix classes

parent ecd659b7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/localfunctions/PolynomialBasisCoefficients.hh
+ 57
− 17
View file @ c1d1f5c3
......@@ -7,26 +7,66 @@ namespace Dune
namespace Impl
{
template <class F, unsigned int size>
struct PolynomialBasisVecMatrix
{
using Field = F;
PolynomialBasisVecMatrix() = delete;
// TODO is this a deep copy? if yes change that
// PolynomialBasisVecMatrix(Dune::FieldMatrix<Field, size, size> mat_) : mat(mat_)
// {
// }
PolynomialBasisVecMatrix(Dune::FieldMatrix<Field, size, size> &&mat_) : mat(mat_)
{
}
unsigned int cols() const { return size; };
unsigned int rows() const { return size; };
template <class Vector>
void row(const unsigned int row, Vector &vec) const
{
const unsigned int N = cols();
assert(vec.size() == N);
for (unsigned int i = 0; i < N; ++i)
{
field_cast(mat[row][i], vec[i]);
}
}
Dune::FieldMatrix<Field, size, size> mat;
};
template <class F, unsigned int dim>
using HermiteVecMatrix = PolynomialBasisVecMatrix<F, (dim == 1) ? 4 : (dim == 2) ? 10
: 20>;
template <class F>
using ArgyrisVecMatrix = PolynomialBasisVecMatrix<F, 21>;
template <class F>
using MorleyVecMatrix = PolynomialBasisVecMatrix<F, 6>;
namespace PolynomialBasisCoefficients
{
/**
* @brief Get the Hermite Coefficients Matrix
*
* @tparam F Field type
* @tparam dim dimesion of domain of Reference triangle
* @return auto FieldMatrix<F,size,size> where size depends on dim
*/
* @brief Get the Hermite Coefficients Matrix
*
* @tparam F Field type
* @tparam dim dimesion of domain of Reference triangle
* @return auto FieldMatrix<F,size,size> where size depends on dim
*/
template <class F, int dim>
auto getHermiteCoefficients()
{
if constexpr (dim == 1)
{
return Dune::FieldMatrix<F, 4, 4>{
{1, 0, -3, 2}, {0, 1, -2, 1}, {0, 0, 3, -2}, {0, 0, -1, 1}};
return HermiteVecMatrix<F, 1>({{1, 0, -3, 2}, {0, 1, -2, 1}, {0, 0, 3, -2}, {0, 0, -1, 1}});
}
else if constexpr (dim == 2)
{
return Dune::FieldMatrix<F, 10, 10>{
return HermiteVecMatrix<F, 2>({
{1, 0, 0, -3, -13, -3, 2, 13, 13, 2},
{0, 1, 0, -2, -3, 0, 1, 3, 2, 0},
{0, 0, 1, 0, -3, -2, 0, 2, 3, 1}, // l_2
......@@ -36,11 +76,11 @@ namespace Dune
{0, 0, 0, 0, -7, 3, 0, 7, 7, -2}, // l_6
{0, 0, 0, 0, -1, 0, 0, 1, 2, 0},
{0, 0, 0, 0, 2, -1, 0, -2, -2, 1},
{0, 0, 0, 0, 27, 0, 0, -27, -27, 0}}; // l_9, inner dof
{0, 0, 0, 0, 27, 0, 0, -27, -27, 0}}); // l_9, inner dof
}
else if constexpr (dim == 3)
{
return Dune::FieldMatrix<F, 20, 20>{
return HermiteVecMatrix<F, 3>({
{1, 0, 0, 0, -3, -13, -3, -13, -13, -3, // deg 0 to 2
2, 13, 13, 2, 13, 33, 13, 13, 13, 2}, // deg 3
{0, 1, 0, 0, -2, -3, 0, -3, 0, 0,
......@@ -80,14 +120,14 @@ namespace Dune
{0, 0, 0, 0, 0, 0, 0, 27, 0, 0,
0, 0, 0, 0, -27, -27, 0, -27, 0, 0},
{0, 0, 0, 0, 0, 27, 0, 0, 0, 0,
0, -27, -27, 0, 0, -27, 0, 0, 0, 0}};
0, -27, -27, 0, 0, -27, 0, 0, 0, 0}});
}
else
{
DUNE_THROW(Dune::NotImplemented,
"Hermite-coefficients not implemented for dim == " +
std::to_string(dim));
return -1;
return 1;
}
}
......@@ -98,7 +138,7 @@ namespace Dune
* @return Dune::FieldMatrix<F,21,21>
*/
template <class F>
FieldMatrix<F, 21, 21> getArgyrisCoefficients()
ArgyrisVecMatrix<F> getArgyrisCoefficients()
{
F sqrt2 = -8. * sqrt(2.);
return Dune::FieldMatrix<F, 21, 21>{
......@@ -277,11 +317,11 @@ namespace Dune
* @return Dune::FieldMatrix<F,6,6>
*/
template <class F>
FieldMatrix<F, 6, 6> getMorleyCoefficients()
auto getMorleyCoefficients()
{
F sqrt2 = 0.5 * sqrt(2.);
return Dune::FieldMatrix<F, 6, 6>{
return MorleyVecMatrix<F>{
{1, -1, -1, 0, 2, 0},
{0, 0.5, 0.5, 0.5, -1, -0.5},
{0, 0.5, 0.5, -0.5, -1, 0.5},
......@@ -289,6 +329,6 @@ namespace Dune
{0, -1., 0, 1., 0, 0},
{0, 0, -1, 0, 0, 1}};
}
} // namespace PolBasisCoeff
} // namespace PolynomialBasisCoefficients
} // namespace Impl
} // namespace Dune
\ No newline at end of file
dune/localfunctions/argyris.hh
+ 2
− 34
View file @ c1d1f5c3
......@@ -23,36 +23,6 @@ namespace Dune
{
namespace Impl
{
/** please doc me*/ // #Todo
template <class F> // dim == 2
struct ArgyrisVecMatrix
{
using Field = F;
static constexpr int size = 21;
Dune::FieldMatrix<Field, size, size> mat;
ArgyrisVecMatrix() = delete;
ArgyrisVecMatrix(Dune::FieldMatrix<Field, size, size> coeffs)
: mat(coeffs) {}
unsigned int cols() const
{
return size;
};
unsigned int rows() const { return size; };
template <class Vector>
void row(const unsigned int row, Vector &vec) const
{
const unsigned int N = cols();
assert(vec.size() == N);
for (unsigned int i = 0; i < N; ++i)
{
field_cast(mat[row][i], vec[i]);
}
}
};
/**
* \brief Implementation of Argyris shape Function using
......@@ -96,9 +66,7 @@ namespace Dune
if (dim == 2)
{
ArgyrisVecMatrix<D> coeffs(
PolynomialBasisCoefficients::getArgyrisCoefficients<D>());
this->fill(coeffs);
this->fill(PolynomialBasisCoefficients::getArgyrisCoefficients<D>());
}
else
DUNE_THROW(Dune::NotImplemented, "only implemented for dim == 1");
......@@ -208,7 +176,7 @@ namespace Dune
DUNE_THROW(NotImplemented, "ArgyrisInterpolation only implemented for dim==2!");
}
};
}
} // namespace Impl
/**
* \brief Argyris finite element for triangles
* \tparam D Type used for domain coordinates
......
dune/localfunctions/hermite.hh
+ 1
− 30
View file @ c1d1f5c3
......@@ -24,34 +24,7 @@ namespace Dune
{
namespace Impl
{
template <class F, unsigned int dim>
struct HermiteVecMatrix
{
using Field = F;
static constexpr int size = (dim == 1) ? 4 : (dim == 2) ? 10
: 20;
HermiteVecMatrix() = delete;
// TODO is this a deep copy? if yes change that
HermiteVecMatrix(Dune::FieldMatrix<Field, size, size> mat_) : mat(mat_) {}
unsigned int cols() const { return size; };
unsigned int rows() const { return size; };
template <class Vector>
void row(const unsigned int row, Vector &vec) const
{
const unsigned int N = cols();
assert(vec.size() == N);
for (unsigned int i = 0; i < N; ++i)
{
field_cast(mat[row][i], vec[i]);
}
}
Dune::FieldMatrix<Field, size, size> mat;
};
/**
* \brief Implementation of hermite Polynomials using PolynomialBasisWithMatrix
* \tparam D Type to represent the field in the domain
......@@ -93,9 +66,7 @@ namespace Dune
if (dim <= 3)
{
HermiteVecMatrix<D, dim> coeffs(
PolynomialBasisCoefficients::getHermiteCoefficients<D, dim>());
this->fill(coeffs);
this->fill(PolynomialBasisCoefficients::getHermiteCoefficients<D, dim>());
}
else
DUNE_THROW(Dune::NotImplemented, "only implemented for dim <= 3");
......
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