Skip to content

GitLab

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

Expressions for gradientOf(vector-values-DOFVector) and trans(matrix) added

parent 8b3e2168
Hide whitespace changes
Inline Side-by-side
Showing with 161 additions and 176 deletions
AMDiS/src/Expressions.h
View file @ c22a8b09
......@@ -94,42 +94,75 @@
namespace AMDiS {
// _______ ZeroOrderTerms ______________________________________________________
template<int Order>
struct GetTerm {
typedef typename boost::mpl::if_c<Order == 0, ZeroOrderTerm,
typename boost::mpl::if_c<Order == 1, FirstOrderTerm,
typename boost::mpl::if_c<Order == 2, SecondOrderTerm,
OperatorTerm
>::type >::type >::type type;
};
template<typename Term>
struct GenericZeroOrderTerm : public ZeroOrderTerm
template<typename Term, int Order = -1>
struct GenericOperatorTerm : public GetTerm<Order>::type
{
typedef typename GetTerm<Order>::type super;
Term term;
GenericZeroOrderTerm(const Term& term_)
: ZeroOrderTerm(term_.getDegree()), term(term_)
GenericOperatorTerm(const Term& term_)
: super(term_.getDegree()), term(term_)
{
term.insertFeSpaces(auxFeSpaces);
term.insertFeSpaces(this->auxFeSpaces);
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad, NULL);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad, NULL);
}
};
template<typename Term>
struct GenericOperatorTerm<Term, -1> : public GenericOperatorTerm<Term, -2>
{
typedef GenericOperatorTerm<Term, -2> super;
GenericOperatorTerm(const Term& term_) : super(term_) { }
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double factor) {};
};
// _______ ZeroOrderTerms ______________________________________________________
template<typename Term>
struct GenericZeroOrderTerm : public GenericOperatorTerm<Term, 0>
{
GenericZeroOrderTerm(const Term& term_)
: GenericOperatorTerm<Term,0>(term_)
{ }
void getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += term(iq);
C[iq] += this->term(iq);
}
void eval( int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
......@@ -138,7 +171,7 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * term(iq) * uhAtQP[iq];
result[iq] += fac * this->term(iq) * uhAtQP[iq];
}
};
......@@ -146,31 +179,11 @@ struct GenericZeroOrderTerm : public ZeroOrderTerm
template<typename Term>
struct GenericFirstOrderTerm_1 : public FirstOrderTerm
struct GenericFirstOrderTerm_1 : public GenericOperatorTerm<Term, 1>
{
Term term;
GenericFirstOrderTerm_1(const Term& term_)
: FirstOrderTerm(term_.getDegree()), term(term_)
{
term.insertFeSpaces(auxFeSpaces);
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
: GenericOperatorTerm<Term, 1>(term_)
{ }
/// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo,
......@@ -180,47 +193,27 @@ struct GenericFirstOrderTerm_1 : public FirstOrderTerm
const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], term(iq));
this->l1(grdLambda, Lb[iq], this->term(iq));
}
};
template<int I, typename Term>
struct GenericFirstOrderTerm_i : public FirstOrderTerm
struct GenericFirstOrderTerm_i : public GenericOperatorTerm<Term, 1>
{
Term term;
GenericFirstOrderTerm_i(const Term& term_)
: FirstOrderTerm(term_.getDegree()), term(term_)
: GenericOperatorTerm<Term, 1>(term_)
{
term.insertFeSpaces(auxFeSpaces);
FirstOrderTerm::bOne = I;
}
GenericFirstOrderTerm_i(const Term& term_, int I0)
: FirstOrderTerm(term_.getDegree()), term(term_)
: GenericOperatorTerm<Term, 1>(term_)
{
term.insertFeSpaces(auxFeSpaces);
FirstOrderTerm::bOne = I0;
TEST_EXIT_DBG( I < 0 && I0 >= 0 )("You yould specify eather template<int I>, or constructor(int I0)\n");
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
/// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo,
std::vector<mtl::dense_vector<double> >& Lb) const
......@@ -229,38 +222,17 @@ struct GenericFirstOrderTerm_i : public FirstOrderTerm
const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], term(iq));
this->lb_one(grdLambda, Lb[iq], this->term(iq));
}
};
template<typename Term>
struct GenericFirstOrderTerm_b : public FirstOrderTerm
struct GenericFirstOrderTerm_b : public GenericOperatorTerm<Term, 1>
{
Term term;
GenericFirstOrderTerm_b(const Term& term_)
: FirstOrderTerm(term_.getDegree()), term(term_)
{
term.insertFeSpaces(auxFeSpaces);
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
: GenericOperatorTerm<Term, 1>(term_)
{ }
/// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo,
......@@ -270,7 +242,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, term(iq), Lb[iq], 1.0);
this->lb(grdLambda, this->term(iq), Lb[iq], 1.0);
}
};
......@@ -278,31 +250,12 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
template<typename Term>
struct GenericSecondOrderTerm_1 : public SecondOrderTerm
struct GenericSecondOrderTerm_1 : public GenericOperatorTerm<Term, 2>
{
Term term;
GenericSecondOrderTerm_1(const Term& term_)
: SecondOrderTerm(term_.getDegree()), term(term_)
{
term.insertFeSpaces(auxFeSpaces);
setSymmetric(true);
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
: GenericOperatorTerm<Term, 2>(term_)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
this->setSymmetric(true);
}
/// Implements SecondOrderTerm::getLALt().
......@@ -312,7 +265,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
const int nPoints = static_cast<int>(LALt.size());
for (int iq = 0; iq < nPoints; iq++)
l1lt(grdLambda, LALt[iq], term(iq));
this->l1lt(grdLambda, LALt[iq], this->term(iq));
}
/// Implemetation of SecondOrderTerm::eval().
......@@ -331,7 +284,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
for (int i = 0; i < dow; i++) {
resultQP += D2UhAtQP[iq][i][i];
}
result[iq] += resultQP * f * term(iq);
result[iq] += resultQP * f * this->term(iq);
}
}
}
......@@ -342,38 +295,19 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
axpy(term(iq), grdUhAtQP[iq], result[iq]);
axpy(this->term(iq), grdUhAtQP[iq], result[iq]);
}
};
template<typename Term, bool symmetric = false>
struct GenericSecondOrderTerm_A : public SecondOrderTerm
{
Term term;
struct GenericSecondOrderTerm_A : public GenericOperatorTerm<Term, 2>
{
GenericSecondOrderTerm_A(const Term& term_)
: SecondOrderTerm(term_.getDegree()), term(term_)
{
term.insertFeSpaces(auxFeSpaces);
setSymmetric(symmetric);
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
: GenericOperatorTerm<Term, 2>(term_)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
this->setSymmetric(symmetric);
}
void getLALt(const ElInfo *elInfo,
......@@ -383,7 +317,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
const int nPoints = static_cast<int>(LALt.size());
for (int iq = 0; iq < nPoints; iq++)
lalt(grdLambda, term(iq), LALt[iq], symmetric, 1.0);
this->lalt(grdLambda, this->term(iq), LALt[iq], symmetric, 1.0);
}
void eval(int nPoints,
......@@ -398,7 +332,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
WorldMatrix<double> A = term(iq);
WorldMatrix<double> A = this->term(iq);
if (num_rows(D2UhAtQP) > 0)
for (int i = 0; i < dow; i++)
......@@ -418,8 +352,7 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = term(iq);
result[iq] += term(iq) * grdUhAtQP[iq];
result[iq] += this->term(iq) * grdUhAtQP[iq];
}
}
};
......@@ -427,43 +360,24 @@ struct GenericSecondOrderTerm_A : public SecondOrderTerm
template<int I, int J, typename Term>
struct GenericSecondOrderTerm_ij : public SecondOrderTerm
struct GenericSecondOrderTerm_ij : public GenericOperatorTerm<Term, 2>
{
Term term;
int row, col;
GenericSecondOrderTerm_ij(const Term& term_)
: SecondOrderTerm(term_.getDegree()), term(term_), row(I), col(J)
: GenericOperatorTerm<Term, 2>(term_), row(I), col(J)
{
term.insertFeSpaces(auxFeSpaces);
setSymmetric(row == col);
this->setSymmetric(row == col);
}
GenericSecondOrderTerm_ij(const Term& term_, int I0, int J0)
: SecondOrderTerm(term_.getDegree()), term(term_), row(I0), col(J0)
: GenericOperatorTerm<Term, 2>(term_), row(I0), col(J0)
{
term.insertFeSpaces(auxFeSpaces);
setSymmetric(row == col);
this->setSymmetric(row == col);
TEST_EXIT_DBG( I < 0 && I0 >= 0 && J < 0 && J0 >= 0 )
("You yould specify eather template<int I, int J>, or constructor(int I0, int J0)\n");
}
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, elInfo, subAssembler, quad);
}
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
}
void getLALt(const ElInfo *elInfo,
std::vector<mtl::dense2D<double> > &LALt) const
{
......@@ -471,7 +385,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
const int nPoints = static_cast<int>(LALt.size());
for (int iq = 0; iq < nPoints; iq++)
lalt_kl(grdLambda, row, col, LALt[iq], term(iq));
this->lalt_kl(grdLambda, row, col, LALt[iq], this->term(iq));
}
void eval(int nPoints,
......@@ -483,7 +397,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
{
if (num_rows(D2UhAtQP) > 0) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += D2UhAtQP[iq][row][col] * term(iq) * fac;
result[iq] += D2UhAtQP[iq][row][col] * this->term(iq) * fac;
}
}
......@@ -492,7 +406,7 @@ struct GenericSecondOrderTerm_ij : public SecondOrderTerm
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
result[iq][row] += grdUhAtQP[iq][col] * term(iq);
result[iq][row] += grdUhAtQP[iq][col] * this->term(iq);
}
};
......
AMDiS/src/Expressions.hh
View file @ c22a8b09
......@@ -24,6 +24,7 @@
namespace AMDiS {
#if 0
namespace detail {
template<typename Term, typename T, typename Enable = void>
......@@ -50,6 +51,7 @@ template<typename Term>
struct GenericOperatorTerm {
typedef typename detail::GenericOperatorTerm<Term, typename Term::value_type>::type type;
};
#endif
template<typename Term>
......@@ -58,7 +60,7 @@ integrate(Term term, Mesh* mesh_)
{
typedef typename Term::value_type TOut;
typename GenericOperatorTerm<Term>::type ot(term);
GenericOperatorTerm<Term> ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
TEST_EXIT(mesh_ || !feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
......@@ -96,7 +98,7 @@ template<typename Term, typename Functor>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
accumulate(Term term, Functor f, typename Term::value_type value0)
{
typename GenericOperatorTerm<Term>::type ot(term);
GenericOperatorTerm<Term> ot(term);
std::set<const FiniteElemSpace*> feSpaces = ot.getAuxFeSpaces();
TEST_EXIT(!feSpaces.empty())("The Expression must contain a DOFVector or FeSpace depended value!\n");
......@@ -148,7 +150,7 @@ transformDOF(Term term, DOFVector<T>* result)
DOFVector<TOut> temp(result->getFeSpace(), "temp");
DOFVector<int> assigned(result->getFeSpace(), "assigned");
typename GenericOperatorTerm<Term>::type ot(term);
GenericOperatorTerm<Term> ot(term);
Mesh* mesh = result->getFeSpace()->getMesh();
const FiniteElemSpace* resultFeSpace = temp.getFeSpace();
......
AMDiS/src/expressions/vec_functors.hpp
View file @ c22a8b09
......@@ -422,6 +422,35 @@ namespace AMDiS
return function_(expressions::MatComponent<typename Term::value_type>(I, J), t);
}
// transpose a matrix
// _____________________________________________________________________________
namespace expressions
{
template<typename Mat>
struct MatTranspose : public FunctorBase
{
typedef Mat result_type;
int getDegree(int d0) const { return d0; }
result_type operator()(const Mat &m) const
{
Mat result;
for (size_t r = 0; r < num_rows(m); r++)
for (size_t c = 0; c < num_cols(m); c++)
result[c][r] = m[r][c];
return result;
}
};
}
template<typename Term>
typename result_of::UnaryExpr<expressions::MatTranspose, Term, traits::is_matrix>::type
trans(const Term& t)
{
return function_(expressions::MatTranspose<typename Term::value_type>(), t);
}
}
#endif // AMDIS_VEC_FUNCTORS_HPP
AMDiS/src/operations/product.hpp
View file @ c22a8b09
......@@ -237,22 +237,46 @@ namespace AMDiS
template <typename T>
struct outer< WorldVector<T>, WorldVector<T>, tag::vector, tag::vector >
: boost::mpl::identity<WorldMatrix<T> > {};
template <typename Vector, typename Scalar>
struct outer< Vector, Scalar, tag::vector, tag::scalar >
: boost::mpl::identity<Vector> {};
template <typename Scalar, typename Vector>
struct outer< Scalar, Vector, tag::scalar, tag::vector >
: boost::mpl::identity<Vector> {};
}
template <typename T>
WorldMatrix<T>
inline outer_dispatch(const WorldVector<T>& v1, const WorldVector<T>& v2, boost::mpl::false_)
inline outer_dispatch(const WorldVector<T>& v1, const WorldVector<T>& v2, tag::vector, tag::vector)
{
WorldMatrix<T> result;
result.vecProduct(v1, v2);
return result;
}
template <typename T>
WorldMatrix<T>
inline outer(const WorldVector<T>& v1, const WorldVector<T>& v2)
template <typename Vector, typename Scalar>
Vector inline outer_dispatch(const Vector& v, const Scalar& s, tag::vector, tag::scalar)
{
Vector result(v);
result *= s;
return result;
}
template <typename Scalar, typename Vector>
Vector inline outer_dispatch(const Scalar& s, const Vector& v, tag::scalar, tag::vector)
{
return outer_dispatch(v1, v2, boost::mpl::false_());
return outer_dispatch(v, s, tag::vector(), tag::scalar());
}
template <typename T1, typename T2>
typename result_of::outer<T1, T2, typename traits::category<T1>::tag, typename traits::category<T2>::tag>::type
inline outer(const T1& v1, const T2& v2)
{
typedef typename traits::category<T1>::tag tag1_;
typedef typename traits::category<T2>::tag tag2_;
return outer_dispatch(v1, v2, tag1_(), tag2_());
}
......
AMDiS/src/parallel/PetscSolver.cc
View file @ c22a8b09
......@@ -103,7 +103,7 @@ namespace AMDiS { namespace Parallel {
destroyMatrixData();
}
return 0;
return getErrorCode();
}
void PetscSolver::fillPetscMatrix(DOFMatrix* mat)
......
AMDiS/src/solver/LinearSolverInterface.h
View file @ c22a8b09
......@@ -126,6 +126,14 @@ namespace AMDiS {
} else if (residual >= 0.0) {
MSG("Residual norm: ||b-Ax|| = %e\n", residual);
}
if (getIterations() > 0) {
MSG("Nr. of iterations needed = %d\n", getIterations());
}
if (error_code != 0) {
MSG("ERROR-Code = %d\n", error_code);
}
// test for absolute tolerance
TEST_EXIT((isNumber(residual) && (residual < 0.0 || tolerance < 1.e-30 || residual <= tolerance))
......
AMDiS/src/traits/category.hpp
View file @ c22a8b09
......@@ -167,6 +167,14 @@ namespace AMDiS
typedef int size_type;
};
template<typename T>
struct category<WorldVector<WorldVector<T> > >
{
typedef tag::matrix tag;
typedef T value_type;
typedef int size_type;
};
template<typename T, typename Parameters>
struct category<mtl::dense2D<T, Parameters> >
{
......