Skip to content
GitLab
Commit 9bbe8339 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

evalAtPoint

parent 7856d906
Hide whitespace changes
Inline Side-by-side
AMDiS/src/BoundaryManager.cc
View file @ 9bbe8339
......@@ -71,16 +71,20 @@ namespace AMDiS {
basisFcts->getBound(elInfo, localBound);
// get dof indices
basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
basisFcts->getLocalIndices(elInfo->getElement(),
feSpace->getAdmin(),
dofVec);
// apply non dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
for (BoundaryIndexMap::iterator it = localBCs.begin();
it != localBCs.end(); ++it)
if ((*it).second && !(*it).second->isDirichlet())
(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
localBound, nBasFcts);
// apply dirichlet boundary conditions
for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
for (BoundaryIndexMap::iterator it = localBCs.begin();
it != localBCs.end(); ++it)
if ((*it).second && (*it).second->isDirichlet())
(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
localBound, nBasFcts);
......@@ -101,7 +105,7 @@ namespace AMDiS {
// get boundaries of all DOFs
basisFcts->getBound(elInfo, localBound);
// get dof indices
// get DOF indices
basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
// apply non dirichlet boundary conditions
......
AMDiS/src/BoundaryManager.h
View file @ 9bbe8339
......@@ -137,12 +137,12 @@ namespace AMDiS {
/** \brief
* For every boundary id we store here all possible boundary object (although
* it's not clear if it is meaningful to have different boundary conditions on the
* same boundary id).
* it's not clear if it is meaningful to have different boundary conditions on
* the same boundary id).
*
* We have to use this global variable, because the mesh traverse interface does
* not provide more information about traversed boundaries at elements than the
* boundary id.
* We have to use this global variable, because the mesh traverse interface
* does not provide more information about traversed boundaries at elements
* than the boundary id.
*
* TODO: Change interface such that mesh traverse returns the boundary objects
* directly and we can remove this global variable. The biggest problem will be
......
AMDiS/src/DOFVector.cc
View file @ 9bbe8339
......@@ -68,9 +68,10 @@ namespace AMDiS {
template<>
const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
const double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
ElInfo *oldElInfo) const
{
FUNCNAME("DOFVector<double>::evalAtCoords()");
FUNCNAME("DOFVector<double>::evalAtPoint()");
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basFcts = feSpace->getBasisFcts();
......@@ -78,7 +79,7 @@ namespace AMDiS {
int dim = mesh->getDim();
int nBasFcts = basFcts->getNumber();
DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
std::vector<DegreeOfFreedom> localIndices(nBasFcts);
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
......@@ -86,7 +87,8 @@ namespace AMDiS {
bool inside = false;
if (oldElInfo && oldElInfo->getMacroElement()) {
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
inside = mesh->findElInfoAtPoint(p, elInfo, lambda,
oldElInfo->getMacroElement(), NULL, NULL);
delete oldElInfo;
} else
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
......@@ -95,7 +97,8 @@ namespace AMDiS {
oldElInfo = elInfo;
if (inside) {
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
localIndices);
ElementVector uh(lambda.size());
for (int i = 0; i < lambda.size(); i++)
uh[i] = operator[](localIndices[i]);
......@@ -103,20 +106,19 @@ namespace AMDiS {
} else
throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
delete [] localIndices;
if (oldElInfo == NULL)
delete elInfo;
if (values != NULL)
*values = value;
return value;
};
}
template<>
const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
const WorldVector<double>
DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
ElInfo *oldElInfo) const
{
FUNCNAME("DOFVector<double>::evalAtCoords()");
FUNCNAME("DOFVector<double>::evalAtPoint()");
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basFcts = feSpace->getBasisFcts();
......@@ -124,18 +126,17 @@ namespace AMDiS {
int dim = mesh->getDim();
int nBasFcts = basFcts->getNumber();
DegreeOfFreedom *localIndices = new DegreeOfFreedom[nBasFcts];
vector<DegreeOfFreedom> localIndices(nBasFcts);
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
WorldVector<double> *val = (NULL != values) ? values : &Values;
WorldVector<double> value(DEFAULT_VALUE, 0.0);
bool inside = false;
if (oldElInfo && oldElInfo->getMacroElement()) {
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
inside = mesh->findElInfoAtPoint(p, elInfo, lambda,
oldElInfo->getMacroElement(), NULL, NULL);
delete oldElInfo;
} else
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
......@@ -144,20 +145,20 @@ namespace AMDiS {
oldElInfo = elInfo;
if (inside) {
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
localIndices);
mtl::dense_vector<WorldVector<double> > uh(lambda.size());
for (int i = 0; i < lambda.size(); i++)
uh[i] = operator[](localIndices[i]);
*val = basFcts->evalUh(lambda, uh);
value = basFcts->evalUh(lambda, uh);
} else
throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
delete [] localIndices;
if (oldElInfo == NULL)
delete elInfo;
return ((*val));
};
return value;
}
template<>
......
AMDiS/src/DOFVector.h
View file @ 9bbe8339
......@@ -65,17 +65,13 @@ namespace AMDiS {
virtual ~DOFVectorBase();
/** \brief
* For the given element, this function returns an array of all DOFs of this
* DOFVector that are defined on this element.
*/
/// For the given element, this function returns an array of all DOFs of this
/// DOFVector that are defined on this element.
virtual void getLocalVector(const Element *el,
mtl::dense_vector<T>& localVec) const;
/** \brief
* Evaluates the DOF vector at a set of quadrature points defined on the
* given element.
*/
/// Evaluates the DOF vector at a set of quadrature points defined on the
/// given element.
void getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
......@@ -87,10 +83,8 @@ namespace AMDiS {
const FastQuadrature *quadFast,
mtl::dense_vector<T>& vecAtQPs) const;
/** \brief
* Evaluates the gradient of a DOF vector at a set of quadrature points defined on the
* given element.
*/
/// Evaluates the gradient of a DOF vector at a set of quadrature points
/// defined on the given element.
void getGrdAtQPs( const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
......@@ -102,10 +96,8 @@ namespace AMDiS {
const FastQuadrature *quadFast,
mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;
/** \brief
* Evaluates the comp'th component of the derivative of a DOF vector at a set of quadrature points defined on the
* given element.
*/
/// Evaluates the comp'th component of the derivative of a DOF vector at a
/// set of quadrature points defined on the given element.
void getDerivativeAtQPs( const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
......@@ -119,10 +111,8 @@ namespace AMDiS {
int comp,
mtl::dense_vector<T> &derivativeAtQPs) const;
/** \brief
* Evaluates the jacobian of a DOF vector at a set of quadrature points defined on the
* given element.
*/
/// Evaluates the jacobian of a DOF vector at a set of quadrature points
/// defined on the given element.
void getD2AtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
......@@ -134,10 +124,8 @@ namespace AMDiS {
return feSpace;
}
/** \brief
* Assembles the element vector for the given ellement and adds the
* element matrix to the current DOF vector.
*/
/// Assembles the element vector for the given ellement and adds the element
/// matrix to the current DOF vector.
void assemble(T factor, ElInfo *elInfo,
const BoundaryType *bound,
Operator *op = NULL);
......@@ -154,11 +142,9 @@ namespace AMDiS {
ElInfo *elInfo,
bool add = true);
/* \brief
* That function must be called after the matrix assembling has been finished.
* This makes it possible to start some cleanup or matrix data compressing
* procedures.
*/
/// That function must be called after the matrix assembling has been
/// finished. This makes it possible to start some cleanup or matrix data
/// compressing procedures.
void finishAssembling();
inline void addOperator(Operator* op,
......@@ -391,10 +377,8 @@ namespace AMDiS {
return vec.end();
}
/** \brief
* Used by DOFAdmin to compress this DOFVector. Implementation of
* DOFIndexedBase::compress()
*/
/// Used by DOFAdmin to compress this DOFVector. Implementation of
/// DOFIndexedBase::compress()
virtual void compressDOFIndexed(int first, int last,
std::vector<DegreeOfFreedom> &newDof);
......@@ -473,22 +457,18 @@ namespace AMDiS {
/// Calculates Integral of this DOFVector
double Int(Quadrature* q = NULL) const;
/** \brief
* Calculates Integral of this DOFVector over parts of the domain
* boundary, indicated by boundaryType. Implemented for DOFVector<double>
**/
/// Calculates Integral of this DOFVector over parts of the domain
/// boundary, indicated by boundaryType. Implemented for DOFVector<double>
T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const
{
FUNCNAME("DOFVector::IntOnBoundary())");
TEST_EXIT(false)("Please implement your integration\n");
return 0.0;
}
/** \brief
* Calculates Integral of this DOFVector times normal vector over parts
* of the domain boundary, indicated by boundaryType. Implemented for
* DOFVector<WorldVector<double> >
**/
/// Calculates Integral of this DOFVector times normal vector over parts
/// of the domain boundary, indicated by boundaryType. Implemented for
/// DOFVector<WorldVector<double> >
double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const
{
FUNCNAME("DOFVector::IntOnBoundaryNormal())");
......@@ -581,25 +561,21 @@ namespace AMDiS {
///
int calcMemoryUsage() const;
/** \brief
* Computes the coefficients of the interpolant of the function fct and
* stores these in the DOFVector
*/
/// Computes the coefficients of the interpolant of the function fct and
/// stores these in the DOFVector
void interpol(AbstractFunction<T, WorldVector<double> > *fct);
void interpol(DOFVector<T> *v, double factor = 1.0);
/** eval DOFVector at given point p. If oldElInfo != NULL the search for the element, where p is inside,
* starts from oldElInfo.
* implemented for: double, WorldVector< double >
*/
inline const T& evalAtPoint(
WorldVector<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const
/// Eval DOFVector at given point p. If oldElInfo != NULL the search for
/// the element, where p is inside, starts from oldElInfo. implemented for:
/// double, WorldVector<double>
inline const T evalAtPoint(WorldVector<double> &p,
ElInfo *oldElInfo = NULL) const
{
FUNCNAME("DOFVector::evalAtPoint())");
TEST_EXIT(false)("Please implement your evaluation\n");
return *value;
}
/// Writes the data of the DOFVector to an output stream.
......@@ -619,13 +595,11 @@ namespace AMDiS {
in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
}
// DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
DOFVector<typename GradientType<T>::type>*
getGradient(DOFVector<typename GradientType<T>::type> *grad) const;
WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;
// DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
DOFVector<typename GradientType<T>::type>*
getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;
......@@ -654,12 +628,13 @@ namespace AMDiS {
BoundaryType boundaryType, Quadrature* q) const;
template<>
const double& DOFVector<double>::evalAtPoint(
WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
const double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
ElInfo *oldElInfo) const;
template<>
const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
const WorldVector<double>
DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p,
ElInfo *oldElInfo) const;
template<>
void DOFVector<double>::refineInterpol(RCNeighbourList&, int);
......@@ -686,10 +661,8 @@ namespace AMDiS {
public DOFContainer
{
public:
/** \brief
* Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
* as DOFContainer at DOFAdmin
*/
/// Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
/// as DOFContainer at DOFAdmin
DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
: DOFVector<DegreeOfFreedom>(feSpace_, name_)
{
......@@ -702,10 +675,8 @@ namespace AMDiS {
feSpace->getAdmin()->removeDOFContainer(this);
}
/** \brief
* Implements DOFContainer::operator[]() by calling
* DOFVector<DegreeOfFreedom>::operator[]()
*/
/// Implements DOFContainer::operator[]() by calling
/// DOFVector<DegreeOfFreedom>::operator[]()
DegreeOfFreedom& operator[](DegreeOfFreedom i)
{
return DOFVector<DegreeOfFreedom>::operator[](i);
......@@ -869,20 +840,22 @@ namespace AMDiS {
std::vector<DOFVector<double>*> *res);
/** \brief
* Computes the integral of a function that includes two different DOFVectors. This
* function works also for the case that the DOFVectors are defined on two different
* meshes.
* Computes the integral of a function that includes two different DOFVectors.
* This function works also for the case that the DOFVectors are defined on
* two different meshes.
*/
double integrate(const DOFVector<double> &vec1,
const DOFVector<double> &vec2,
BinaryAbstractFunction<double, double, double> *fct);
/// Computes the integral of a function with two DOFVectors defined on the same mesh.
/// Computes the integral of a function with two DOFVectors defined on the
/// same mesh.
double intSingleMesh(const DOFVector<double> &vec1,
const DOFVector<double> &vec2,
BinaryAbstractFunction<double, double, double> *fct);
/// Computes the integral of a function with two DOFVectors defined on different meshes.
/// Computes the integral of a function with two DOFVectors defined on
/// different meshes.
double intDualMesh(const DOFVector<double> &vec1,
const DOFVector<double> &vec2,
BinaryAbstractFunction<double, double, double> *fct);
......
AMDiS/src/Mesh.cc
View file @ 9bbe8339
......@@ -589,7 +589,7 @@ namespace AMDiS {
{
const DOFAdmin *localAdmin = NULL;
for (int i = 0; i < static_cast<int>(admin.size()); i++) {
for (unsigned int i = 0; i < admin.size(); i++) {
if (admin[i]->getNumberOfDofs(VERTEX)) {
if (!localAdmin)
localAdmin = admin[i];
......
AMDiS/src/Mesh.h
View file @ 9bbe8339
......@@ -223,7 +223,8 @@ namespace AMDiS {
/// Returns a DOFAdmin which at least manages vertex DOFs
const DOFAdmin* getVertexAdmin() const;
/// Allocates a array of DOF pointers. The array holds one pointer for each node.
/// Allocates an array of DOF pointers. The array holds one pointer for
/// each node.
DegreeOfFreedom **createDofPtrs();
/// Returns \ref preserveCoarseDOFs of the mesh
......@@ -850,7 +851,8 @@ namespace AMDiS {
*/
Element* elementPrototype;
/// Prototype for leaf data. Used for creation of new leaf data while refinement.
/// Prototype for leaf data. Used for creation of new leaf data while
/// refinement.
ElementData* elementDataPrototype;
/// Used for enumeration of all mesh elements
......@@ -862,23 +864,20 @@ namespace AMDiS {
/// Map of managed periodic vertex associations.
map<BoundaryType, VertexVector*> periodicAssociations;
/** \brief
* If the mesh has been created by reading a macro file, here the information are
* are stored about the content of the file.
*/
/// If the mesh has been created by reading a macro file, here the information
/// are stored about the content of the file.
MacroInfo *macroFileInfo;
/** \brief
* This index is incremented every time the mesh is changed, e.g. by the refinement
* or the coarsening manager. It can be used by other object if the mesh has been
* changed by first copying this variable elsewhere and comparing its values.
*/
/// This index is incremented every time the mesh is changed, e.g. by the
/// refinement or the coarsening manager. It can be used by other object if
/// the mesh has been changed by first copying this variable elsewhere and
/// comparing its values.
long changeIndex;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
/// In parallel computations the mesh may be globally prerefined to achieve a fine
/// enought starting mesh for the given number of ranks. The value of the variable
/// will be defined in function \ref checkParallelMacroFile.
/// In parallel computations the mesh may be globally prerefined to achieve a
/// fine enought starting mesh for the given number of ranks. The value of the
/// variable will be defined in function \ref checkParallelMacroFile.
int nParallelPreRefinements;
#endif
......
AMDiS/src/PeriodicBC.cc
View file @ 9bbe8339
......@@ -122,8 +122,8 @@ namespace AMDiS {
Mesh *mesh = matrix->getRowFeSpace()->getMesh();
associated = mesh->getPeriodicAssociations()[boundaryType];
TEST_EXIT(associated)("No associations for periodic boundary condition %d!\n",
boundaryType);
TEST_EXIT(associated)
("No associations for periodic boundary condition %d!\n", boundaryType);
}
}
......@@ -136,65 +136,64 @@ namespace AMDiS {
{
FUNCNAME("PeriodicBC::fillBoundaryCondition()");
if (matrix == masterMatrix) {
if (matrix != masterMatrix)
return;
int dim = rowFeSpace->getMesh()->getDim();
if (dim > 1) {
DOFAdmin *admin = rowFeSpace->getAdmin();
FixVec<int, WORLD> elFace(dim, NO_INIT);
FixVec<int, WORLD> neighFace(dim, NO_INIT);
DimVec<int> vertexPermutation(dim, NO_INIT);
const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
int num = basFcts->getNumber();
Element *element = elInfo->getElement();
DimVec<DegreeOfFreedom> periodicDOFs(dim - 1, NO_INIT);
GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
std::vector<DegreeOfFreedom> neighIndices(num);
for (int side = 0; side < dim + 1; side++) {
if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
for (int vertex = 0; vertex < dim; vertex++) {
int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
periodicDOFs[vertex] = (*associated)[dofIndices[index]];
}
Element *neigh = elInfo->getNeighbour(side);
TEST_EXIT_DBG(neigh)("Wrong neighbour information at side %d!\n", side);
basFcts->getLocalIndices(neigh, admin, neighIndices);
int oppVertex = 0;
for (int i = 0; i < dim + 1; i++) {
// get vertex permutation
if (i == side) {
vertexPermutation[i] = 0;
} else {
DegreeOfFreedom periodicDOF =
periodicDOFs[element->getPositionOfVertex(side, i)];
int j = 0;
for (; j < dim + 1; j++)
if (neigh->getDof(j, 0) == periodicDOF)
break;
vertexPermutation[i] = j;
}
oppVertex += i - vertexPermutation[i];
}
vertexPermutation[side] = oppVertex;
// get DOF permutation
const DegreeOfFreedom *dofPermutation =
periodicDOFMapping->getDOFPermutation(vertexPermutation);
int dim = rowFeSpace->getMesh()->getDim();
if (dim == 1)
return;
DOFAdmin *admin = rowFeSpace->getAdmin();
FixVec<int, WORLD> elFace(dim, NO_INIT);
FixVec<int, WORLD> neighFace(dim, NO_INIT);
DimVec<int> vertexPermutation(dim, NO_INIT);
const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
int num = basFcts->getNumber();
Element *element = elInfo->getElement();
DimVec<DegreeOfFreedom> periodicDOFs(dim - 1, NO_INIT);
GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
std::vector<DegreeOfFreedom> neighIndices(num);
for (int side = 0; side <= dim; side++) {
if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
for (int vertex = 0; vertex < dim; vertex++) {
int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
periodicDOFs[vertex] = (*associated)[dofIndices[index]];
}
Element *neigh = elInfo->getNeighbour(side);
TEST_EXIT_DBG(neigh)("Wrong neighbour information at side %d!\n", side);
basFcts->getLocalIndices(neigh, admin, neighIndices);
// set associated dofs
for (int i = 0; i < num; i++)
if ((*(basFcts->getCoords(i)))[side] == 0)
(*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
int oppVertex = 0;
for (int i = 0; i < dim + 1; i++) {
// get vertex permutation
if (i == side) {
vertexPermutation[i] = 0;
} else {
DegreeOfFreedom periodicDOF =
periodicDOFs[element->getPositionOfVertex(side, i)];
int j = 0;
for (; j < dim + 1; j++)
if (neigh->getDof(j, 0) == periodicDOF)
break;
vertexPermutation[i] = j;
}
oppVertex += i - vertexPermutation[i];
}
vertexPermutation[side] = oppVertex;
// get DOF permutation
const DegreeOfFreedom *dofPermutation =
periodicDOFMapping->getDOFPermutation(vertexPermutation);
// set associated dofs
for (int i = 0; i < num; i++)
if ((*(basFcts->getCoords(i)))[side] == 0)
(*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
}
}
}
......
AMDiS/src/PeriodicBC.h
View file @ 9bbe8339
......@@ -31,6 +31,8 @@
namespace AMDiS {
using namespace std;
template<typename T>
class DimVecLess
{
......@@ -61,13 +63,18 @@ namespace AMDiS {
const DegreeOfFreedom *getDOFPermutation(const DimVec<int> &vertexPermutation);
protected: