diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index a9e695f74272e34dfef84dcbde37cf71a21653d7..711aec0a60874c3ea0ffb96b98e906b08a5dfe67 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -275,7 +275,8 @@ namespace AMDiS {
}
ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat,
- const ElInfo *elInfo)
+ const ElInfo *rowElInfo,
+ const ElInfo *colElInfo)
{
if (!elMat) {
elMat = NEW ElementMatrix(nRow, nCol);
@@ -283,17 +284,23 @@ namespace AMDiS {
elMat->set(0.0);
- rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
+ rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
rowFESpace->getAdmin(),
&(elMat->rowIndices));
-
if (rowFESpace == colFESpace) {
elMat->colIndices = elMat->rowIndices;
} else {
- colFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
- colFESpace->getAdmin(),
- &(elMat->colIndices));
+ if (colElInfo) {
+ colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
+ colFESpace->getAdmin(),
+ &(elMat->colIndices));
+ } else {
+ // If there is no colElInfo pointer, use rowElInfo the get the indices.
+ colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
+ colFESpace->getAdmin(),
+ &(elMat->colIndices));
+ }
}
return elMat;
@@ -309,7 +316,6 @@ namespace AMDiS {
elVec->set(0.0);
DOFAdmin *rowAdmin = rowFESpace->getAdmin();
-
Element *element = elInfo->getElement();
rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index f097a8b83ae0cc1b5d48cb32ec7adb3536dda6ac..2ea2f69a7c31bea27d2680240d84bbfa3424eaf3 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -70,7 +70,8 @@ namespace AMDiS {
virtual ~Assembler() {};
ElementMatrix *initElementMatrix(ElementMatrix *elMat,
- const ElInfo *elInfo);
+ const ElInfo *rowElInfo,
+ const ElInfo *colElInfo = NULL);
ElementVector *initElementVector(ElementVector *elVec,
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 6f705d1ba0badbe7050e88d91b78e3e10dd5b825..bb255e2cfe35b8133328150bab438d4ce95effc1 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -245,7 +245,6 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) { // for all rows of element matrix
row = elMat.rowIndices[i];
-
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
@@ -402,30 +401,28 @@ namespace AMDiS {
void DOFMatrix::freeDOFContent(int index)
{
- int i, j, col=0, col2;
+ int col = 0;
if (0 < matrix[index].size()) {
// for all columns in this row
int size = static_cast<int>(matrix[index].size());
- for (i=0; i<size; i++) {
+ for (int i = 0; i < size; i++) {
// if entry is used
- if (entryUsed(index,i)) {
+ if (entryUsed(index, i)) {
// get column of this entry
col = matrix[index][i].col;
if (col != index) { // remove symmetric entry if exists
int colsize = static_cast<int>(matrix[col].size());
- for (j=0; j< colsize; j++) {
- col2 = matrix[col][j].col;
+ for (int j = 0; j < colsize; j++) {
+ int col2 = matrix[col][j].col;
if (col2 == index) {
matrix[col][j].col = DOFMatrix::UNUSED_ENTRY;
- }
- else if (col2 == DOFMatrix::NO_MORE_ENTRIES) {
+ } else if (col2 == DOFMatrix::NO_MORE_ENTRIES) {
break;
}
}
}
- }
- else if (col == DOFMatrix::NO_MORE_ENTRIES) {
+ } else if (col == DOFMatrix::NO_MORE_ENTRIES) {
break;
}
}
@@ -473,7 +470,8 @@ namespace AMDiS {
}
Operator *operat = op ? op : operators[0];
- operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, rowElInfo);
+ operat->getAssembler(omp_get_thread_num())->
+ initElementMatrix(elementMatrix, rowElInfo, colElInfo);
if (op) {
op->getElementMatrix(rowElInfo, colElInfo, elementMatrix);
@@ -488,7 +486,7 @@ namespace AMDiS {
*factorIt ? **factorIt : 1.0);
}
}
-
+
addElementMatrix(factor, *elementMatrix, bound);
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 9858c2df926f9ed40dc4339c5f888b5061fbec08..4e068a33e5adbe00a43e65c411da3168c65e6c04 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -655,7 +655,9 @@ namespace AMDiS {
}
- DOFVectorDOF::DOFVectorDOF() : DOFVector<DegreeOfFreedom>() {};
+ DOFVectorDOF::DOFVectorDOF() :
+ DOFVector<DegreeOfFreedom>()
+ {};
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 572887ab594191b3ee3f7210ab03ca317de9eae0..aca25f463960f579e6b6cd6b4c5d805509b567b2 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -316,8 +316,8 @@ namespace AMDiS {
*/
DOFVector()
: DOFVectorBase<T>(),
- refineInter(false),
feSpace(NULL),
+ refineInter(false),
coarsenOperation(NO_OPERATION)
{};
@@ -531,6 +531,11 @@ namespace AMDiS {
inline double l1norm() const {
return asum();
};
+
+ /** \brief
+ * Calculates doublewell of this DOFVector
+ */
+ double DoubleWell(Quadrature* q = NULL) const;
/** \brief
* Calculates the sum of this DOFVector
@@ -637,6 +642,11 @@ namespace AMDiS {
*/
static int H1Norm_fct(ElInfo* elinfo);
+ /** \brief
+ * Used by DoubleWell while mesh traversal
+ */
+ static int DoubleWell_fct(ElInfo* elinfo);
+
protected:
/** \brief
* Name of this DOFVector
@@ -859,6 +869,15 @@ namespace AMDiS {
const DOFVector<T>& v2,
DOFVector<T>& result);
+ template<typename T>
+ inline const DOFVector<T>& mod(const DOFVector<T>& v,
+ DOFVector<T>& result);
+
+ template<typename T>
+ inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
+ DOFVector<T>& result);
+
+
template<typename T>
inline void set(DOFVector<T>& vec, T d)
{
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index bc7179a9785d05c8f1a78e52eba961185dfc25a6..3fe022e434ffc7cf7a78d3f12b102579899546db 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -1,5 +1,6 @@
#include <list>
#include <algorithm>
+#include <math.h>
#include "FixVec.h"
#include "Boundary.h"
@@ -16,6 +17,7 @@
#include "Assembler.h"
#include "OpenMP.h"
#include "Operator.h"
+#include "Parameters.h"
namespace AMDiS {
@@ -1161,9 +1163,43 @@ namespace AMDiS {
*rIterator = (*v1Iterator) + (*v2Iterator);
};
return result;
+ }
+
+ template<typename T>
+ inline const DOFVector<T>& mod(const DOFVector<T>& v,
+ DOFVector<T>& result)
+ {
+ typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
+ typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
+ for (vIterator.reset(), rIterator.reset();
+ !vIterator.end();
+ ++vIterator, ++rIterator) {
+ *rIterator = fmod((*vIterator), 1.0);
+ }
+
+ return result;
+ }
+
+ template<typename T>
+ inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
+ DOFVector<T>& result)
+ {
+ typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
+ typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
+
+ double eps;
+ GET_PARAMETER(1, "vecphase->epsilon", "%f", &eps);
+ for (vIterator.reset(), rIterator.reset();
+ !vIterator.end();
+ ++vIterator, ++rIterator) {
+ *rIterator = 0.5 * (1.0 - tanh(3.0 * (*vIterator) / eps));
+ }
+
+ return result;
}
+
template<typename T>
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{
@@ -1248,4 +1284,47 @@ namespace AMDiS {
return const_cast<const T*>(result);
}
+ // integral u^2(1-u)^2
+ template<typename T>
+ double DOFVector<T>::DoubleWell(Quadrature* q) const
+ {
+ FUNCNAME("DOFVector::DoubleWell()");
+
+ Mesh* mesh = feSpace->getMesh();
+
+ int deg = 0;
+ dim = mesh->getDim();
+
+ if (!q) {
+ deg = 2 * feSpace->getBasisFcts()->getDegree();
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+
+ quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
+ *q, INIT_PHI);
+
+ norm = 0.0;
+ traverseVector = const_cast<DOFVector<T>*>(this);
+
+ mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, DoubleWell_fct);
+
+ return norm;
+ }
+
+ template<typename T>
+ int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
+ {
+ double det = elinfo->getDet();
+ const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL,
+ quad_fast, NULL);
+ int numPoints = quad_fast->getNumPoints();
+ double normT = 0.0;
+ for (int iq = 0; iq < numPoints; iq++) {
+ normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq])*sqr(1.0 - uh_vec[iq]);
+ }
+ norm += det * normT;
+
+ return 0;
+ }
+
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index d2ff25732e4ee565b0ea22460e538f7515408631..a10e8c12c4d4dc9edb1e59c4089aa86c83557289 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -251,7 +251,7 @@ namespace AMDiS {
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0)
- {
+ {
if (!assembler[rank])
if (optimized) {
assembler[rank] =
@@ -280,6 +280,14 @@ namespace AMDiS {
{
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
}
+
+ void Vec2AtQP_SOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ }
void CoordsAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
@@ -310,6 +318,14 @@ namespace AMDiS {
gradAtQPs = subAssembler->getGradientsAtQPs(vec, elInfo, quad);
}
+ void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
+ gradAtQPs = subAssembler->getGradientsAtQPs(vec, elInfo, quad);
+ }
+
void VecAndCoordsAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -400,6 +416,15 @@ namespace AMDiS {
vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
}
+ void Vec3AtQP_ZOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ vecAtQPs3 = subAssembler->getVectorAtQPs(vec3, elInfo, quad);
+ }
+
void VecAndCoordsAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -407,7 +432,16 @@ namespace AMDiS {
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
-
+
+ void Vec2AndGradAtQP_ZOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
+ vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
+ gradAtQPs = subAssembler->getGradientsAtQPs(vec1, elInfo, quad);
+ }
+
void FctGradientCoords_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -476,7 +510,14 @@ namespace AMDiS {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
}
-
+
+ void Vec2AtQP_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)(vecAtQPs1[iq], vecAtQPs2[iq]));
+ }
+ }
+
void CoordsAtQP_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++) {
@@ -506,6 +547,14 @@ namespace AMDiS {
}
}
+ void VecMatrixGradientAtQP_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++) {
+ lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ }
+ }
+
void VecAndCoordsAtQP_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++) {
@@ -644,6 +693,28 @@ namespace AMDiS {
}
}
+ void Vec3AtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
+ for (int iq = 0; iq < numPoints; iq++) {
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
+ }
+ }
+
+ void Vec3AtQP_ZOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ for (int iq = 0; iq < numPoints; iq++) {
+ result[iq] +=
+ fac *
+ (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) *
+ uhAtQP[iq];
+ }
+ }
+
+
void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
@@ -704,6 +775,28 @@ namespace AMDiS {
}
}
+ void Vec2AndGradAtQP_ZOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const
+ {
+ for (int iq = 0; iq < numPoints; iq++) {
+ result[iq] +=
+ fac *
+ (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]) *
+ uhAtQP[iq];
+ }
+ }
+
+ void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const
+ {
+ for (int iq = 0; iq < numPoints; iq++) {
+ C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
+ }
+ }
+
void VecAndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -959,6 +1052,39 @@ namespace AMDiS {
}
}
+ void Vec2AtQP_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)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ double resultQP = 0.0;
+ for (int i = 0; i < dow; i++) {
+ resultQP += D2UhAtQP[iq][i][i];
+ }
+ result[iq] += fac * factor * resultQP;
+ }
+ }
+ }
+
+ void Vec2AtQP_SOT::weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const
+ {
+ if (grdUhAtQP) {
+ for (int iq = 0; iq < numPoints; iq++) {
+ double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ axpy(factor, grdUhAtQP[iq], result[iq]);
+ }
+ }
+ }
+
void CoordsAtQP_SOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1058,7 +1184,48 @@ namespace AMDiS {
}
}
+ void VecMatrixGradientAtQP_SOT::eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const
+ {
+ int dow = Global::getGeo(WORLD);
+
+ for (int iq = 0; iq < numPoints; iq++) {
+ double resultQP = 0.0;
+
+ WorldMatrix<double> A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+
+ if (D2UhAtQP) {
+ for (int i = 0; i < dow; i++) {
+ for (int j = 0; j < dow; j++) {
+ resultQP += A[i][j] * D2UhAtQP[iq][j][i];
+ }
+ }
+ }
+
+ if (grdUhAtQP) {
+ resultQP += (*divFct)(A) * grdUhAtQP[iq];
+ }
+
+ result[iq] += resultQP * factor;
+ }
+ }
+ void VecMatrixGradientAtQP_SOT::weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const
+ {
+ if (grdUhAtQP) {
+ WorldMatrix<double> A;
+ for (int iq = 0; iq < numPoints; iq++) {
+ A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
+ result[iq] += A * grdUhAtQP[iq];
+ }
+ }
+ }
void VecAndCoordsAtQP_SOT::eval(int numPoints,
const double *uhAtQP,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 28cc15c31e02c65d2f72cb0e6c74203ed9d69c58..73540d4cbadee09efea8a719d2ab533443b03ea2 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -718,6 +718,77 @@ namespace AMDiS {
AbstractFunction<double, double> *f;
};
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * Laplace operator multiplied with a function evaluated at the quadrature
+ * points of a given DOFVector:
+ * \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$
+ */
+ class Vec2AtQP_SOT : public SecondOrderTerm {
+ public:
+ /** \brief
+ * Constructor.
+ */
+ Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *f_)
+ : SecondOrderTerm(f_->getDegree()),
+ vec1(dv1),
+ vec2(dv2),
+ f(f_)
+ {
+ setSymmetric(true);
+ };
+
+ /** \brief
+ * Implementation of \ref OperatorTerm::initElement().
+ */
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /** \brief
+ * Implements SecondOrderTerm::getLALt().
+ */
+ void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
+
+ /** \brief
+ * Implenetation of SecondOrderTerm::eval().
+ */
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ /** \brief
+ * Implenetation of SecondOrderTerm::weakEval().
+ */
+ void weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const;
+
+ protected:
+ /** \brief
+ * DOFVector to be evaluated at quadrature points.
+ */
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /** \brief
+ * Pointer to an array containing the DOFVector evaluated at quadrature
+ * points.
+ */
+ double* vecAtQPs1;
+ double* vecAtQPs2;
+
+ /** \brief
+ * Function evaluated at \ref vecAtQPs.
+ */
+ BinaryAbstractFunction<double, double, double> *f;
+ };
+
// ============================================================================
/**
@@ -989,6 +1060,82 @@ namespace AMDiS {
WorldVector<double>* gradAtQPs;
};
+ // ============================================================================
+
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * Laplace multiplied with a function which maps the gradient of a DOFVector
+ * at each quadrature point to a double:
+ * \f$ \nabla \cdot A(v(\vec{x}), \nabla v(\vec{x})) \nabla u(\vec{x}) \f$
+ */
+ class VecMatrixGradientAtQP_SOT : public SecondOrderTerm
+ {
+ public:
+ /** \brief
+ * Constructor.
+ */
+ VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
+ BinaryAbstractFunction<WorldMatrix<double>,
+ double, WorldVector<double> > *f_,
+ AbstractFunction<WorldVector<double>,
+ WorldMatrix<double> > *divFct_,
+ bool symm = false)
+ : SecondOrderTerm(f_->getDegree()), vec(dv), f(f_),
+ divFct(divFct_), symmetric(symm)
+ {
+ setSymmetric(symmetric);
+ };
+
+ /** \brief
+ * Implementation of \ref OperatorTerm::initElement().
+ */
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /** \brief
+ * Implements SecondOrderTerm::getLALt().
+ */
+ void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
+
+ /** \brief
+ * Implenetation of SecondOrderTerm::eval().
+ */
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ /** \brief
+ * Implenetation of SecondOrderTerm::weakEval().
+ */
+ void weakEval(int numPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const;
+
+ protected:
+ DOFVectorBase<double>* vec;
+
+ /** \brief
+ * Function wich maps \ref gradAtQPs to a double.
+ */
+ BinaryAbstractFunction<WorldMatrix<double>,
+ double, WorldVector<double> > *f;
+
+ AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;
+
+ /** \brief
+ * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
+ * at each quadrature point.
+ */
+ double* vecAtQPs;
+ WorldVector<double>* gradAtQPs;
+
+ bool symmetric;
+};
// ============================================================================
@@ -2226,9 +2373,6 @@ namespace AMDiS {
AbstractFunction<double, double> *f;
};
- // Andreas ergaenzt:
-
-
/**
* \ingroup Assembler
*
@@ -2355,6 +2499,72 @@ namespace AMDiS {
BinaryAbstractFunction<double, double, double> *f;
};
+
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * Zero order term: \f$ f(v(\vec{x}), w(\vec{x})) u(\vec{x})\f$
+ */
+ class Vec3AtQP_ZOT : public ZeroOrderTerm
+ {
+ public:
+ /** \brief
+ * Constructor.
+ */
+ Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ DOFVectorBase<double> *dv3,
+ TertiaryAbstractFunction<double, double, double, double> *f_)
+ : ZeroOrderTerm(f_->getDegree()),
+ vec1(dv1),
+ vec2(dv2),
+ vec3(dv3),
+ f(f_)
+ {};
+
+ /** \brief
+ * Implementation of \ref OperatorTerm::initElement().
+ */
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /** \brief
+ * Implements ZeroOrderTerm::getC().
+ */
+ void getC(const ElInfo *, int numPoints, double *C) const;
+
+ /** \brief
+ * Implements ZeroOrderTerm::eval().
+ */
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const;
+
+ protected:
+ /** \brief
+ * DOFVector to be evaluated at quadrature points.
+ */
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+ DOFVectorBase<double>* vec3;
+
+ /** \brief
+ * Vector v at quadrature points.
+ */
+ double *vecAtQPs1;
+ double *vecAtQPs2;
+ double *vecAtQPs3;
+
+ /** \brief
+ * Function for c.
+ */
+ TertiaryAbstractFunction<double, double, double, double> *f;
+ };
+
// ============================================================================
/**
@@ -2542,7 +2752,68 @@ namespace AMDiS {
};
- // Ende Andreas Ergaenzungen
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ */
+ class Vec2AndGradAtQP_ZOT : public ZeroOrderTerm
+ {
+ public:
+ /** \brief
+ * Constructor.
+ */
+ Vec2AndGradAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ TertiaryAbstractFunction<double, double, WorldVector<double>,
+ double > *f_)
+ : ZeroOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_)
+ {
+ };
+
+ /** \brief
+ * Implementation of \ref OperatorTerm::initElement().
+ */
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /** \brief
+ * Implements ZeroOrderTerm::getC().
+ */
+ void getC(const ElInfo *, int numPoints, double *C) const;
+
+ /** \brief
+ * Implements ZeroOrderTerm::eval().
+ */
+ void eval(int numPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double fac) const;
+
+ protected:
+ /** \brief
+ * DOFVector to be evaluated at quadrature points.
+ */
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /** \brief
+ * Vector v at quadrature points.
+ */
+ double *vecAtQPs1;
+ double *vecAtQPs2;
+
+ /** \brief
+ * Gradient at quadrature points.
+ */
+ WorldVector<double> *gradAtQPs;
+
+ /** \brief
+ * Function for c.
+ */
+ TertiaryAbstractFunction<double, double, WorldVector<double>, double > *f;
+ };
// ============================================================================
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index aa3bb9936d5367d42c3624357055cbd73c1d6758..f11a7f1c2a303a62714331fb000920be6ef818b2 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -1004,10 +1004,12 @@ namespace AMDiS {
Element *colElem = colElInfo->getElement();
if (rowElInfo->getLevel() != colElInfo->getLevel()) {
- if (smallElInfo == rowElInfo)
+ if (smallElInfo == rowElInfo) {
std::cout << "Row = small\n";
- else
+ ERROR_EXIT("NOCH EIN PAAR GEDANKEN MACHEN!\n");
+ } else {
std::cout << "Row = large\n";
+ }
if (largeElInfo == colElInfo)
std::cout << "Col = large\n";
@@ -1019,21 +1021,21 @@ namespace AMDiS {
}
if (matrix) {
- matrix->assemble(1.0, colElInfo, rowElInfo, bound);
+ matrix->assemble(1.0, rowElInfo, colElInfo, bound);
+ }
+ } else {
+ if (useGetBound_) {
+ basisFcts->getBound(rowElInfo, bound);
}
- }
-
- if (useGetBound_) {
- basisFcts->getBound(rowElInfo, bound);
- }
- if (matrix) {
- matrix->assemble(1.0, rowElInfo, bound);
-
- if (matrix->getBoundaryManager()) {
- matrix->getBoundaryManager()->
- fillBoundaryConditions(rowElInfo, matrix);
- }
+ if (matrix) {
+ matrix->assemble(1.0, rowElInfo, bound);
+
+ if (matrix->getBoundaryManager()) {
+ matrix->getBoundaryManager()->
+ fillBoundaryConditions(rowElInfo, matrix);
+ }
+ }
}
if (vector) {
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 7ffe128bfe208fdd705f922675d9db292d86f85b..6ff161b385143d76a6b28bf80503daab209412a0 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -193,7 +193,7 @@ namespace AMDiS {
FREE_MEMORY(phival, double, nCol);
FREE_MEMORY(c, double, nPoints);
- ERROR_EXIT("SO, HIER GEHTS WEITER\n");
+ // ERROR_EXIT("SO, HIER GEHTS WEITER\n");
}