Commit 4dea0e15 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Better solution to create element matrices for assembling on different meshes.

parent 6caf6e9c
...@@ -16,7 +16,6 @@ if USE_PARALLEL_AMDIS ...@@ -16,7 +16,6 @@ if USE_PARALLEL_AMDIS
$(PARALLEL_DIR)/MeshStructure_ED.h \ $(PARALLEL_DIR)/MeshStructure_ED.h \
$(PARALLEL_DIR)/ParallelError.h $(PARALLEL_DIR)/ParallelError.hh \ $(PARALLEL_DIR)/ParallelError.h $(PARALLEL_DIR)/ParallelError.hh \
$(PARALLEL_DIR)/ParallelProblem.h $(PARALLEL_DIR)/ParallelProblem.cc \ $(PARALLEL_DIR)/ParallelProblem.h $(PARALLEL_DIR)/ParallelProblem.cc \
$(PARALLEL_DIR)/ParallelDomainProblem.h $(PARALLEL_DIR)/ParallelDomainProblem.cc \
$(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \ $(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \
$(PARALLEL_DIR)/PartitionElementData.h $(PARALLEL_DIR)/PartitionElementData.h
PARALLEL_INCLUDES = -I/work/home7/witkowsk/local/include \ PARALLEL_INCLUDES = -I/work/home7/witkowsk/local/include \
......
...@@ -109,22 +109,35 @@ namespace AMDiS { ...@@ -109,22 +109,35 @@ namespace AMDiS {
ElementMatrix& mat = rememberElMat ? elementMatrix : userMat; ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
TEST_EXIT(zeroOrderAssembler)("Not yet implemented for not zero order assembler!\n"); if (secondOrderAssembler) {
if (secondOrderAssembler)
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat); secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
if (firstOrderAssemblerGrdPsi)
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
mat = tmpMat;
}
if (firstOrderAssemblerGrdPsi) {
ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat); firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
if (firstOrderAssemblerGrdPhi) }
if (firstOrderAssemblerGrdPhi) {
ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat); firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
if (zeroOrderAssembler) }
zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat(); if (zeroOrderAssembler) {
ElementMatrix tmpMat(nRow, nRow); zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
tmpMat = m * mat;
mat = tmpMat;
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
mat = tmpMat;
}
if (rememberElMat) if (rememberElMat)
userMat += factor * elementMatrix; userMat += factor * elementMatrix;
...@@ -240,23 +253,20 @@ namespace AMDiS { ...@@ -240,23 +253,20 @@ namespace AMDiS {
FUNCNAME("Assembler::matVecAssemble()"); FUNCNAME("Assembler::matVecAssemble()");
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
const BasisFunction *basFcts = rowFESpace->getBasisFcts(); double *uhOldLoc = new double[nRow];
int nBasFcts = basFcts->getNumber();
double *uhOldLoc = new double[nBasFcts];
operat->uhOld->getLocalVector(el, uhOldLoc); operat->uhOld->getLocalVector(el, uhOldLoc);
if (el != lastMatEl) if (el != lastMatEl)
calculateElementMatrix(elInfo, elementMatrix); calculateElementMatrix(elInfo, elementMatrix);
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nRow; i++) {
double val = 0.0; double val = 0.0;
for (int j = 0; j < nBasFcts; j++) { for (int j = 0; j < nRow; j++) {
val += elementMatrix[i][j] * uhOldLoc[j]; val += elementMatrix[i][j] * uhOldLoc[j];
} }
vec[i] += val; vec[i] += val;
} }
delete [] uhOldLoc; delete [] uhOldLoc;
} }
...@@ -283,7 +293,8 @@ namespace AMDiS { ...@@ -283,7 +293,8 @@ namespace AMDiS {
operat->uhOld->getLocalVector(auxEl, uhOldLoc); operat->uhOld->getLocalVector(auxEl, uhOldLoc);
const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat(); mtl::dense2D<double> m(nRow, nRow);
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
uhOldLoc2[i] = 0.0; uhOldLoc2[i] = 0.0;
...@@ -306,7 +317,6 @@ namespace AMDiS { ...@@ -306,7 +317,6 @@ namespace AMDiS {
vec[i] += val; vec[i] += val;
} }
delete [] uhOldLoc; delete [] uhOldLoc;
delete [] uhOldLoc2; delete [] uhOldLoc2;
} }
......
...@@ -72,14 +72,13 @@ namespace AMDiS { ...@@ -72,14 +72,13 @@ namespace AMDiS {
virtual void reset() { virtual void reset() {
position = 0; position = 0;
dofFreeIterator = dofFree->begin(); dofFreeIterator = dofFree->begin();
if(dofFreeIterator == dofFree->end()) { if (dofFreeIterator == dofFree->end())
return; return;
}
goToBeginOfIteratedObject(); goToBeginOfIteratedObject();
if(type != ALL_DOFS) { if (type != ALL_DOFS)
if(*dofFreeIterator == (type == USED_DOFS)) if (*dofFreeIterator == (type == USED_DOFS))
operator++(); operator++();
}
} }
/** \brief /** \brief
...@@ -196,14 +195,20 @@ namespace AMDiS { ...@@ -196,14 +195,20 @@ namespace AMDiS {
virtual void decObjectIterator() {} virtual void decObjectIterator() {}
protected: protected:
DOFAdmin *dofAdmin; /**< DOFAdmin which contains /// DOFAdmin which contains the dofFree vector.
* the dofFree vector DOFAdmin *dofAdmin;
*/
int position; /**< current position index */ /// Current position index.
std::vector<bool> *dofFree; /**< stores which DOFs are used int position;
*/
std::vector<bool>::iterator dofFreeIterator; /**< iterator for dofFree */ /// Stores which DOFs are used.
const DOFIteratorType type; /**< type of this iterator */ std::vector<bool> *dofFree;
/// Iterator for dofFree.
std::vector<bool>::iterator dofFreeIterator;
/// Type of this iterator.
const DOFIteratorType type;
}; };
......
#include <boost/numeric/mtl/mtl.hpp> #include <boost/numeric/mtl/mtl.hpp>
#include "DOFVector.h" #include "DOFVector.h"
#include "Traverse.h" #include "Traverse.h"
#include "DualTraverse.h" #include "DualTraverse.h"
...@@ -314,10 +313,8 @@ namespace AMDiS { ...@@ -314,10 +313,8 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n"); TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast) { if (quad && quadFast)
TEST_EXIT_DBG(quad == quadFast->getQuadrature()) TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts()) TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions"); ("invalid basis functions");
...@@ -335,29 +332,29 @@ namespace AMDiS { ...@@ -335,29 +332,29 @@ namespace AMDiS {
delete [] grd; delete [] grd;
grd = new WorldVector<double>[nPoints]; grd = new WorldVector<double>[nPoints];
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++)
grd[i].set(0.0); grd[i].set(0.0);
}
result = grd; result = grd;
} }
double *localVec = localVectors[myRank]; double *localVec = localVectors[myRank];
getLocalVector(largeElInfo->getElement(), localVec); getLocalVector(largeElInfo->getElement(), localVec);
const mtl::dense2D<double> m = smallElInfo->getSubElemCoordsMat();
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> m(nBasFcts, nBasFcts);
smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
DimVec<double> &grd1 = *grdTmp[myRank]; DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim); int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda(); const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();
const BasisFunction *basFcts = feSpace->getBasisFcts();
DimVec<double>* grdPhi = grdPhis[myRank]; DimVec<double>* grdPhi = grdPhis[myRank];
DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0); DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < dim + 1; j++) { for (int j = 0; j < dim + 1; j++)
grd1[j] = 0.0; grd1[j] = 0.0;
}
for (int j = 0; j < nBasFcts; j++) { for (int j = 0; j < nBasFcts; j++) {
...@@ -368,16 +365,14 @@ namespace AMDiS { ...@@ -368,16 +365,14 @@ namespace AMDiS {
*grdPhi += tmp; *grdPhi += tmp;
} }
for (int k = 0; k < parts; k++) { for (int k = 0; k < parts; k++)
grd1[k] += (*grdPhi)[k] * localVec[j]; grd1[k] += (*grdPhi)[k] * localVec[j];
}
} }
for (int l = 0; l < dow; l++) { for (int l = 0; l < dow; l++) {
result[i][l] = 0.0; result[i][l] = 0.0;
for (int k = 0; k < parts; k++) { for (int k = 0; k < parts; k++)
result[i][l] += grdLambda[k][l] * grd1[k]; result[i][l] += grdLambda[k][l] * grd1[k];
}
} }
} }
...@@ -804,7 +799,8 @@ namespace AMDiS { ...@@ -804,7 +799,8 @@ namespace AMDiS {
double *localVec = localVectors[omp_get_thread_num()]; double *localVec = localVectors[omp_get_thread_num()];
getLocalVector(largeElInfo->getElement(), localVec); getLocalVector(largeElInfo->getElement(), localVec);
const mtl::dense2D<double> m = smallElInfo->getSubElemCoordsMat(); mtl::dense2D<double> m(nBasFcts, nBasFcts);
smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++) {
result[i] = 0.0; result[i] = 0.0;
......
...@@ -290,11 +290,11 @@ namespace AMDiS { ...@@ -290,11 +290,11 @@ namespace AMDiS {
{} {}
DOFVector<T> *create() { DOFVector<T> *create() {
return NEW DOFVector<T>(feSpace, ""); return new DOFVector<T>(feSpace, "");
} }
void free(DOFVector<T> *vec) { void free(DOFVector<T> *vec) {
DELETE vec; delete vec;
} }
private: private:
...@@ -493,12 +493,14 @@ namespace AMDiS { ...@@ -493,12 +493,14 @@ namespace AMDiS {
void interpol(DOFVector<T> *v, double factor); void interpol(DOFVector<T> *v, double factor);
/// Writes the data of the DOFVector to an output stream.
void serialize(std::ostream &out) { void serialize(std::ostream &out) {
unsigned int size = vec.size(); unsigned int size = vec.size();
out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int)); out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T)); out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
} }
/// Reads data of an DOFVector from an input stream.
void deserialize(std::istream &in) { void deserialize(std::istream &in) {
unsigned int size; unsigned int size;
in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int)); in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
......
...@@ -170,12 +170,19 @@ namespace AMDiS { ...@@ -170,12 +170,19 @@ namespace AMDiS {
if (!fillSubElemMat) if (!fillSubElemMat)
return; return;
mtl::dense2D<double>& subCoordsMat = // mtl::dense2D<double>& subCoordsMat =
const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat()); // const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
// mtl::dense2D<double>& subCoordsMat_so =
// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
if (elInfo1 == elInfoSmall) { if (elInfo1 == elInfoSmall) {
stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat); // stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
// stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
} else { } else {
stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat); // stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
// stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
} }
} }
} }
...@@ -29,7 +29,8 @@ namespace AMDiS { ...@@ -29,7 +29,8 @@ namespace AMDiS {
neighbourCoord_(mesh_->getDim(), NO_INIT), neighbourCoord_(mesh_->getDim(), NO_INIT),
oppVertex_(mesh_->getDim(), NO_INIT), oppVertex_(mesh_->getDim(), NO_INIT),
grdLambda(mesh_->getDim(), NO_INIT), grdLambda(mesh_->getDim(), NO_INIT),
subElemCoordsMat(2, 2) refinementPath(0),
refinementPathLength(0)
{ {
projection_.set(NULL); projection_.set(NULL);
...@@ -46,14 +47,12 @@ namespace AMDiS { ...@@ -46,14 +47,12 @@ namespace AMDiS {
WorldVector<double>& w) const WorldVector<double>& w) const
{ {
int dim = l.getSize() - 1;
double c = l[0]; double c = l[0];
for (int j = 0; j < dimOfWorld; j++) for (int j = 0; j < dimOfWorld; j++)
w[j] = c * coord_[0][j]; w[j] = c * coord_[0][j];
int vertices = Global::getGeo(VERTEX, dim); int vertices = Global::getGeo(VERTEX, l.getSize() - 1);
for (int i = 1; i < vertices; i++) { for (int i = 1; i < vertices; i++) {
c = l[i]; c = l[i];
...@@ -62,10 +61,6 @@ namespace AMDiS { ...@@ -62,10 +61,6 @@ namespace AMDiS {
} }
} }
/****************************************************************************/
/* compute volume of an element */
/****************************************************************************/
double ElInfo::calcDet() const double ElInfo::calcDet() const
{ {
testFlag(Mesh::FILL_COORDS); testFlag(Mesh::FILL_COORDS);
......
...@@ -205,13 +205,11 @@ namespace AMDiS { ...@@ -205,13 +205,11 @@ namespace AMDiS {
return parametric_; return parametric_;
} }
/** \brief virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
* Returns the element transformation matrix \ref subElemCoordsMat . int basisFctsDegree) const {}
* This value is set only during dual traverse.
*/ virtual void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
inline const mtl::dense2D<double>& getSubElemCoordsMat() const { int basisFctsDegree) const {}
return subElemCoordsMat;
}
/** \} */ /** \} */
...@@ -279,14 +277,18 @@ namespace AMDiS { ...@@ -279,14 +277,18 @@ namespace AMDiS {
parametric_ = param; parametric_ = param;
} }
/** \brief /// Set \ref refinementPath
* Sets the element transformation matrix \ref subElemCoordsMat . inline void setRefinementPath(unsigned long rPath)
* This value is used only during dual traverse. {
*/ refinementPath = rPath;
// inline void setSubElemCoordsMat(DimMat<double> *mat) { }
// subElemCoordsMat = mat;
// } /// Set \ref refinementPathLength
inline void setRefinementPathLength(int length)
{
refinementPathLength = length;
}
/** \} */ /** \} */
...@@ -372,14 +374,6 @@ namespace AMDiS { ...@@ -372,14 +374,6 @@ namespace AMDiS {
return(0.0); return(0.0);
} }
virtual void getRefSimplexCoords(const BasisFunction *basisFcts,
mtl::dense2D<double>& coords) const = 0;
virtual void getSubElementCoords(const BasisFunction *basisFcts,
int iChild,
mtl::dense2D<double>& coords) const = 0;
protected: protected:
/// Pointer to the current mesh /// Pointer to the current mesh
Mesh *mesh_; Mesh *mesh_;
...@@ -471,6 +465,10 @@ namespace AMDiS { ...@@ -471,6 +465,10 @@ namespace AMDiS {
/// Stores the world dimension. /// Stores the world dimension.
int dimOfWorld; int dimOfWorld;
unsigned long refinementPath;
int refinementPathLength;
/** \brief /** \brief
* This is a transformation matrix used during dual traverse. It is set, if * This is a transformation matrix used during dual traverse. It is set, if
* the current element is the smaller element of an element pair in the traverse. * the current element is the smaller element of an element pair in the traverse.
...@@ -480,6 +478,8 @@ namespace AMDiS { ...@@ -480,6 +478,8 @@ namespace AMDiS {
*/ */
mtl::dense2D<double> subElemCoordsMat; mtl::dense2D<double> subElemCoordsMat;
mtl::dense2D<double> subElemCoordsMat_so;
public: public:
/** \brief /** \brief
* child_vertex[el_type][child][i] = father's local vertex index of new * child_vertex[el_type][child][i] = father's local vertex index of new
......
...@@ -40,6 +40,10 @@ namespace AMDiS { ...@@ -40,6 +40,10 @@ namespace AMDiS {
{1.0, 0.75, 0.0}, {1.0, 0.75, 0.0},
{0.0, 0.375, 1.0}}; {0.0, 0.375, 1.0}};
mtl::dense2D<double> ElInfo1d::mat_d2_right(mat_d2_right_val); mtl::dense2D<double> ElInfo1d::mat_d2_right(mat_d2_right_val);
double ElInfo1d::test1_val[2][2] = {{1.0, 0.0},
{0.0, 1.0}};
mtl::dense2D<double> ElInfo1d::test2_val(test1_val);
void ElInfo1d::fillMacroInfo(const MacroElement * mel) void ElInfo1d::fillMacroInfo(const MacroElement * mel)
...@@ -122,11 +126,10 @@ namespace AMDiS { ...@@ -122,11 +126,10 @@ namespace AMDiS {
if (adet2 < 1.0E-15) { if (adet2 < 1.0E-15) {
MSG("det*det = %lf\n", adet2); MSG("det*det = %lf\n", adet2);
for (int i = 0; i <= 1; i++) grd_lam[0] = grd_lam[1] = 0.0;
grd_lam[i] = 0.0;
} else { } else {
grd_lam[1] = e * (1.0 / adet2); grd_lam[1] = e * (1.0 / adet2);
grd_lam[0] = grd_lam[1] * (-1.0); grd_lam[0] = grd_lam[1] * -1.0;
} }
return sqrt(adet2); return sqrt(adet2);
...@@ -296,48 +299,52 @@ namespace AMDiS { ...@@ -296,48 +299,52 @@ namespace AMDiS {
} }
} }
void ElInfo1d::getRefSimplexCoords(const BasisFunction *basisFcts, void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int basisFctsDegree) const
mtl::dense2D<double>& coords) const
{ {
FUNCNAME("ElInfo1d::getRefSimplexCoords()"); switch (basisFctsDegree) {
switch (basisFcts->getDegree()) {
case 1: case 1:
coords = mat_d1; mat = mat_d1;
for (int i = 0; i < refinementPathLength; i++) {
if (refinementPath & (1 << i))
mat *= mat_d1_left;
else
mat *= mat_d1_right;
}
break; break;
case 2: case 2:
coords = mat_d2; mat = mat_d2;