Commit 5e18632e authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

No commit message

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