Commit 55d33f24 authored by Thomas Witkowski's avatar Thomas Witkowski

Several bugfixes in multi mesh codes.

parent 0aba6ca3
This diff is collapsed.
......@@ -73,6 +73,7 @@ namespace AMDiS {
userMat += factor * elementMatrix;
}
void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
......@@ -110,9 +111,10 @@ namespace AMDiS {
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
mat = tmpMat;
mat = tmpMat;
}
if (firstOrderAssemblerGrdPsi) {
......@@ -130,8 +132,17 @@ namespace AMDiS {
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
#if 0
std::cout << "ASM MAT:" << std::endl;
std::cout << mat << std::endl;
std::cout << "MULT MAT:" << std::endl;
std::cout << m << std::endl;
#endif
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
//tmpMat = m * mat;
tmpMat = mat * m;
mat = tmpMat;
}
......@@ -139,6 +150,7 @@ namespace AMDiS {
userMat += factor * elementMatrix;
}
void Assembler::calculateElementVector(const ElInfo *elInfo,
ElementVector& userVec,
double factor)
......
......@@ -185,22 +185,6 @@ namespace AMDiS {
return vectorComponents[row].getAuxFESpace(0);
}
void fakeFESpace(const FiniteElemSpace *feOld,
const FiniteElemSpace *feNew)
{
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (matrixComponents[i][j].getAuxFESpace(0) == feOld) {
matrixComponents[i][j].setAuxFESpace(feNew, 0);
}
}
if (vectorComponents[i].getAuxFESpace(0) == feOld) {
vectorComponents[i].setAuxFESpace(feNew, 0);
}
}
}
protected:
int nComponents;
......
......@@ -186,6 +186,21 @@ namespace AMDiS {
}
}
using namespace mtl;
#if 0
std::cout << "----- PRINT MAT --------" << std::endl;
std::cout << elMat << std::endl;
std::cout << "rows: ";
for (int i = 0; i < rowIndices.size(); i++)
std::cout << rowIndices[i] << " ";
std::cout << std::endl;
std::cout << "cols: ";
for (int i = 0; i < colIndices.size(); i++)
std::cout << colIndices[i] << " ";
std::cout << std::endl;
#endif
for (int i = 0; i < nRow; i++) {
DegreeOfFreedom row = rowIndices[i];
......@@ -212,14 +227,17 @@ namespace AMDiS {
}
}
double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
{
return matrix[a][b];
}
void DOFMatrix::freeDOFContent(int index)
{}
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
{
FUNCNAME("DOFMatrix::assemble()");
......@@ -235,9 +253,10 @@ namespace AMDiS {
if (factor != 1.0)
elementMatrix *= factor;
addElementMatrix(elementMatrix, bound, elInfo, NULL);
addElementMatrix(elementMatrix, bound, elInfo, NULL);
}
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
Operator *op)
{
......@@ -254,6 +273,7 @@ namespace AMDiS {
addElementMatrix(elementMatrix, bound, elInfo, NULL);
}
void DOFMatrix::assemble(double factor,
ElInfo *rowElInfo, ElInfo *colElInfo,
ElInfo *smallElInfo, ElInfo *largeElInfo,
......@@ -272,12 +292,11 @@ namespace AMDiS {
} else {
std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin();
for (; it != operators.end(); ++it, ++factorIt) {
for (; it != operators.end(); ++it, ++factorIt)
(*it)->getElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
*factorIt ? **factorIt : 1.0);
}
if (factor != 1.0)
......@@ -286,10 +305,11 @@ namespace AMDiS {
addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);
}
void DOFMatrix::assemble2(double factor,
ElInfo *mainElInfo, ElInfo *auxElInfo,
ElInfo *smallElInfo, ElInfo *largeElInfo,
const BoundaryType *bound, Operator *op)
ElInfo *smallElInfo, ElInfo *largeElInfo,
const BoundaryType *bound, Operator *op)
{
FUNCNAME("DOFMatrix::assemble2()");
......@@ -419,27 +439,32 @@ namespace AMDiS {
}
}
void DOFMatrix::copy(const DOFMatrix& rhs)
{
matrix = rhs.matrix;
}
void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
{
inserter_type &ins= *inserter;
FUNCNAME("DOFMatrix::removeRowsWithDBC()");
inserter_type &ins = *inserter;
for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
ins[*it][*it] = 1.0;
ins[*it][*it] = 1.0;
rows->clear();
}
int DOFMatrix::memsize()
{
return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
+ matrix.nnz() * sizeof(base_matrix_type::value_type);
}
void DOFMatrix::startInsertion(int nnz_per_row)
{
if (inserter) {
......
......@@ -307,8 +307,9 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast)
if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
......@@ -708,16 +709,10 @@ namespace AMDiS {
DOFVector<DegreeOfFreedom>()
{}
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof)
{
/*
std::vector<DegreeOfFreedom>::iterator it;
std::vector<DegreeOfFreedom>::iterator end = vec.end();
DegreeOfFreedom pos = 0;
for (it = vec.begin(); it != end; ++it, ++pos)
if (*it == dof) *it = pos;
*/
}
{}
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *res)
......@@ -745,6 +740,7 @@ namespace AMDiS {
return r;
}
template<>
const double *DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
......@@ -765,7 +761,8 @@ namespace AMDiS {
("invalid basis functions");
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
int nPoints =
quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static double *localvec = NULL;
double *result;
......@@ -800,6 +797,7 @@ namespace AMDiS {
return const_cast<const double*>(result);
}
template<>
void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
const BoundaryType *bound,
......
......@@ -987,9 +987,10 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast)
if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
......
......@@ -44,12 +44,39 @@ namespace AMDiS {
int myRank = MPI::COMM_WORLD.Get_rank();
if (rank == -1 || myRank == rank) {
DOFVector<double> tmp(feSpace, "tmp");
VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu", false);
VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu",
false);
}
}
#endif
void writeDofIndexMesh(FiniteElemSpace *feSpace)
{
DOFVector<double> tmp(feSpace, "tmp");
DOFIterator<double> it(&tmp, USED_DOFS);
for (it.reset(); !it.end(); ++it)
*it = it.getDOFIndex();
VtkWriter::writeFile(tmp, "dofindex.vtu");
}
void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename)
{
std::map<int, double> vec;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
int index = elInfo->getElement()->getIndex();
vec[index] = index;
elInfo = stack.traverseNext(elInfo);
}
ElementFileWriter::writeFile(vec, feSpace, filename);
}
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el)
{
// === Get local indices of the given element. ===
......@@ -187,22 +214,6 @@ namespace AMDiS {
}
void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename)
{
std::map<int, double> vec;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
int index = elInfo->getElement()->getIndex();
vec[index] = index;
elInfo = stack.traverseNext(elInfo);
}
ElementFileWriter::writeFile(vec, feSpace, filename);
}
void writeMatlabMatrix(DOFMatrix &mat, std::string filename)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
......
......@@ -43,9 +43,36 @@ namespace AMDiS {
void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename);
/** \brief
* Writes a vtu file with the mesh, where all DOFs are set to zero, and only
* one given DOF is set to one. This can be used to easily identify DOFs in
* a mesh.
*
* \param[in] rank If set to -1, the vtu files are written on all ranks.
* Otherwise, only on the given rank the mesh is written.
* \param[in] dof Defines the DOF, which value is set to one in the mesh file.
* \param[in] feSpace The FE space to be used.
*/
void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace);
#endif
/** \brief
* Create a vtu file with name 'dofindex.vtu'. All nodes in the mesh are colored
* by the global dof index.
*
* \param[in] feSpace The FE space to be used.
*/
void writeDofIndexMesh(FiniteElemSpace *feSpace);
/** \brief
* Creates a vtu file where all elements in the mesh are colored by the global
* element indices.
*
* \param[in] feSpace The FE space to be used.
* \param[in] filename Name of the file.
*/
void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename);
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec,
......@@ -60,8 +87,6 @@ namespace AMDiS {
void printMatValuesStatistics(Matrix<DOFMatrix*> *mat);
void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename);
/** \brief
* Creates a text file storing the value of a sparse matrix. Each line of the file
* has three columns:
......
......@@ -294,10 +294,12 @@ namespace AMDiS {
}
}
void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat,
int basisFctsDegree) const
void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
{
switch (basisFctsDegree) {
FUNCNAME("ElInfo1d::getSubElemCoordsMat()");
switch (degree) {
case 1:
mat = mat_d1;
......@@ -323,14 +325,16 @@ namespace AMDiS {
break;
default:
ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
}
void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
int basisFctsDegree) const
void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
{
switch (basisFctsDegree) {
FUNCNAME("ElInfo1d::getSubElemCoordsMat_so()");
switch (degree) {
case 1:
mat = test2_val;
......@@ -339,7 +343,7 @@ namespace AMDiS {
break;
default:
ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
}
......
......@@ -62,11 +62,9 @@ namespace AMDiS {
return (i + 1) % 2;
}
void getSubElemCoordsMat(mtl::dense2D<double>& mat,
int basisFctsDegree) const;
void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
int basisFctsDegree) const;
void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
protected:
static double mat_d1_val[2][2];
......
......@@ -19,16 +19,31 @@ namespace AMDiS {
{0.0, 0.0, 1.0}};
mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val);
/*
double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5},
{0.0, 0.0, 0.5},
{1.0, 0.0, 0.0}};
{0.0, 0.0, 0.5},
{1.0, 0.0, 0.0}};
*/
double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 0.0, 1.0},
{1.0, 0.0, 0.0},
{0.5, 0.5, 0.0}};
mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);
/*
double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5},
{1.0, 0.0, 0.5},
{0.0, 1.0, 0.0}};
{1.0, 0.0, 0.5},
{0.0, 1.0, 0.0}};
*/
double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0},
{0.5, 0.5, 0.0}};
mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);
ElInfo2d::ElInfo2d(Mesh *aMesh)
: ElInfo(aMesh)
{
......@@ -37,6 +52,7 @@ namespace AMDiS {
normal = new WorldVector<double>;
}
ElInfo2d::~ElInfo2d()
{
delete e1;
......@@ -44,6 +60,7 @@ namespace AMDiS {
delete normal;
}
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
{
FUNCNAME("ElInfo::fillMacroInfo()");
......@@ -571,7 +588,8 @@ namespace AMDiS {
return adet;
}
const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
DimVec<double>* lambda) const
{
FUNCNAME("ElInfo::worldToCoord()");
......@@ -654,6 +672,7 @@ namespace AMDiS {
return det;
}
/****************************************************************************/
/* calculate the normal of the element for dim of world = dim + 1 */
/* return the absulute value of the determinant from the */
......@@ -680,4 +699,43 @@ namespace AMDiS {
return det;
}
void ElInfo2d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
{
FUNCNAME("ElInfo2d::getSubElemCoordsMat()");
using namespace mtl;
switch (degree) {
case 1:
{
mat = mat_d1;
dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
for (int i = 0; i < refinementPathLength; i++) {
if (refinementPath & (1 << i)) {
tmpMat = mat_d1_right * mat;
mat = tmpMat;
// mat *= mat_d1_right;
} else {
tmpMat = mat_d1_left * mat;
mat = tmpMat;
// mat *= mat_d1_left;
}
}
}
break;
default:
ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
}
void ElInfo2d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
{
FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");
ERROR_EXIT("Not yet implemented!\n");
}
}
......@@ -58,6 +58,10 @@ namespace AMDiS {
/// 2-dimensional realisation of ElInfo's getElementNormal method.
double getElementNormal(WorldVector<double> &normal) const;
void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
protected:
/// Temp vectors for function \ref calcGrdLambda.
WorldVector<double> *e1, *e2, *normal;
......
......@@ -590,67 +590,20 @@ namespace AMDiS {
}
}
#if 0
void ElInfo3d::getRefSimplexCoords(const BasisFunction *basisFcts,
mtl::dense2D<double>& coords) const
void ElInfo3d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
{
// (*coords)[0][0] = 1.0;
// (*coords)[1][0] = 0.0;
// (*coords)[2][0] = 0.0;
// (*coords)[3][0] = 0.0;
// (*coords)[0][1] = 0.0;
// (*coords)[1][1] = 1.0;
// (*coords)[2][1] = 0.0;
// (*coords)[3][1] = 0.0;
// (*coords)[0][2] = 0.0;
// (*coords)[1][2] = 0.0;
// (*coords)[2][2] = 1.0;
// (*coords)[3][2] = 0.0;
// (*coords)[0][3] = 0.0;
// (*coords)[1][3] = 0.0;
// (*coords)[2][3] = 0.0;
// (*coords)[3][3] = 1.0;
FUNCNAME("ElInfo3d::getSubElemCoordsMat()");
ERROR_EXIT("Not yet implemented!\n");
}
void ElInfo3d::getSubElementCoords(const BasisFunction *basisFcts,
int iChild,
mtl::dense2D<double>& coords) const
void ElInfo3d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
{
FUNCNAME("ElInfo3d::getSubElementCoords()");
ERROR_EXIT("Not yet\n");
// double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
// double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
// double c2 = ((*coords)[2][0] + (*coords)[2][1]) * 0.5;
// double c3 = ((*coords)[3][0] + (*coords)[3][1]) * 0.5;
// if (iChild == 0) {
// (*coords)[0][1] = (*coords)[0][0];
// (*coords)[1][1] = (*coords)[1][0];
// (*coords)[2][1] = (*coords)[2][0];
// (*coords)[0][0] = (*coords)[0][2];
// (*coords)[1][0] = (*coords)[1][2];
// (*coords)[2][0] = (*coords)[2][2];
// } else {
// (*coords)[0][1] = (*coords)[0][2];
// (*coords)[1][1] = (*coords)[1][2];
// (*coords)[2][1] = (*coords)[2][2];
// (*coords)[0][0] = (*coords)[0][1];
// (*coords)[1][0] = (*coords)[1][1];
// (*coords)[2][0] = (*coords)[2][1];
// }
// (*coords)[0][3] = c0;
// (*coords)[1][3] = c1;
// (*coords)[2][3] = c2;
// (*coords)[3][3] = c3;
FUNCNAME("ElInfo3d::getSubElemCoordsMat_so()");
ERROR_EXIT("Not yet implemented!\n");
}
#endif
}
......@@ -80,6 +80,10 @@ namespace AMDiS {
orientation = o;
}
void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
protected:
/** \brief
* +/- 1: sign of the determinant of the transformation to the reference
......
......@@ -91,6 +91,7 @@ namespace AMDiS {
}
}