Commit 78f007ce authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* New interpol-function for DOFVector<WorldVector<double>>

parent 196751d4
...@@ -46,6 +46,27 @@ namespace AMDiS { ...@@ -46,6 +46,27 @@ namespace AMDiS {
} }
const WorldVector<double>& BasisFunction::evalUh(const DimVec<double>& lambda,
const WorldVector<double> *uh_loc,
WorldVector<double>* values) const
{
static WorldVector<double> Values(DEFAULT_VALUE, 0.);
WorldVector<double> *val = (NULL != values) ? values : &Values;
int dow = Global::getGeo(WORLD);
for (int n = 0; n < dow; n++)
(*val)[n] = 0;
for (int i = 0; i < nBasFcts; i++) {
double phil = (*(*phi)[i])(lambda);
for (int n = 0; n < dow; n++)
(*val)[n] += uh_loc[i][n] * phil;
}
return((*val));
}
const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda, const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda, const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, WorldVector<double>* grd_uh) const const double *uh_loc, WorldVector<double>* grd_uh) const
......
...@@ -325,6 +325,15 @@ namespace AMDiS { ...@@ -325,6 +325,15 @@ namespace AMDiS {
*/ */
double evalUh(const DimVec<double>& lambda, const double* uh) const; double evalUh(const DimVec<double>& lambda, const double* uh) const;
/** \brief
* Evaluates elements value at barycentric coordinates lambda with local
* coefficient vector uh. If val is not NULL the result will be stored
* there, otherwise a pointer to a static local variable is returned which
* will be overwritten after the next call.
*/
const WorldVector<double>& evalUh(const DimVec<double>& lambda,
const WorldVector<double>* uh, WorldVector<double>* val) const;
/** \brief /** \brief
* Evaluates the gradient at barycentric coordinates lambda. Lambda is the * Evaluates the gradient at barycentric coordinates lambda. Lambda is the
* Jacobian of the barycentric coordinates. uh is the local coefficient * Jacobian of the barycentric coordinates. uh is the local coefficient
......
...@@ -422,9 +422,8 @@ namespace AMDiS { ...@@ -422,9 +422,8 @@ namespace AMDiS {
template<> template<>
void DOFVector<double>::interpol(DOFVector<double> *source, double factor) void DOFVector<double>::interpol(DOFVector<double> *source, double factor)
{ {
FUNCNAME("DOFVector<double>::interpol"); FUNCNAME("DOFVector<double>::interpol()");
int i, j;
const FiniteElemSpace *sourceFeSpace = source->getFESpace(); const FiniteElemSpace *sourceFeSpace = source->getFESpace();
const BasisFunction *basisFcts = feSpace->getBasisFcts(); const BasisFunction *basisFcts = feSpace->getBasisFcts();
...@@ -438,22 +437,22 @@ namespace AMDiS { ...@@ -438,22 +437,22 @@ namespace AMDiS {
DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, nBasisFcts); DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, nBasisFcts);
double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts); double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts);
if(feSpace->getMesh() == sourceFeSpace->getMesh()) { if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
DimVec<double> *coords = NULL; DimVec<double> *coords = NULL;
TraverseStack stack; TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
Mesh::CALL_LEAF_EL | Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS); Mesh::FILL_COORDS);
while(elInfo) { while (elInfo) {
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
source->getLocalVector(el, sourceLocalCoeffs); source->getLocalVector(el, sourceLocalCoeffs);
for(i = 0; i < nBasisFcts; i++) { for (int i = 0; i < nBasisFcts; i++) {
if(vec[localIndices[i]] == 0.0) { if (vec[localIndices[i]] == 0.0) {
coords = basisFcts->getCoords(i); coords = basisFcts->getCoords(i);
vec[localIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor; vec[localIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
} }
...@@ -475,28 +474,28 @@ namespace AMDiS { ...@@ -475,28 +474,28 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS,
&elInfo1, &elInfo2, &elInfo1, &elInfo2,
&elInfoSmall, &elInfoLarge); &elInfoSmall, &elInfoLarge);
while(nextTraverse) { while (nextTraverse) {
basisFcts->getLocalIndices(elInfo1->getElement(), basisFcts->getLocalIndices(elInfo1->getElement(),
feSpace->getAdmin(), feSpace->getAdmin(),
localIndices); localIndices);
source->getLocalVector(elInfo2->getElement(), source->getLocalVector(elInfo2->getElement(),
sourceLocalCoeffs); sourceLocalCoeffs);
for(i = 0; i < nBasisFcts; i++) { for (int i = 0; i < nBasisFcts; i++) {
if(vec[localIndices[i]] == 0.0) { if (vec[localIndices[i]] == 0.0) {
coords1 = basisFcts->getCoords(i); coords1 = basisFcts->getCoords(i);
elInfo1->coordToWorld(*coords1, &worldVec); elInfo1->coordToWorld(*coords1, &worldVec);
elInfo2->worldToCoord(worldVec, &coords2); elInfo2->worldToCoord(worldVec, &coords2);
bool isPositive = true; bool isPositive = true;
for(j = 0; j < coords2.size(); j++) { for (int j = 0; j < coords2.size(); j++) {
if (coords2[j] < 0) { if (coords2[j] < 0) {
isPositive = false; isPositive = false;
break; break;
} }
} }
if(isPositive) { if (isPositive) {
vec[localIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs); vec[localIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);
} }
} }
...@@ -511,6 +510,63 @@ namespace AMDiS { ...@@ -511,6 +510,63 @@ namespace AMDiS {
FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts); FREE_MEMORY(sourceLocalCoeffs, double, nSourceBasisFcts);
} }
template<>
void DOFVector<WorldVector<double> >::interpol(DOFVector<WorldVector<double> > *v, double factor)
{
WorldVector<double> nul(DEFAULT_VALUE,0.0);
this->set(nul);
DimVec<double> *coords = NULL;
const FiniteElemSpace *vFESpace = v->getFESpace();
if (feSpace == vFESpace) {
WARNING("same FE-spaces\n");
}
const BasisFunction *basFcts = feSpace->getBasisFcts();
const BasisFunction *vBasFcts = vFESpace->getBasisFcts();
int numBasFcts = basFcts->getNumber();
int vNumBasFcts = vBasFcts->getNumber();
if (feSpace->getMesh() == vFESpace->getMesh()) {
DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts);
WorldVector<double> *vLocalCoeffs = NEW WorldVector<double>[vNumBasFcts];
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS);
while (elInfo) {
Element *el = elInfo->getElement();
basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
v->getLocalVector(el, vLocalCoeffs);
for (int i = 0; i < numBasFcts; i++) {
if (vec[localIndices[i]] == nul) {
coords = basFcts->getCoords(i);
// for(j = 0; j < vNumBasFcts; j++) {
vec[localIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs,NULL) * factor;
// }
}
}
elInfo = stack.traverseNext(elInfo);
}
FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
DELETE [] vLocalCoeffs;
} else {
ERROR_EXIT("not yet for dual traverse\n");
}
}
template<> template<>
WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
{ {
...@@ -531,15 +587,15 @@ namespace AMDiS { ...@@ -531,15 +587,15 @@ namespace AMDiS {
// define result vector // define result vector
static WorldVector<DOFVector<double>*> *result = NULL; static WorldVector<DOFVector<double>*> *result = NULL;
if(grad) { if (grad) {
result = grad; result = grad;
} else { } else {
if(!result) { if (!result) {
result = NEW WorldVector<DOFVector<double>*>; result = NEW WorldVector<DOFVector<double>*>;
result->set(NULL); result->set(NULL);
} }
for(i = 0; i < dow; i++) { for (i = 0; i < dow; i++) {
if((*result)[i] && (*result)[i]->getFESpace() != feSpace) { if ((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
DELETE (*result)[i]; DELETE (*result)[i];
(*result)[i] = NEW DOFVector<double>(feSpace, "gradient"); (*result)[i] = NEW DOFVector<double>(feSpace, "gradient");
} }
...@@ -554,11 +610,11 @@ namespace AMDiS { ...@@ -554,11 +610,11 @@ namespace AMDiS {
int numNodes = 0; int numNodes = 0;
int numDOFs = 0; int numDOFs = 0;
for(i = 0; i < dim + 1; i++) { for (i = 0; i < dim + 1; i++) {
GeoIndex geoIndex = INDEX_OF_DIM(i, dim); GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
int numPositionNodes = mesh->getGeo(geoIndex); int numPositionNodes = mesh->getGeo(geoIndex);
int numPreDOFs = admin->getNumberOfPreDOFs(i); int numPreDOFs = admin->getNumberOfPreDOFs(i);
for(j = 0; j < numPositionNodes; j++) { for (j = 0; j < numPositionNodes; j++) {
int dofs = basFcts->getNumberOfDOFs(geoIndex); int dofs = basFcts->getNumberOfDOFs(geoIndex);
numNodeDOFs.push_back(dofs); numNodeDOFs.push_back(dofs);
numDOFs += dofs; numDOFs += dofs;
...@@ -573,7 +629,7 @@ namespace AMDiS { ...@@ -573,7 +629,7 @@ namespace AMDiS {
TEST_EXIT(numDOFs == nBasFcts) TEST_EXIT(numDOFs == nBasFcts)
("number of dofs != number of basis functions\n"); ("number of dofs != number of basis functions\n");
for(i = 0; i < numDOFs; i++) { for (i = 0; i < numDOFs; i++) {
bary.push_back(basFcts->getCoords(i)); bary.push_back(basFcts->getCoords(i));
} }
...@@ -588,7 +644,7 @@ namespace AMDiS { ...@@ -588,7 +644,7 @@ namespace AMDiS {
WorldVector<double> grd; WorldVector<double> grd;
while(elInfo) { while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
...@@ -597,10 +653,10 @@ namespace AMDiS { ...@@ -597,10 +653,10 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
int localDOFNr = 0; int localDOFNr = 0;
for(i = 0; i < numNodes; i++) { // for all nodes for (i = 0; i < numNodes; i++) { // for all nodes
for(j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node for (j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j]; DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
if(!visited[dofIndex]) { if (!visited[dofIndex]) {
grd = basFcts->evalGrdUh(*(bary[localDOFNr]), grd = basFcts->evalGrdUh(*(bary[localDOFNr]),
grdLambda, grdLambda,
...@@ -623,85 +679,10 @@ namespace AMDiS { ...@@ -623,85 +679,10 @@ namespace AMDiS {
return result; return result;
} }
// template<>
// void interpolGrd(DOFVector<double> *v,
// WorldVector<DOFVector<double>*> *grad,
// double factor, bool add)
// {
// int i, j, k;
// int dow = Global::getGeo(WORLD);
// TEST_EXIT(grad)("no gradient\n");
// for(i = 0; i < dow; i++) {
// TEST_EXIT((*grad)[i])("missing (*grad)[%d]\n", i);
// }
// if(!add) {
// for(i = 0; i < dow; i++) {
// (*grad)[i]->set(0.0);
// }
// }
// DimVec<double> *coords = NULL;
// const FiniteElemSpace *feSpace = (*grad)[0]->getFESpace();
// const FiniteElemSpace *vFESpace = v->getFESpace();
// if(feSpace == vFESpace) {
// WARNING("same FE-spaces\n");
// }
// const BasisFunction *basFcts = feSpace->getBasisFcts();
// const BasisFunction *vBasFcts = vFESpace->getBasisFcts();
// int numBasFcts = basFcts->getNumber();
// int vNumBasFcts = vBasFcts->getNumber();
// if(feSpace->getMesh() == vFESpace->getMesh()) {
// DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts);
// double *vLocalCoeffs = GET_MEMORY(double, vNumBasFcts);
// Mesh *mesh = feSpace->getMesh();
// TraverseStack stack;
// ElInfo *elInfo = stack.traverseFirst(mesh, -1,
// Mesh::CALL_LEAF_EL |
// Mesh::FILL_COORDS |
// Mesh::FILL_GRD_LAMBDA);
// while(elInfo) {
// const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
// Element *el = elInfo->getElement();
// basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
// v->getLocalVector(el, vLocalCoeffs);
// for(i = 0; i < numBasFcts; i++) {
// coords = basFcts->getCoords(i);
// for(j = 0; j < vNumBasFcts; j++) {
// WorldVector<double> g;
// vBasFcts->evalGrdUh(*coords, grdLambda, vLocalCoeffs, &g);
// for(k = 0; k < dow; k++) {
// (*((*grad)[k]))[localIndices[i]] += g[k] * factor;
// }
// }
// }
// elInfo = stack.traverseNext(elInfo);
// }
// FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
// FREE_MEMORY(vLocalCoeffs, double, vNumBasFcts);
// } else {
// ERROR_EXIT("not yet for dual traverse\n");
// }
// }
DOFVectorDOF::DOFVectorDOF() : DOFVector<DegreeOfFreedom>() {}; DOFVectorDOF::DOFVectorDOF() : DOFVector<DegreeOfFreedom>() {};
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) { void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
::std::vector<DegreeOfFreedom>::iterator it; ::std::vector<DegreeOfFreedom>::iterator it;
::std::vector<DegreeOfFreedom>::iterator end = vec.end(); ::std::vector<DegreeOfFreedom>::iterator end = vec.end();
...@@ -711,6 +692,7 @@ namespace AMDiS { ...@@ -711,6 +692,7 @@ namespace AMDiS {
} }
}; };
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec, WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *res) WorldVector<DOFVector<double>*> *res)
{ {
...@@ -741,26 +723,5 @@ namespace AMDiS { ...@@ -741,26 +723,5 @@ namespace AMDiS {
return r; return r;
} }
#if 0
void transform2(DOFVector<WorldVector<double> > *vec,
::std::vector<DOFVector<double>*> *res)
{
FUNCNAME("transform2()");
TEST_EXIT(vec && res)("no vector\n");
int i, j, dow = Global::getGeo(WORLD);
int vecSize = vec->getSize();
for(i = 0; i < vecSize; i++) {
for(j = 0; j < dow; j++) {
(*((*res)[j]))[i] = (*vec)[i][j];
}
}
return r;
}
#endif
} }
...@@ -350,7 +350,9 @@ namespace AMDiS { ...@@ -350,7 +350,9 @@ namespace AMDiS {
/** \brief /** \brief
* Returns \ref coarsenOperation * Returns \ref coarsenOperation
*/ */
inline CoarsenOperation getCoarsenOperation() { return coarsenOperation; }; inline CoarsenOperation getCoarsenOperation() {
return coarsenOperation;
};
/** \brief /** \brief
* Restriction after coarsening. Implemented for DOFVector<double> * Restriction after coarsening. Implemented for DOFVector<double>
...@@ -360,7 +362,9 @@ namespace AMDiS { ...@@ -360,7 +362,9 @@ namespace AMDiS {
/** \brief /** \brief
* Returns \ref refineInter * Returns \ref refineInter
*/ */
inline bool refineInterpol() { return refineInter; }; inline bool refineInterpol() {
return refineInter;
};
/** \brief /** \brief
* Interpolation after refinement. * Interpolation after refinement.
...@@ -371,7 +375,9 @@ namespace AMDiS { ...@@ -371,7 +375,9 @@ namespace AMDiS {
* Returns \ref vec * Returns \ref vec
*/ */
// ::std::vector<T>& getVector() const {return (::std::vector<T>&) vec;} ; // ::std::vector<T>& getVector() const {return (::std::vector<T>&) vec;} ;
::std::vector<T>& getVector() { return vec;} ; ::std::vector<T>& getVector() {
return vec;
};
// /** \brief // /** \brief
// * Returns \ref feSpace // * Returns \ref feSpace
...@@ -593,83 +599,12 @@ namespace AMDiS { ...@@ -593,83 +599,12 @@ namespace AMDiS {
* Used by interpol while mesh traversal * Used by interpol while mesh traversal
*/ */
static int interpolFct(ElInfo* elinfo); static int interpolFct(ElInfo* elinfo);
/** \brief
* Generalized matrix-vector multiplication
*/
//virtual void gemv(MatrixTranspose transpose, T alpha,
// const DOFMatrix &a, const DOFVector<T>& x,
// T beta);
/** \brief
* Matrix-vector multiplication
*/
//virtual void mv(MatrixTranspose transpose,
// const DOFMatrix&a, const DOFVector<T>&x);
/** \brief /** \brief
* Prints \ref vec to stdout * Prints \ref vec to stdout
*/ */
void print() const; void print() const;
// ElementVector *assemble(T factor, ElInfo *elInfo,
// const BoundaryType *bound,
// Operator *op = NULL);
// /** \brief
// * Adds element contribution to the vector
// */
// void addElementVector(T sign,
// const ElementVector &elVec,
// const BoundaryType *bound,
// bool add = true);
// inline void addOperator(Operator* op,
// double *factor = NULL,
// double *estFactor = NULL)
// {
// operators.push_back(op);
// operatorFactor.push_back(factor);
// operatorEstFactor.push_back(estFactor);
// };
// inline ::std::vector<double*>::iterator getOperatorFactorBegin() {
// return operatorFactor.begin();
// };
// inline ::std::vector<double*>::iterator getOperatorFactorEnd() {
// return operatorFactor.end();
// };
// inline ::std::vector<double*>::iterator getOperatorEstFactorBegin() {
// return operatorEstFactor.begin();
// };
// inline ::std::vector<double*>::iterator getOperatorEstFactorEnd() {
// return operatorEstFactor.end();
// };
// inline ::std::vector<Operator*>::iterator getOperatorsBegin() {
// return operators.begin();
// };
// inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
// return operators.end();
// };
// Flag getAssembleFlag();
// /** \brief
// * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] *
// * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
// * \f$ x(\lambda) \f$ are the world coordinates of lambda and
// * \f$ m \f$ is the number of basis functions
// */
// T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);
/** \brief /** \brief
* Computes the coefficients of the interpolant of the function fct and * Computes the coefficients of the interpolant of the function fct and
* stores these in the DOFVector * stores these in the DOFVector
...@@ -678,18 +613,9 @@ namespace AMDiS { ...@@ -678,18 +613,9 @@ namespace AMDiS {
void interpol(DOFVector<T> *v, double factor); void interpol(DOFVector<T> *v, double factor);
// inline ::std::vector<Operator*>& getOperators() { return operators; };
// inline ::std::vector<double*>& getOperatorFactor() { return operatorFactor; };
// inline ::std::vector<double*>& getOperatorEstFactor() {
// return operatorEstFactor;
// };
// ===== Serializable implementation ===== // ===== Serializable implementation =====
void serialize(::std::ostream &out) void serialize(::std::ostream &out) {
{
unsigned int size = vec.size(); unsigned int size = vec.size();