diff --git a/AMDiS/src/BasisFunction.cc b/AMDiS/src/BasisFunction.cc
index 2578159099ba572e7b4f309e9a46faac834f380f..e66b3c5aa1cd20d0c8d244a31a1c7f37aef92112 100644
--- a/AMDiS/src/BasisFunction.cc
+++ b/AMDiS/src/BasisFunction.cc
@@ -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 DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, WorldVector<double>* grd_uh) const
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index a9db74ffe899314524477c59dae49d5c373679ab..0cd8a50ab1cafe3e067cdbed407f4ce962bd07f8 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -325,6 +325,15 @@ namespace AMDiS {
*/
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
* Evaluates the gradient at barycentric coordinates lambda. Lambda is the
* Jacobian of the barycentric coordinates. uh is the local coefficient
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 90af161cbc4371ce79dd91f873ca169dc2c29b64..7428c46bce7369bc4f5a92bdf5d50f5efb1e7dd5 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -422,9 +422,8 @@ namespace AMDiS {
template<>
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 BasisFunction *basisFcts = feSpace->getBasisFcts();
@@ -438,22 +437,22 @@ namespace AMDiS {
DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, nBasisFcts);
double *sourceLocalCoeffs = GET_MEMORY(double, nSourceBasisFcts);
- if(feSpace->getMesh() == sourceFeSpace->getMesh()) {
+ if (feSpace->getMesh() == sourceFeSpace->getMesh()) {
DimVec<double> *coords = NULL;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS);
- while(elInfo) {
+ while (elInfo) {
Element *el = elInfo->getElement();
basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
source->getLocalVector(el, sourceLocalCoeffs);
- for(i = 0; i < nBasisFcts; i++) {
- if(vec[localIndices[i]] == 0.0) {
+ for (int i = 0; i < nBasisFcts; i++) {
+ if (vec[localIndices[i]] == 0.0) {
coords = basisFcts->getCoords(i);
vec[localIndices[i]] = sourceBasisFcts->evalUh(*coords, sourceLocalCoeffs) * factor;
}
@@ -475,28 +474,28 @@ namespace AMDiS {
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS,
&elInfo1, &elInfo2,
&elInfoSmall, &elInfoLarge);
- while(nextTraverse) {
+ while (nextTraverse) {
basisFcts->getLocalIndices(elInfo1->getElement(),
feSpace->getAdmin(),
localIndices);
source->getLocalVector(elInfo2->getElement(),
sourceLocalCoeffs);
- for(i = 0; i < nBasisFcts; i++) {
- if(vec[localIndices[i]] == 0.0) {
+ for (int i = 0; i < nBasisFcts; i++) {
+ if (vec[localIndices[i]] == 0.0) {
coords1 = basisFcts->getCoords(i);
elInfo1->coordToWorld(*coords1, &worldVec);
elInfo2->worldToCoord(worldVec, &coords2);
bool isPositive = true;
- for(j = 0; j < coords2.size(); j++) {
+ for (int j = 0; j < coords2.size(); j++) {
if (coords2[j] < 0) {
isPositive = false;
break;
}
}
- if(isPositive) {
+ if (isPositive) {
vec[localIndices[i]] = sourceBasisFcts->evalUh(coords2, sourceLocalCoeffs);
}
}
@@ -511,6 +510,63 @@ namespace AMDiS {
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<>
WorldVector<DOFVector<double>*> *DOFVector<double>::getGradient(WorldVector<DOFVector<double>*> *grad) const
{
@@ -531,15 +587,15 @@ namespace AMDiS {
// define result vector
static WorldVector<DOFVector<double>*> *result = NULL;
- if(grad) {
+ if (grad) {
result = grad;
} else {
- if(!result) {
+ if (!result) {
result = NEW WorldVector<DOFVector<double>*>;
result->set(NULL);
}
- for(i = 0; i < dow; i++) {
- if((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
+ for (i = 0; i < dow; i++) {
+ if ((*result)[i] && (*result)[i]->getFESpace() != feSpace) {
DELETE (*result)[i];
(*result)[i] = NEW DOFVector<double>(feSpace, "gradient");
}
@@ -554,11 +610,11 @@ namespace AMDiS {
int numNodes = 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);
int numPositionNodes = mesh->getGeo(geoIndex);
int numPreDOFs = admin->getNumberOfPreDOFs(i);
- for(j = 0; j < numPositionNodes; j++) {
+ for (j = 0; j < numPositionNodes; j++) {
int dofs = basFcts->getNumberOfDOFs(geoIndex);
numNodeDOFs.push_back(dofs);
numDOFs += dofs;
@@ -573,7 +629,7 @@ namespace AMDiS {
TEST_EXIT(numDOFs == nBasFcts)
("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));
}
@@ -588,7 +644,7 @@ namespace AMDiS {
WorldVector<double> grd;
- while(elInfo) {
+ while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
@@ -597,10 +653,10 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
int localDOFNr = 0;
- for(i = 0; i < numNodes; i++) { // for all nodes
- for(j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
+ for (i = 0; i < numNodes; i++) { // for all nodes
+ for (j = 0; j < numNodeDOFs[i]; j++) { // for all dofs at this node
DegreeOfFreedom dofIndex = dof[i][numNodePreDOFs[localDOFNr] + j];
- if(!visited[dofIndex]) {
+ if (!visited[dofIndex]) {
grd = basFcts->evalGrdUh(*(bary[localDOFNr]),
grdLambda,
@@ -623,85 +679,10 @@ namespace AMDiS {
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>() {};
+
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
::std::vector<DegreeOfFreedom>::iterator it;
::std::vector<DegreeOfFreedom>::iterator end = vec.end();
@@ -711,6 +692,7 @@ namespace AMDiS {
}
};
+
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *res)
{
@@ -741,26 +723,5 @@ namespace AMDiS {
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
-
}
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index ea2b9376adace17e9dfd2a4dcdb04415fbc929be..28b3927ad84c216740a09a276866a4832f80df48 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -350,7 +350,9 @@ namespace AMDiS {
/** \brief
* Returns \ref coarsenOperation
*/
- inline CoarsenOperation getCoarsenOperation() { return coarsenOperation; };
+ inline CoarsenOperation getCoarsenOperation() {
+ return coarsenOperation;
+ };
/** \brief
* Restriction after coarsening. Implemented for DOFVector<double>
@@ -360,7 +362,9 @@ namespace AMDiS {
/** \brief
* Returns \ref refineInter
*/
- inline bool refineInterpol() { return refineInter; };
+ inline bool refineInterpol() {
+ return refineInter;
+ };
/** \brief
* Interpolation after refinement.
@@ -371,7 +375,9 @@ namespace AMDiS {
* Returns \ref vec
*/
// ::std::vector<T>& getVector() const {return (::std::vector<T>&) vec;} ;
- ::std::vector<T>& getVector() { return vec;} ;
+ ::std::vector<T>& getVector() {
+ return vec;
+ };
// /** \brief
// * Returns \ref feSpace
@@ -593,83 +599,12 @@ namespace AMDiS {
* Used by interpol while mesh traversal
*/
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
* Prints \ref vec to stdout
*/
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
* Computes the coefficients of the interpolant of the function fct and
* stores these in the DOFVector
@@ -678,18 +613,9 @@ namespace AMDiS {
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 =====
- void serialize(::std::ostream &out)
- {
+ void serialize(::std::ostream &out) {
unsigned int size = vec.size();
out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
@@ -703,49 +629,12 @@ namespace AMDiS {
};
DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
- // {
- // FUNCNAME("DOFVector<T>::getGradient()");
- // ERROR_EXIT("only specialized for DOFVector<double>\n");
- // return NULL;
- // };
WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;
DOFVector<WorldVector<T> >*
- getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
- // {
- // FUNCNAME("DOFVector<T>::getRecoveryGradient()");
- // ERROR_EXIT("only specialized for DOFVector<double>\n");
- // return NULL;
- // };
-
- //const T *getLocalVector(const Element *el, T* localVec);
-
- // const T *getVecAtQPs(const ElInfo *elInfo,
- // const Quadrature *quad,
- // const FastQuadrature *quadFast,
- // T *vecAtQPs) const;
-
- // const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo,
- // const Quadrature *quad,
- // const FastQuadrature *quadFast,
- // WorldVector<T> *grdAtQPs) const;
- // {
- // FUNCNAME("DOFVector<T>::getGrdAtQPs()");
- // ERROR_EXIT("only specialized for DOFVector<double>\n");
- // return NULL;
- // };
-
- // const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo,
- // const Quadrature *quad,
- // const FastQuadrature *quadFast,
- // WorldMatrix<T> *d2AtQPs) const;
- // {
- // FUNCNAME("DOFVector<T>::getD2AtQPs()");
- // ERROR_EXIT("only specialized for DOFVector<double>\n");
- // return NULL;
- // };
+ getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
protected:
@@ -1015,18 +904,8 @@ namespace AMDiS {
return vec->getUsedSize();
};
- //template<typename T>
- //void interpolGrd(DOFVector<T> *v,
- // WorldVector<DOFVector<T>*> *grad,
- // double factor, bool add = false);
-
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *result);
-#if 0
- void transform2(DOFVector<WorldVector<double> > *vec,
- ::std::vector<DOFVector<double>*> *result);
-#endif
-
}
#include "DOFVector.hh"
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 77e9a26a863a243575c1b1865fe17524be694c5a..4cc8e2411398c5da092ac74ebf9e4767e146f10b 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -647,37 +647,33 @@ namespace AMDiS {
template<typename T>
void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
{
- FUNCNAME("interpol");
+ FUNCNAME("DOFVector<T>::interpol()");
TEST_EXIT(interFct = fct)("no function to interpolate\n");
- if (!this->getFESpace())
- {
- MSG("no dof admin in vec %s, skipping interpolation\n",
- this->getName().c_str());
- return;
- }
+ if (!this->getFESpace()) {
+ MSG("no dof admin in vec %s, skipping interpolation\n",
+ this->getName().c_str());
+ return;
+ }
- if (!(this->getFESpace()->getAdmin()))
- {
- MSG("no dof admin in feSpace %s, skipping interpolation\n",
- this->getFESpace()->getName().c_str());
- return;
- }
+ if (!(this->getFESpace()->getAdmin())) {
+ MSG("no dof admin in feSpace %s, skipping interpolation\n",
+ this->getFESpace()->getName().c_str());
+ return;
+ }
- if (!(this->getFESpace()->getBasisFcts()))
- {
- MSG("no basis functions in admin of vec %s, skipping interpolation\n",
- this->getName().c_str());
- return;
- }
+ if (!(this->getFESpace()->getBasisFcts())) {
+ MSG("no basis functions in admin of vec %s, skipping interpolation\n",
+ this->getName().c_str());
+ return;
+ }
- if (!(fct))
- {
- MSG("function that should be interpolated only pointer to NULL, ");
- Msg::print("skipping interpolation\n");
- return;
- }
+ if (!(fct)) {
+ MSG("function that should be interpolated only pointer to NULL, ");
+ Msg::print("skipping interpolation\n");
+ return;
+ }
traverseVector = this;
this->getFESpace()->getMesh()->traverse(-1,