From 633d6024153943d1b57bc7ede566c8797b1bb9bd Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 14 May 2009 09:20:46 +0000
Subject: [PATCH] Simons changes.
---
AMDiS/src/ElementFunction.h | 4 +-
AMDiS/src/Operator.cc | 78 +++++++++++++++++++++++++++
AMDiS/src/Operator.h | 103 ++++++++++++++++++++++++++++++++++++
3 files changed, 183 insertions(+), 2 deletions(-)
diff --git a/AMDiS/src/ElementFunction.h b/AMDiS/src/ElementFunction.h
index b764d7f6..7d4c6b77 100644
--- a/AMDiS/src/ElementFunction.h
+++ b/AMDiS/src/ElementFunction.h
@@ -63,7 +63,7 @@ namespace AMDiS {
{}
/// evaluation at given coordinates.
- const T& operator()(const DimVec<double>& bary) const {
+ T operator()(const DimVec<double>& bary) const {
WorldVector<double> worldCoords;
this->elInfo_->coordToWorld(bary, &worldCoords);
return (*fct_)(worldCoords);
@@ -87,7 +87,7 @@ namespace AMDiS {
{}
/// evaluation at given coordinates.
- const T& operator()(const DimVec<double>& bary) const
+ T operator()(const DimVec<double>& bary) const
{
static T t;
const T* localVec =
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 79ca035a..910ea062 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -2982,4 +2982,82 @@ namespace AMDiS {
}
}
+ void VecGrad_SOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ int dow = Global::getGeo(WORLD);
+
+ if (D2UhAtQP) {
+ for (int iq = 0; iq < numPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ double resultQP = 0.0;
+ for (int i = 0; i < dow; i++) {
+ resultQP += D2UhAtQP[iq][i][i];
+ }
+ result[iq] += fac * resultQP * factor;
+ }
+ }
+ }
+
+ void VecGrad_SOT::weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const
+ {
+ if (grdUhAtQP) {
+ for (int iq = 0; iq < numPoints; iq++) {
+ double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
+ }
+ }
+ }
+
+ void VecGrad_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ for (int iq = 0; iq < numPoints; iq++) {
+ l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
+ }
+ }
+
+ void VecGrad_SOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
+ }
+
+ void FctGrad2_FOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ grad1AtQPs = subAssembler->getGradientsAtQPs(vec1, elInfo, quad);
+ grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
+ }
+
+ void FctGrad2_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ for (int iq = 0; iq < numPoints; iq++) {
+ lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
+ }
+ }
+
+ void FctGrad2_FOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const
+ {
+ if (grdUhAtQP) {
+ for (int iq = 0; iq < numPoints; iq++) {
+ WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
+ result[iq] += b * grdUhAtQP[iq] * factor;
+ }
+ }
+ }
+
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 6378e98b..ae1c4b1f 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -1870,6 +1870,56 @@ namespace AMDiS {
BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
};
+
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * First order term: \f$ b(\nabla v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
+ */
+ class FctGrad2_FOT : public FirstOrderTerm
+ {
+ public:
+ /// Constructor
+ FctGrad2_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct_)
+ : FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_)
+ {}
+
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /// Implements FirstOrderTerm::getLb().
+ void getLb(const ElInfo *elInfo, int qPoint,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
+
+ /// Implements FirstOrderTerm::eval().
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ protected:
+ /// DOFVector to be evaluated at quadrature points.
+ DOFVectorBase<double>* vec1;
+
+ /// DOFVector to be evaluated at quadrature points.
+ DOFVectorBase<double>* vec2;
+
+ /// Gradient v at quadrature points.
+ WorldVector<double> *grad1AtQPs;
+
+ /// Gradient v at quadrature points.
+ WorldVector<double> *grad2AtQPs;
+
+ /// Function for b.
+ BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct;
+ };
+
+
class General_FOT : public FirstOrderTerm
{
public:
@@ -2657,6 +2707,59 @@ namespace AMDiS {
BinaryAbstractFunction<double, double, WorldVector<double> > *f;
};
+
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * Second order term: \f$ f(v(\vec{x}),\nabla w(\vec{x})) \Delta u(\vec{x})\f$
+ */
+ class VecGrad_SOT : public SecondOrderTerm
+ {
+ public:
+ /// Constructor.
+ VecGrad_SOT( DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
+ : SecondOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_)
+ {}
+
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /// Implements SecondOrderTerm::getLALt().
+ inline void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
+
+ /// Implements SecondOrderTerm::eval().
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ /// Implements SecondOrderTerm::weakEval().
+ void weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const;
+
+ protected:
+ /// DOFVectors to be evaluated at quadrature points.
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /// Vector v at quadrature points.
+ double *vecAtQPs;
+
+ /// Gradient at quadrature points.
+ WorldVector<double> *gradAtQPs;
+
+ /// Function for c.
+ BinaryAbstractFunction<double, double, WorldVector<double> > *f;
+ };
+
+
/*
* \ingroup Assembler
*
--
GitLab