diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h
index d24a2c4a890c7cd411c2df55659aef2fddee37b5..ef70750889b0efb5edff8f1f6b68e77e489de72d 100644
--- a/AMDiS/src/MatrixVector.h
+++ b/AMDiS/src/MatrixVector.h
@@ -36,7 +36,8 @@ namespace AMDiS {
public:
/// Constructor.
Vector(int i = 0)
- : size(i)
+ : size(i),
+ valArray(NULL)
{
if (size == 0)
valArray = NULL;
@@ -53,7 +54,7 @@ namespace AMDiS {
inline void resize(int newSize)
{
if (size != newSize) {
- if (valArray)
+ if (valArray != NULL)
delete [] valArray;
valArray = new T[newSize];
size = newSize;
@@ -71,7 +72,10 @@ namespace AMDiS {
/// Destructor.
virtual ~Vector()
{
- delete [] valArray;
+ if (valArray != NULL) {
+ delete [] valArray;
+ valArray = NULL;
+ }
}
/// Assignement operator
@@ -81,9 +85,8 @@ namespace AMDiS {
T *rhsIt, *thisIt;
for (rhsIt = rhs.begin(), thisIt = this->begin();
rhsIt != rhs.end();
- ++rhsIt, ++thisIt) {
- *thisIt = *rhsIt;
- }
+ ++rhsIt, ++thisIt)
+ *thisIt = *rhsIt;
return *this;
}
@@ -91,10 +94,7 @@ namespace AMDiS {
/// Assignement operator
inline const Vector<T>& operator=(const T& scal)
{
- T *thisIt;
- for (thisIt = this->begin();
- thisIt != this->end();
- ++thisIt)
+ for (T *thisIt = this->begin(); thisIt != this->end(); ++thisIt)
*thisIt = scal;
return *this;
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 1c8658d6b81d06528786ecf5178386b3a7250aab..f5d4f6ccb031b89733b4698dc2e4e22c45b727b9 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -183,5 +183,16 @@ namespace AMDiS {
{
assembler[rank] = a;
}
+
+ void Operator::weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result) const
+ {
+ int myRank = omp_get_thread_num();
+
+ for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin();
+ termIt != secondOrder[myRank].end();
+ ++termIt)
+ static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
+ }
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index bff990b68963c86f11394254683bba0a22587662..d0bba2771635ac5e98f78705ced1c66e84e7a6fb 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -280,15 +280,7 @@ namespace AMDiS {
/// Weak evaluation of all terms in \ref secondOrder.
void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const
- {
- int myRank = omp_get_thread_num();
-
- for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin();
- termIt != secondOrder[myRank].end();
- ++termIt)
- static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
- }
+ std::vector<WorldVector<double> > &result) const;
/// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 9df83ac1838735235480daf5b74b8970979b2f93..e1c42b352e658a7e10699d23bbe012a5896ac4c8 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -61,6 +61,7 @@ namespace AMDiS {
}
}
+
RobinBC::RobinBC(BoundaryType type,
DOFVectorBase<double> *j,
DOFVectorBase<double> *alpha,
@@ -114,6 +115,7 @@ namespace AMDiS {
}
}
+
RobinBC::RobinBC(BoundaryType type,
Operator* jOp, Operator* alphaOp,
FiniteElemSpace *rowFESpace_,
@@ -152,6 +154,7 @@ namespace AMDiS {
}
}
+
void RobinBC::fillBoundaryCondition(DOFVectorBase<double>* vector,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
@@ -163,15 +166,13 @@ namespace AMDiS {
int dim = elInfo->getMesh()->getDim();
- if (neumannOperators) {
- for (int i = 0; i < dim + 1; i++) {
- if (elInfo->getBoundary(i) == boundaryType) {
+ if (neumannOperators)
+ for (int i = 0; i < dim + 1; i++)
+ if (elInfo->getBoundary(i) == boundaryType)
vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
- }
- }
- }
}
+
void RobinBC::fillBoundaryCondition(DOFMatrix* matrix,
ElInfo* elInfo,
const DegreeOfFreedom* dofIndices,
@@ -187,6 +188,7 @@ namespace AMDiS {
}
}
+
double RobinBC::boundResidual(ElInfo *elInfo,
DOFMatrix *matrix,
const DOFVectorBase<double> *dv)
@@ -196,13 +198,11 @@ namespace AMDiS {
TEST_EXIT(matrix->getColFESpace() == colFESpace)("invalid col fe space\n");
int dim = elInfo->getMesh()->getDim();
- int iq;
DimVec<double> lambda(dim, NO_INIT);
double n_A_grdUh, val = 0.0;
WorldVector<double> normal;
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
double det = elInfo->getDet();
- int numPoints;
bool neumannQuad = false;
const BasisFunction *basFcts = dv->getFESpace()->getBasisFcts();
@@ -229,7 +229,7 @@ namespace AMDiS {
std::vector<Operator*>::iterator op;
for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
- (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
+ (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
for (int face = 0; face < dim + 1; face++) {
@@ -242,43 +242,41 @@ namespace AMDiS {
(*robinOperators)[face]->getAssembler(omp_get_thread_num())->
getZeroOrderAssembler()->getQuadrature();
- numPoints = quadrature->getNumPoints();
+ int nPoints = quadrature->getNumPoints();
if (elInfo->getBoundary(face) == boundaryType) {
- (*neumannOperators)[face]->getAssembler(omp_get_thread_num())->getZeroOrderAssembler()->
- initElement(elInfo);
+ (*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
+ getZeroOrderAssembler()->initElement(elInfo);
- const double *uhAtQp = dv->getVecAtQPs(elInfo,
- quadrature,
- NULL,
- NULL);
+ const double *uhAtQp = dv->getVecAtQPs(elInfo, quadrature, NULL, NULL);
- std::vector<double> f(numPoints, 0.0);
+ std::vector<double> f(nPoints, 0.0);
if (robinOperators)
- (*robinOperators)[face]->evalZeroOrder(numPoints,
+ (*robinOperators)[face]->evalZeroOrder(nPoints,
uhAtQp,
NULL,
NULL,
&f[0],
1.0);
- std::vector<WorldVector<double> > grdUh(numPoints);
- std::vector<WorldVector<double> > A_grdUh(numPoints);
+ std::vector<WorldVector<double> > grdUh(nPoints);
+ std::vector<WorldVector<double> > A_grdUh(nPoints);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
A_grdUh[iq].set(0.0);
lambda = quadrature->getLambda(iq);
basFcts->evalGrdUh(lambda, Lambda, uhEl, &grdUh[iq]);
}
-
+
for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
- (*op)->weakEvalSecondOrder(grdUh, A_grdUh);
+ (*op)->weakEvalSecondOrder(grdUh, A_grdUh);
if (neumannOperators)
- (*neumannOperators)[face]->getC(elInfo, numPoints, f);
+ (*neumannOperators)[face]->getC(elInfo, nPoints, f);
- for (val = iq = 0; iq < numPoints; iq++) {
+ val = 0.0;
+ for (int iq = 0; iq < nPoints; iq++) {
n_A_grdUh = (normal*A_grdUh[iq]) - f[iq];
val += quadrature->getWeight(iq)*sqr(n_A_grdUh);
}
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index 2e41ca904073685a397688e2e27561497854f83f..951fe1c839ab8eab35de4fe33599bcdede8eeea3 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -166,13 +166,13 @@ namespace AMDiS {
Quadrature *localQuad = quad ? quad : quadrature;
if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
- return valuesAtQPs[vec]->values.getValArray();
+ return &(valuesAtQPs[vec]->values[0]);
if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
- double *values = valuesAtQPs[vec]->values.getValArray();
+ double *values = &(valuesAtQPs[vec]->values[0]);
bool sameFESpaces =
(vec->getFESpace() == owner->rowFESpace) ||
(vec->getFESpace() == owner->colFESpace);
@@ -225,7 +225,7 @@ namespace AMDiS {
valuesAtQPs[vec] = new ValuesAtQPs;
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
- double *values = valuesAtQPs[vec]->values.getValArray();
+ double *values = &(valuesAtQPs[vec]->values[0]);
valuesAtQPs[vec]->valid = true;
vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
@@ -244,7 +244,7 @@ namespace AMDiS {
TEST_EXIT_DBG(vec)("no dof vector!\n");
if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid)
- return gradientsAtQPs[vec]->values.getValArray();
+ return &(gradientsAtQPs[vec]->values[0]);
Quadrature *localQuad = quad ? quad : quadrature;
@@ -253,7 +253,7 @@ namespace AMDiS {
gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());
- WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();
+ WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]);
const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
const BasisFunction *phi = owner->colFESpace->getBasisFcts();
@@ -299,7 +299,7 @@ namespace AMDiS {
TEST_EXIT_DBG(vec)("no dof vector!\n");
if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid)
- return gradientsAtQPs[vec]->values.getValArray();
+ return &(gradientsAtQPs[vec]->values[0]);
Quadrature *localQuad = quad ? quad : quadrature;
@@ -307,7 +307,7 @@ namespace AMDiS {
gradientsAtQPs[vec] = new GradientsAtQPs;
gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());
- WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();
+ WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]);
gradientsAtQPs[vec]->valid = true;
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index 8e66d1ab1b433f8f84304c58ddf5e53ff919a8af..55ddf5abc3a7d5a41a6363ca52813ab6fd3f9187 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -162,14 +162,24 @@ namespace AMDiS {
/// Used for \ref getVectorAtQPs().
class ValuesAtQPs {
public:
- Vector<double> values;
+ ValuesAtQPs()
+ : values(0),
+ valid(false)
+ {}
+
+ std::vector<double> values;
bool valid;
};
/// Used for \ref getGradientsAtQPs().
class GradientsAtQPs {
public:
- Vector<WorldVector<double> > values;
+ GradientsAtQPs()
+ : values(0),
+ valid(false)
+ {}
+
+ std::vector<WorldVector<double> > values;
bool valid;
};