Commit f81fe982 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

removed some bugs in operator evaluation with dim not dow

parent 27ef2c6a
Pipeline #522 failed with stage
in 2 minutes and 23 seconds
......@@ -275,7 +275,8 @@ namespace AMDiS
using Element = typename RowView::Element;
auto element = rowView.element();
auto geometry = element.geometry();
const int dim = Element::dimension;
const int dim = RowView::GridView::dimension;
const int dow = RowView::GridView::dimensionworld;
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
......@@ -306,7 +307,7 @@ namespace AMDiS
colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > colGradients(colReferenceGradients.size());
std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());
for (size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
......@@ -335,7 +336,8 @@ namespace AMDiS
using Element = typename RowView::Element;
auto element = rowView.element();
auto geometry = element.geometry();
const int dim = Element::dimension;
const int dim = RowView::GridView::dimension;
const int dow = RowView::GridView::dimensionworld;
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
......@@ -361,7 +363,7 @@ namespace AMDiS
rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > rowGradients(rowReferenceGradients.size());
std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
for (size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
......@@ -392,7 +394,8 @@ namespace AMDiS
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
const int dim = RowView::GridView::dimension;
const int dow = RowView::GridView::dimensionworld;
auto geometry = element.geometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
......@@ -418,7 +421,7 @@ namespace AMDiS
rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > rowGradients(rowReferenceGradients.size());
std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
for (size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
......@@ -442,7 +445,8 @@ namespace AMDiS
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
const int dim = RowView::GridView::dimension;
const int dow = RowView::GridView::dimensionworld;
auto geometry = element.geometry();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
......@@ -472,7 +476,7 @@ namespace AMDiS
rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
std::vector<Dune::FieldVector<double,dow> > gradients(referenceGradients.size());
for (size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
......
......@@ -15,73 +15,112 @@ K sum(Dune::FieldVector<K, SIZE> const& vector)
result += vector[i];
return result;
}
template <class T, int N>
Dune::FieldMatrix<T, N, N> outer(Dune::FieldVector<T,N> const& vec1, Dune::FieldVector<T,N> const& vec2)
{
Dune::FieldMatrix<T, N, N> mat;
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j)
mat[i][j] = vec1[i] * vec2[j];
return mat;
}
template <class T, int N>
T inner(Dune::FieldVector<T,N> const& vec1, Dune::FieldVector<T,N> const& vec2)
{
T result = 0;
for (int i = 0; i < N; ++i)
result += vec1[i] * vec2[i];
return result;
}
template <class T, int N>
Dune::FieldVector<T,N> operator*(Dune::FieldVector<T,N> vec1, T scalar)
{
return vec1 *= scalar;
}
template <class T, int N>
Dune::FieldVector<T,N> operator*(T scalar, Dune::FieldVector<T,N> vec1)
{
return vec1 *= scalar;
}
template <class T, int N, int M>
Dune::FieldVector<T,N> operator*(Dune::FieldMatrix<T,N,M> const& mat, Dune::FieldVector<T,M> const& vec)
{
return Dune::FMatrixHelp::mult(mat, vec);
}
class OperatorEvaluation
{
using Scalar = Dune::FieldVector<double, 1>;
protected:
template <class... Args>
auto evalZotImpl(Args...) const { assert( false ); return 0; }
template <class ValueType>
auto evalZotImpl(tag::scalar, tag::none, ValueType const& v,
auto evalZotImpl(tag::scalar, tag::none, ValueType const& v,
Scalar const& test, Scalar const& trial) const
{
return v * test * trial;
}
template <class... Args>
auto evalFotImpl(Args...) const { assert( false ); return 0; }
template <size_t I, class ValueType, class Gradient>
auto evalFotImpl(tag::scalar, VectorComponent<I>, ValueType const& v,
auto evalFotImpl(tag::scalar, VectorComponent<I>, ValueType const& v,
Gradient const& grad_test, Scalar const& trial) const
{
return v * grad_test[I] * trial;
}
template <class ValueType, class Gradient>
auto evalFotImpl(tag::vector, tag::none, ValueType const& v,
auto evalFotImpl(tag::vector, tag::none, ValueType const& v,
Gradient const& grad_test, Scalar const& trial) const
{
return (v * grad_test) * trial;
}
template <class ValueType, class Gradient>
auto evalFotImpl(tag::scalar, tag::none, ValueType const& v,
auto evalFotImpl(tag::scalar, tag::none, ValueType const& v,
Gradient const& grad_test, Scalar const& trial) const
{
return v * sum(grad_test) * trial;
}
template <class... Args>
auto evalSotImpl(Args...) const { assert( false ); return 0; }
template <size_t R, size_t C, class ValueType, class Gradient>
auto evalSotImpl(tag::scalar, MatrixComponent<R, C>, ValueType const& v,
auto evalSotImpl(tag::scalar, MatrixComponent<R, C>, ValueType const& v,
Gradient const& grad_test, Gradient const& grad_trial) const
{
return v * (grad_test[R] * grad_trial[C]);
}
template <class ValueType, class Gradient>
auto evalSotImpl(tag::matrix, tag::none, ValueType const& M,
auto evalSotImpl(tag::matrix, tag::none, ValueType const& M,
Gradient const& grad_test, Gradient const& grad_trial) const
{
return grad_test * (M * grad_trial);
}
template <class ValueType, class Gradient>
auto evalSotImpl(tag::scalar, tag::none, ValueType const& v,
auto evalSotImpl(tag::scalar, tag::none, ValueType const& v,
Gradient const& grad_test, Gradient const& grad_trial) const
{
return v * (grad_test * grad_trial);
}
};
} // end namespace AMDiS
\ No newline at end of file
} // end namespace AMDiS
......@@ -23,7 +23,8 @@ namespace AMDiS
using Codim0 = typename MeshView::template Codim<0>;
using Element = typename Codim0::Entity;
static constexpr int dim = Element::dimension;
static constexpr int dim = MeshView::dimension;
static constexpr int dow = MeshView::dimensionworld;
using QuadratureRule = Dune::QuadratureRule<double, dim>;
using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
......@@ -38,15 +39,15 @@ namespace AMDiS
virtual double evalFot1(size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,dim> const& grad_trial) const = 0;
Dune::FieldVector<double,dow> const& grad_trial) const = 0;
virtual double evalFot2(size_t iq,
Dune::FieldVector<double,dim> const& grad_test,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,1> const trial = 1.0) const = 0;
virtual double evalSot(size_t iq,
Dune::FieldVector<double,dim> const& grad_test,
Dune::FieldVector<double,dim> const& grad_trial) const = 0;
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,dow> const& grad_trial) const = 0;
virtual int getDegree(Dune::GeometryType const& t) const = 0;
};
......@@ -63,7 +64,8 @@ namespace AMDiS
using Element = typename Super::Element;
using PointList = typename Super::PointList;
static constexpr int dim = Element::dimension;
static constexpr int dim = MeshView::dimension;
static constexpr int dow = MeshView::dimensionworld;
public:
GenericOperatorTerm(Term const& term)
......@@ -90,21 +92,21 @@ namespace AMDiS
virtual double evalFot1(size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,dim> const& grad_trial) const override
Dune::FieldVector<double,dow> const& grad_trial) const override
{
return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
}
virtual double evalFot2(size_t iq,
Dune::FieldVector<double,dim> const& grad_test,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,1> const trial = 1.0) const override
{
return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
}
virtual double evalSot(size_t iq,
Dune::FieldVector<double,dim> const& grad_test,
Dune::FieldVector<double,dim> const& grad_trial) const override
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,dow> const& grad_trial) const override
{
return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
}
......
......@@ -26,6 +26,15 @@ namespace AMDiS
template <size_t i, size_t j>
constexpr index_<i+j> operator+(index_<i>, index_<j>) { return {}; }
template <size_t i, size_t j>
constexpr index_<i-j> operator-(index_<i>, index_<j>) { return {}; }
template <size_t i, size_t j>
constexpr index_<i*j> operator*(index_<i>, index_<j>) { return {}; }
template <size_t i, size_t j>
constexpr index_<i/j> operator/(index_<i>, index_<j>) { return {}; }
namespace Impl
{
......
......@@ -10,6 +10,7 @@
using namespace AMDiS;
using PolarParam = DefaultTraitsMesh<Dune::AlbertaGrid<2, 3>, 1, 1, 1, 1>;
using PolarProblem = ProblemStat<PolarParam>;
using PolarProblemInstat = ProblemInstat<PolarParam>;
......@@ -32,7 +33,7 @@ int main(int argc, char** argv)
Parameters::get("polar->K", K);
Parameters::get("polar->R", R);
double H = 1.0/R;
double H = 0.0;
std::srand(std::time(0));
......@@ -51,18 +52,30 @@ int main(int argc, char** argv)
prob.addMatrixOperator(opL, _i, 2+_i);
// <f(p1,p2) * grad(p_i), grad(theta_i)>
auto opL2 = Op::sot( func([wn](WorldVector const& x, WorldVector const& grd_p0, WorldVector const& grd_p1)
auto opL2 = Op::sot( func([wn](WorldVector const& grd_p0, WorldVector const& grd_p1)
{
auto n = x.dot(grd_p0);
auto f = grd_p0; f.axpy(n,x);
f[0] += x[1]*grd_p1[2] - x[2]-grd_p1[1];
f[1] += x[2]*grd_p1[0] - x[0]-grd_p1[2];
f[2] += x[0]*grd_p1[1] - x[1]-grd_p1[0];
return wn * (1.0 - f.dot(f));
}, X<3>(), gradientOf(prob.getSolution(0_c)), gradientOf(prob.getSolution(1_c)) )
return wn * (1.0 - inner(grd_p0, grd_p0) - inner(grd_p1, grd_p1));
}, gradientOf(prob.getSolution(_0)), gradientOf(prob.getSolution(_1)) )
);
prob.addMatrixOperator(opL2, _i, _i);
auto opL2Rhs = Op::fot( func([wn](WorldVector const& grd_pi, WorldVector const& grd_pj)
{
return ((-2.0*wn) * (inner(grd_pi, grd_pi) + inner(grd_pj, grd_pj))) * grd_pi;
}, gradientOf(prob.getSolution(_i)), gradientOf(prob.getSolution(_1 - _i)) ),
GRD_PSI);
prob.addVectorOperator(opL2Rhs, _i);
for_<0,2>([&](auto _j)
{
auto opL3 = Op::sot( func([wn](WorldVector const& grd_pi, WorldVector const& grd_pj)
{
return outer(grd_pj, (-2.0*wn)*grd_pi);
}, gradientOf(prob.getSolution(_i)), gradientOf(prob.getSolution(_j)) )
);
prob.addMatrixOperator(opL3, _i, _j);
});
// <q_i, psi_i>
auto opM = Op::zot( constant(-1.0/K) );
......@@ -70,12 +83,12 @@ int main(int argc, char** argv)
// <-K*grad(p_i), grad(psi_i)>
auto opL3 = Op::sot( constant(1.0) );
prob.addMatrixOperator(opL3, 2+_i, _i);
prob.addMatrixOperator(opL3, _2 + _i, _i);
for (auto& p : prob.getSolution(_i).getVector())
p = 2.0*(std::rand() / double(RAND_MAX)) - 1.0;
for (auto& q : prob.getSolution(2_c + _i).getVector())
for (auto& q : prob.getSolution(_2 + _i).getVector())
q = 0.0;
});
......
......@@ -30,7 +30,7 @@
},
"adapt": {
"timestep": 1.0e-3,
"timestep": 1.0e-4,
"start time": 0.0,
"end time": 1.0
}
......
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