-
Thomas Witkowski authoredThomas Witkowski authored
SystemVector.h 13.76 KiB
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// ============================================================================
// == ==
// == crystal growth group ==
// == ==
// == Stiftung caesar ==
// == Ludwig-Erhard-Allee 2 ==
// == 53175 Bonn ==
// == germany ==
// == ==
// ============================================================================
// == ==
// == http://www.caesar.de/cg/AMDiS ==
// == ==
// ============================================================================
/** \file SystemVector.h */
#ifndef AMDIS_SYSTEMVECTOR_H
#define AMDIS_SYSTEMVECTOR_H
#include "MemoryManager.h"
#include "MatrixVector.h"
#include "DOFVector.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "OpenMP.h"
namespace AMDiS {
// ============================================================================
// ===== class SystemVector ===================================================
// ============================================================================
/** \brief
* A system vector is a vector of dof vectors used for vector valued
* problems.
*/
class SystemVector : public Serializable
{
public:
MEMORY_MANAGED(SystemVector);
/** \brief
* Creator.
*/
class Creator
: public CreatorInterface<SystemVector> {
public:
MEMORY_MANAGED(Creator);
/** \brief
* Creators constructor.
*/
Creator(const ::std::string &name_,
::std::vector<FiniteElemSpace*> feSpace_,
int size_)
: name(name_),
feSpace(feSpace_),
size(size_)
{};
/** \brief
* Destructor
*/
virtual ~Creator() {};
/** \brief
* Implementation of CreatorInterface::create().
*/
SystemVector *create() {
char number[10];
::std::string numberedName;
SystemVector *result = NEW SystemVector(name, feSpace, size);
for (int i = 0; i < size; i++) {
sprintf(number, "[%d]", i);
numberedName = name + ::std::string(number);
result->setDOFVector(i, NEW DOFVector<double>(feSpace[i], numberedName));
}
return result;
};
/** \brief
* Implementation of CreatorInterface::free().
*/
void free(SystemVector *vec) {
for (int i = 0; i < size; i++) {
DELETE vec->getDOFVector(i);
vec->setDOFVector(i, NULL);
}
DELETE vec;
};
private:
/** \brief
* Name of the system vector
*/
::std::string name;
/** \brief
* Finite element space used for creation.
*/
::std::vector<FiniteElemSpace*> feSpace;
/** \brief
* Number of component vectors to be created
*/
int size;
};
public:
/** \brief
* Constructor.
*/
SystemVector(const ::std::string& name_,
::std::vector<FiniteElemSpace*> feSpace_,
int size)
: name(name_),
feSpace(feSpace_),
vectors(size)
{
vectors.set(NULL);
};
/** \brief
* Copy Constructor.
*/
SystemVector(const SystemVector& rhs)
: name(rhs.name),
feSpace(rhs.feSpace),
vectors(rhs.vectors.getSize())
{
for (int i = 0; i < vectors.getSize(); i++) {
vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
} };
virtual ~SystemVector() {};
/** \brief
* Sets \ref vectors[index] = vec.
*/
inline void setDOFVector(int index, DOFVector<double> *vec) {
TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
vectors[index] = vec;
};
/** \brief
* Returns \ref vectors[index].
*/
inline DOFVector<double> *getDOFVector(int index) {
TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
return vectors[index];
};
/** \brief
* Returns \ref vectors[index].
*/
inline const DOFVector<double> *getDOFVector(int index) const {
TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
return vectors[index];
};
/** \brief
* Returns sum of used vector sizes.
*/
inline int getUsedSize() const {
int totalSize = 0;
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
totalSize += vectors[i]->getUsedSize();
}
return totalSize;
};
/** \brief
* Returns number of contained vectors.
*/
inline int getNumVectors() const {
return vectors.getSize();
};
inline int getSize() const {
return vectors.getSize();
};
/** \brief
* Returns \ref feSpace.
*/
inline FiniteElemSpace *getFESpace(int i) const {
return feSpace[i];
};
/** \brief
* Here the system vector is interpreted as one large vector. The given
* is used as a global index which indicates a local vector number and
* a local index on this vector. The operator returns this local vector
* at the local index.
*/
inline double& operator[](int index) {
int localIndex = index;
int vectorIndex = 0;
while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
localIndex -= vectors[vectorIndex++]->getUsedSize(); }
return (*(vectors[vectorIndex]))[localIndex];
};
/** \brief
* For const access.
*/
inline double operator[](int index) const {
int localIndex = index;
int vectorIndex = 0;
while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
localIndex -= vectors[vectorIndex++]->getUsedSize();
}
return (*(vectors[vectorIndex]))[localIndex];
};
/** \brief
* Sets all entries in all vectors to value.
*/
inline void set(double value) {
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
vectors[i]->set(value);
}
};
/** \brief
* Sets all entries in all vectors to value.
*/
inline SystemVector& operator=(double value) {
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
(*(vectors[i])) = value;
}
return *this;
};
/** \brief
* Assignement operator.
*/
inline SystemVector& operator=(const SystemVector& rhs) {
TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
(*(vectors[i])) = (*(rhs.getDOFVector(i)));
}
return *this;
};
void serialize(::std::ostream &out) {
int size = vectors.getSize();
out.write(reinterpret_cast<const char*>(&size), sizeof(int));
for (int i = 0; i < size; i++) {
vectors[i]->serialize(out);
}
};
void deserialize(::std::istream &in) {
int size, oldSize = vectors.getSize();
in.read(reinterpret_cast<char*>(&size), sizeof(int));
vectors.resize(size);
for (int i = oldSize; i < size; i++) {
vectors[i] = NEW DOFVector<double>(feSpace[i], "");
}
for (int i = 0; i < size; i++) {
vectors[i]->deserialize(in);
} };
void copy(const SystemVector& rhs) {
int size = vectors.getSize();
TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
for (int i = 0; i < size; i++) {
vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
}
};
void interpol(::std::vector<AbstractFunction<double, WorldVector<double> >*> *f) {
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
vectors[i]->interpol((*f)[i]);
}
};
void print() {
int size = vectors.getSize();
for (int i = 0; i < size; i++) {
vectors[i]->print();
}
};
protected:
/** \brief
* Name of the system vector
*/
::std::string name;
/** \brief
* Finite element space.
*/
::std::vector<FiniteElemSpace*> feSpace;
/** \brief
* Local dof vectors.
*/
Vector<DOFVector<double>*> vectors;
};
/** \brief
* multiplication with scalar
*/
inline const SystemVector& operator*=(SystemVector& x, double d) {
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
*(x.getDOFVector(i)) *= d;
}
return x;
};
/** \brief
* scalar product
*/
inline double operator*(SystemVector& x, SystemVector& y) {
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
double result = 0.0;
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
};
return result;
};
/** \brief
* addition of two system vectors */
inline const SystemVector& operator+=(SystemVector& x,
const SystemVector& y)
{
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
(*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
}
return x;
};
/**
* subtraction of two system vectors.
*/
inline const SystemVector& operator-=(SystemVector& x,
SystemVector& y)
{
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
(*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
}
return x;
};
/** \brief
* multiplication with a scalar
*/
inline SystemVector operator*(SystemVector& x, double d) {
SystemVector result = x;
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
(*(result.getDOFVector(i))) *= d;
}
return result;
};
/** \brief
* multiplication with a scalar
*/
inline SystemVector operator*(double d, SystemVector& x) {
SystemVector result = x;
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
(*(result.getDOFVector(i))) *= d;
}
return result;
};
/** \brief
* addition of two system vectors
*/
inline SystemVector operator+(const SystemVector& x,
const SystemVector& y)
{
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
SystemVector result = x;
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
(*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
}
return result;
};
/** \brief
* Calls SystemVector::set(). Used for solving. */
inline void set(SystemVector& x, double value) {
x.set(value);
};
/** \brief
* Calls SystemVector::set(). Used for solving.
*/
inline void setValue(SystemVector& x, double value) {
x.set(value);
};
/** \brief
* Norm of system vector.
*/
inline double norm(SystemVector* x) {
double result = 0.0;
int size = x->getNumVectors();
for (int i = 0; i < size; i++) {
result += x->getDOFVector(i)->squareNrm2();
}
return sqrt(result);
};
/** \brief
* L2 norm of system vector.
*/
inline double L2Norm(SystemVector* x) {
double result = 0.0;
int size = x->getNumVectors();
for (int i = 0; i < size; i++) {
result += x->getDOFVector(i)->L2NormSquare();
}
return sqrt(result);
};
/** \brief
* H1 norm of system vector.
*/
inline double H1Norm(SystemVector* x) {
double result = 0.0;
int size = x->getNumVectors();
for (int i = 0; i < size; i++) {
result += x->getDOFVector(i)->H1NormSquare();
}
return sqrt(result);
};
inline void mv(Matrix<DOFMatrix*> &matrix,
const SystemVector &x,
SystemVector &result,
bool add = false)
{
int size = x.getNumVectors();
int i;
TEST_EXIT_DBG(size == result.getNumVectors())("incompatible sizes\n");
TEST_EXIT_DBG(size == matrix.getNumRows())("incompatible sizes\n");
TEST_EXIT_DBG(size == matrix.getNumCols())("incompatible sizes\n");
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1) num_threads(size)
#endif
for (i = 0; i < size; i++) {
if (!add)
result.getDOFVector(i)->set(0.0);
for (int j = 0; j < size; j++) {
if (matrix[i][j]) {
mv<double>(NoTranspose, *(matrix[i][j]),
*(x.getDOFVector(j)),
*(result.getDOFVector(i)),
true);
}
}
}
};
/** \brief
* y = a*x + y;
*/
inline void axpy(double a, SystemVector& x, SystemVector& y)
{
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
int size = x.getNumVectors();
int i;
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1) num_threads(size)
#endif
for (i = 0; i < size; i++) {
axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
}
};
/** \brief
* y = x + a*y
*/
inline void xpay(double a, SystemVector& x, SystemVector& y)
{
TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
("invalid size\n");
int size = x.getNumVectors();
for (int i = 0; i < size; i++) {
xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
}
}
/** \brief
* Returns SystemVector::getUsedSize().
*/
inline int size(SystemVector* vec) {
return vec->getUsedSize();
};
inline void print(SystemVector* vec) {
vec->print();
};
}
#endif // AMDIS_SYSTEMVECTOR_H