Commit c1eca4e2 authored by Thomas Witkowski's avatar Thomas Witkowski

UAH, IST DAS GROSSSSSSS

parent 56c306c7
......@@ -125,6 +125,7 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
MeshDistributor::globalMeshDistributor->exitParallelization();
delete MeshDistributor::globalMeshDistributor;
#ifdef HAVE_PARALLEL_MTL4
if (mtl_environment)
delete mtl_environment;
......
......@@ -177,11 +177,9 @@ namespace AMDiS {
return vectorComponents[row].getStatus();
}
/** \brief
* Returns true, if for the given matrix component the row and the col FE spaces
* are equal. Note that this is also the case, if there is another aux FE space,
* that may be different from both, the row and the col FE spaces.
*/
/// Returns true, if for the given matrix component the row and the col FE spaces
/// are equal. Note that this is also the case, if there is another aux FE space,
/// that may be different from both, the row and the col FE spaces.
inline bool eqSpaces(int row, int col) const
{
int status = matrixComponents[row][col].getStatus();
......
......@@ -1495,11 +1495,12 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
nullify(grd1);
for (int j = 0; j < nBasFcts; j++) // #BasisFunctions
for (int j = 0; j < nBasFcts; j++) { // #BasisFunctions
for (int k = 0; k < parts; k++) // #edges (2d) or #faces (3d)
grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j];
}
for (int l=0; l < dow; l++) {
for (int l = 0; l < dow; l++) {
nullify(grdAtQPs[iq][l]);
for (int k = 0; k < parts; k++)
grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
......
......@@ -51,20 +51,16 @@ namespace AMDiS {
/// Protected constructor. Avoids instatiation of the basis class
ElInfo();
/** \brief
* Protected constructor. Avoids instatiation of the basis class.
* \param mesh pointer to the corresponding mesh.
*/
/// Protected constructor. Avoids instatiation of the basis class.
/// \param mesh pointer to the corresponding mesh.
ElInfo(Mesh *mesh);
public:
/// Virtual destructor because ElInfo is pure virtual.
virtual ~ElInfo();
/** \brief
* Assignement operator.
* \param rhs right hand side.
*/
/// Assignement operator.
/// \param rhs right hand side.
ElInfo& operator=(const ElInfo& rhs)
{
mesh = rhs.mesh;
......@@ -137,37 +133,29 @@ namespace AMDiS {
return iChild;
}
/** \brief
* Get ElInfo's \ref coord[i]. This is a WorldVector<double> filled with the
* coordinates of the i-th vertex of element \ref el.
*/
/// Get ElInfo's \ref coord[i]. This is a WorldVector<double> filled with
/// the coordinates of the i-th vertex of element \ref el.
inline WorldVector<double>& getCoord(int i)
{
return coord[i];
}
/** \brief
* Get ElInfo's \ref coord[i]. This is a WorldVector<double> filled with the
* coordinates of the i-th vertex of element \ref el.
*/
/// Get ElInfo's \ref coord[i]. This is a WorldVector<double> filled with the
/// coordinates of the i-th vertex of element \ref el.
inline const WorldVector<double>& getCoord(int i) const
{
return coord[i];
}
/** \brief
* Get ElInfo's \ref coord. This is a FixVec<WorldVector<double> > filled with the
* coordinates of the all vertice of element \ref el.
*/
/// Get ElInfo's \ref coord. This is a FixVec<WorldVector<double> > filled
/// with the coordinates of the all vertice of element \ref el.
inline FixVec<WorldVector<double>, VERTEX>& getCoords()
{
return coord;
}
/** \brief
* Get ElInfo's \ref coord. This is a FixVec<WorldVector<double> > filled with the
* coordinates of the all vertice of element \ref el.
*/
/// Get ElInfo's \ref coord. This is a FixVec<WorldVector<double> > filled
/// with the coordinates of the all vertice of element \ref el.
inline const FixVec<WorldVector<double>, VERTEX>& getCoords() const
{
return coord;
......@@ -370,83 +358,67 @@ namespace AMDiS {
/** \} */
/** \brief
* Returns the absolute value of the determinant of the affine linear
* parametrization's Jacobian
*/
/// Returns the absolute value of the determinant of the affine linear
/// parametrization's Jacobian
virtual double calcDet() const;
/** \brief
* Used by non static method \ref calcDet(). Calculates the determinant
* for a given vector of vertex coordinates.
*/
/// Used by non static method \ref calcDet(). Calculates the determinant
/// for a given vector of vertex coordinates.
double calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const;
/// from CFE_Integration
double calcSurfaceDet(VectorOfFixVecs<DimVec<double> > &surfVert) const;
/// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the program exits.
/// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the
/// program exits.
void testFlag(const Flag& flag) const;
/** \brief
* Transforms local barycentric coordinates of a point defined on this element
* to global world coordinates.
*/
/// Transforms local barycentric coordinates of a point defined on this
/// element to global world coordinates.
void coordToWorld(const DimVec<double>& lambda,
WorldVector<double>& world) const;
/// Fills ElInfo's \ref det and \ref grdLambda entries.
virtual void fillDetGrdLambda();
/** \brief
* Returns a pointer to a vector, which contains the barycentric coordinates
* with respect to \ref element of a point with world coordinates world.
* The barycentric coordinates are stored in lambda.
* pure virtual => must be overriden in sub-class.
*/
/// Returns a pointer to a vector, which contains the barycentric coordinates
/// with respect to \ref element of a point with world coordinates world.
/// The barycentric coordinates are stored in lambda.
/// pure virtual => must be overriden in sub-class.
virtual const int worldToCoord(const WorldVector<double>& world,
DimVec<double>* lambda) const = 0;
/** \brief
* Fills this ElInfo with macro element information of mel.
* pure virtual => must be overriden in sub-class.
*/
/// Fills this ElInfo with macro element information of mel.
/// pure virtual => must be overriden in sub-class.
virtual void fillMacroInfo(const MacroElement *mel) = 0;
/** \brief
* Fills this ElInfo for the child ichild using hierarchy information and
* parent data parentInfo.
* pure virtual => must be overriden in sub-class.
*/
/// Fills this ElInfo for the child ichild using hierarchy information and
/// parent data parentInfo.
/// pure virtual => must be overriden in sub-class.
virtual void fillElInfo(int ichild, const ElInfo *parentInfo) = 0;
void fillElInfo(const MacroElement *mel, int refinementPathLength, unsigned long refinementPath);
void fillElInfo(const MacroElement *mel,
int refinementPathLength,
unsigned long refinementPath);
/** \brief
* Calculates the Jacobian of the barycentric coordinates on \element and
* stores the matrix in grd_lam. The return value of the function is the
* absolute value of the determinant of the affine linear paraetrization's
* Jacobian.
* pure virtual => must be overriden in sub-class.
*/
/// Calculates the Jacobian of the barycentric coordinates on \element and
/// stores the matrix in grd_lam. The return value of the function is the
/// absolute value of the determinant of the affine linear paraetrization's
/// Jacobian.
/// pure virtual => must be overriden in sub-class.
virtual double calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) = 0;
/** \brief
* calculates a normal of the given side (1d, 2d: edge, 3d: face) of \ref element.
* Returns the absolute value of the determinant of the
* transformation to the reference element.
* pure virtual => must be overriden in sub-class.
*/
/// calculates a normal of the given side (1d, 2d: edge, 3d: face) of \ref element.
/// Returns the absolute value of the determinant of the
/// transformation to the reference element.
/// pure virtual => must be overriden in sub-class.
virtual double getNormal(int side, WorldVector<double> &normal) const = 0;
/** \brief
* calculates a normal of the element in dim of world = dim + 1.
* Returns the absolute value of the determinant of the
* transformation to the reference element.
* pure virtual => must be overriden in sub-class.
*/
/// calculates a normal of the element in dim of world = dim + 1.
/// Returns the absolute value of the determinant of the
/// transformation to the reference element.
/// pure virtual => must be overriden in sub-class.
virtual double getElementNormal(WorldVector<double> &elementNormal) const
{
FUNCNAME("ElInfo::getElementNormal()");
......
......@@ -63,7 +63,6 @@ namespace AMDiS {
}
FiniteElemSpace::FiniteElemSpace()
{}
......@@ -104,6 +103,14 @@ namespace AMDiS {
}
void FiniteElemSpace::destroyFeSpaces()
{
for (unsigned int i = 0; i < feSpaces.size(); i++)
delete feSpaces[i];
feSpaces.resize(0);
}
#if DEBUG
FiniteElemSpace* FiniteElemSpace::provideFeSpace(Mesh *mesh)
{
......
......@@ -51,6 +51,8 @@ namespace AMDiS {
Mesh *mesh,
string name = "");
static void destroyFeSpaces();
#if DEBUG
/// For debugging it may be useful to get some FE space for a given mesh at a
/// position in code where it is not possible to access the FE space directly. The
......@@ -104,10 +106,8 @@ namespace AMDiS {
getHighest(vector<const FiniteElemSpace*>& feSpaces);
protected:
/** \brief
* Constructs a FiniteElemSpace with name name_ and the given DOFAdmin,
* BasisFunction and Mesh.
*/
/// Constructs a FiniteElemSpace with name name_ and the given DOFAdmin,
/// BasisFunction and Mesh.
FiniteElemSpace(DOFAdmin* admin,
const BasisFunction* basisFcts,
Mesh* mesh,
......
......@@ -99,7 +99,7 @@ namespace AMDiS {
SubAssembler* subAssembler,
Quadrature *quad)
{
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
}
void CoordsAtQP_FOT::getLb(const ElInfo *elInfo,
......@@ -140,7 +140,7 @@ namespace AMDiS {
SubAssembler* subAssembler,
Quadrature *quad)
{
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
}
void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo,
......@@ -281,7 +281,7 @@ namespace AMDiS {
SubAssembler* subAssembler,
Quadrature *quad)
{
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
}
......@@ -477,7 +477,7 @@ namespace AMDiS {
Quadrature *quad)
{
getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
}
void FctVecAtQP_FOT::getLb(const ElInfo *elInfo,
......@@ -686,7 +686,7 @@ namespace AMDiS {
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
for (int i = 0; i < nVecs; i++)
getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
......@@ -713,7 +713,7 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
lb(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), Lb[iq], 1.0);
lb(grdLambda, (*f_)(coordsAtQPs[iq], vecsArg, gradsArg), Lb[iq], 1.0);
}
}
......@@ -741,7 +741,7 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
const WorldVector<double> &b = (*f_)(coordsAtQPs[iq], vecsArg, gradsArg);
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i] * b[i];
......@@ -788,7 +788,7 @@ namespace AMDiS {
int nGrads = static_cast<int>(grads_.size());
elInfo->getElementNormal(elementNormal_);
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
for (int i = 0; i < nVecs; i++)
getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
......@@ -815,18 +815,18 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
lb(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
lb(grdLambda, (*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg),
Lb[iq], 1.0);
}
}
void GeneralParametric_FOT::eval( int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double factor)
void GeneralParametric_FOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double factor)
{
int dow = Global::getGeo(WORLD);
int nVecs = static_cast<int>(vecs_.size());
......@@ -845,7 +845,7 @@ namespace AMDiS {
gradsArg[i] = gradsAtQPs_[i][iq];
const WorldVector<double> &b =
(*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
(*f_)(coordsAtQPs[iq], elementNormal_, vecsArg, gradsArg);
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i] * b[i];
......
......@@ -336,7 +336,7 @@ namespace AMDiS {
protected:
/// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
/// Function avaluated at world coordinates.
AbstractFunction<double, WorldVector<double> > *g;
......@@ -375,7 +375,7 @@ namespace AMDiS {
protected:
/// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
/// Function evaluated at world coordinates.
AbstractFunction<double, WorldVector<double> > *g;
......@@ -495,7 +495,7 @@ namespace AMDiS {
protected:
/// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
/// Function avaluated at world coordinates.
AbstractFunction<WorldVector<double>, WorldVector<double> > *g;
......@@ -657,9 +657,13 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
mtl::dense_vector<double> vecAtQPs;
WorldVector<double> *coordsAtQPs;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
BinaryAbstractFunction<double, WorldVector<double>, double> *f;
WorldVector<double> *b;
};
......@@ -772,7 +776,7 @@ namespace AMDiS {
vector<double>,
vector<WorldVector<double> > > *f_;
WorldVector<double> *coordsAtQPs_;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
vector<mtl::dense_vector<double> > vecsAtQPs;
......@@ -818,7 +822,8 @@ namespace AMDiS {
vector<double>,
vector<WorldVector<double> > > *f_;
WorldVector<double> *coordsAtQPs_;
mtl::dense_vector<WorldVector<double> > coordsAtQPs;
WorldVector<double> elementNormal_;
vector<mtl::dense_vector<double> > vecsAtQPs;
......
......@@ -53,20 +53,16 @@ namespace AMDiS {
{
public:
/** \brief
* constructor without initialisation. initType must be NO_INIT. If dim is
* not spezified, a FixVec for DIM_OF_WORLD is created.
*/
/// Constructor without initialisation. initType must be NO_INIT. If dim is
/// not spezified, a FixVec for DIM_OF_WORLD is created.
FixVec(int dim = -1, InitType initType = NO_INIT)
: Vector<T>(calcSize(dim))
{
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or missing initializer\n");
}
/** \brief
* constructor with value list initialisation. initType must be VALUE_LIST.
* ini is an array which contains the initialisation values.
*/
/// constructor with value list initialisation. initType must be VALUE_LIST.
/// ini is an array which contains the initialisation values.
FixVec(int dim, InitType initType, const T* ini)
: Vector<T>(calcSize(dim))
{
......@@ -74,10 +70,8 @@ namespace AMDiS {
setValues(ini);
}
/** \brief
* constructor with default value initialisation. initType must be
* DEFAULT_VALUE. All vector entries are set to ini.
*/
/// constructor with default value initialisation. initType must be
/// DEFAULT_VALUE. All vector entries are set to ini.
FixVec(int dim, InitType initType, const T& ini)
: Vector<T>(calcSize(dim))
{
......@@ -131,11 +125,9 @@ namespace AMDiS {
class VectorOfFixVecs
{
public:
/** \brief
* constructs a VectorOfFixVecs without initialisation. dim is passed to
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be NO_INIT.
*/
/// constructs a VectorOfFixVecs without initialisation. dim is passed to
/// FixVec's constructors. size_ is the number of contained FixVecs. initType
/// must be NO_INIT.
VectorOfFixVecs(int d, int s, InitType initType)
: size(s),
dim(d)
......@@ -147,11 +139,9 @@ namespace AMDiS {
vec[i] = new FixVecType(dim, NO_INIT);
}
/** \brief
* constructs a VectorOfFixVecs via an value list. dim is passed to
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be VALUE_LIST. ini contains the initialisation values.
*/
/// constructs a VectorOfFixVecs via an value list. dim is passed to
/// FixVec's constructors. size_ is the number of contained FixVecs. initType
/// must be VALUE_LIST. ini contains the initialisation values.
VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini)
: size(s),
dim(d)
......@@ -163,11 +153,9 @@ namespace AMDiS {
vec[i] = new FixVecType(ini[i]);
}
/** \brief
* constructs a VectorOfFixVecs with an default value. dim is passed to
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be DEFAULT_VALUE. All entries are set to ini.
*/
/// constructs a VectorOfFixVecs with an default value. dim is passed to
/// FixVec's constructors. size_ is the number of contained FixVecs. initType
/// must be DEFAULT_VALUE. All entries are set to ini.
VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini)
: size(s),
dim(d)
......@@ -273,11 +261,9 @@ namespace AMDiS {
class MatrixOfFixVecs
{
public:
/** \brief
* Constructs the matrix without initialisation. r is the number of rows,
* c is the number of columns. The other parameters are analog to the
* VectorOfFixVecs constructors.
*/
/// Constructs the matrix without initialisation. r is the number of rows,
/// c is the number of columns. The other parameters are analog to the
/// VectorOfFixVecs constructors.
MatrixOfFixVecs(int dim, int r, int c, InitType initType)
: rows(r),
columns(c)
......
......@@ -55,6 +55,7 @@ namespace AMDiS {
setFunctionPointer();
}
Lagrange::~Lagrange()
{
for (int i = 0; i < static_cast<int>(bary->size()); i++)
......@@ -64,6 +65,7 @@ namespace AMDiS {
}
}
Lagrange* Lagrange::getLagrange(int dim, int degree)
{
std::list<Lagrange*>::iterator it;
......@@ -76,6 +78,7 @@ namespace AMDiS {
return newLagrange;
}
void Lagrange::clear()
{
for (std::list<Lagrange*>::iterator it = allBasFcts.begin();
......@@ -86,6 +89,7 @@ namespace AMDiS {
}
}
void Lagrange::setFunctionPointer()
{
if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
......@@ -193,6 +197,7 @@ namespace AMDiS {
}
}
void Lagrange::setNDOF()
{
if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
......@@ -781,9 +786,9 @@ namespace AMDiS {
return NULL;
}
void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
void Lagrange::getBound(const ElInfo* elInfo, BoundaryType* bound) const
{
el_info->testFlag(Mesh::FILL_BOUND);
elInfo->testFlag(Mesh::FILL_BOUND);
// boundaries
int index = 0;
......@@ -795,7 +800,7 @@ namespace AMDiS {
for (int i = 0; i < dim; i++) {
int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
for (int j = offset; j < jto; j++) {
boundaryType = el_info->getBoundary(j);
boundaryType = elInfo->getBoundary(j);
int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
for (int k = 0; k < kto; k++)
bound[index++] = boundaryType;
......@@ -811,7 +816,7 @@ namespace AMDiS {
}
void Lagrange::interpol(const ElInfo *el_info,
void Lagrange::interpol(const ElInfo *elInfo,
int no,
const int *b_no,
AbstractFunction<double, WorldVector<double> > *f,
......@@ -821,32 +826,32 @@ namespace AMDiS {
WorldVector<double> x;
el_info->testFlag(Mesh::FILL_COORDS);
elInfo->testFlag(Mesh::FILL_COORDS);
if (b_no) {
TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n");
for (int i = 0; i < no; i++) {
if (b_no[i] < Global::getGeo(VERTEX, dim)) {
rvec[i] = (*f)(el_info->getCoord(b_no[i]));
rvec[i] = (*f)(elInfo->getCoord(b_no[i]));
} else {
el_info->coordToWorld(*(*bary)[b_no[i]], x);
elInfo->coordToWorld(*(*bary)[b_no[i]], x);
rvec[i] = (*f)(x);
}
}
} else {
int vertices = Global::getGeo(VERTEX, dim);
for (int i = 0; i < vertices; i++)
rvec[i] = (*f)(el_info->getCoord(i));
rvec[i] = (*f)(elInfo->getCoord(i));
for (int i = vertices; i < nBasFcts; i++) {
el_info->coordToWorld(*(*bary)[i], x);
elInfo->coordToWorld(*(*bary)[i], x);
rvec[i] = (*f)(x);
}
}
}
void Lagrange::interpol(const ElInfo *el_info,
void Lagrange::interpol(const ElInfo *elInfo,
int no,
const int *b_no,