Commit 70789436 authored by Thomas Witkowski's avatar Thomas Witkowski

Work on periodic boundary conditions.

parent 085e86f9
......@@ -51,15 +51,13 @@ namespace AMDiS {
{
FUNCNAME("PeriodicDOFMapping::getDOFPermutation()");
if(dofPermutation_[vertexPermutation] == NULL) {
if (dofPermutation_[vertexPermutation] == NULL) {
//MSG("new dof permutation needed\n");
int i, j;
int dim = basFcts_->getDim();
int num = basFcts_->getNumber();
int sum = 0;
for(i = 0; i < dim + 1; i++) {
for (int i = 0; i < dim + 1; i++) {
sum += i - vertexPermutation[i];
TEST_EXIT(vertexPermutation[i] < dim + 1)
("invalid vertexPermuation\n");
......@@ -72,11 +70,11 @@ namespace AMDiS {
DegreeOfFreedom *mapping = new DegreeOfFreedom[num];
for(i = 0; i < num; i++) {
for (int i = 0; i < num; i++) {
lambda = basFcts_->getCoords(i);
for(j = 0; j < dim + 1; j++) {
for (int j = 0; j < dim + 1; j++)
newLambda[vertexPermutation[j]] = (*lambda)[j];
}
mapping[i] = indexOfCoords_[newLambda];
}
......@@ -91,7 +89,7 @@ namespace AMDiS {
: BoundaryCondition(type, rowFESpace_, NULL),
masterMatrix_(NULL)
{
if(rowFESpace->getMesh()->getDim() > 1) {
if (rowFESpace->getMesh()->getDim() > 1) {
periodicDOFMapping_ =
PeriodicDOFMapping::providePeriodicDOFMapping(rowFESpace_->getBasisFcts());
} else {
......@@ -111,10 +109,10 @@ namespace AMDiS {
Mesh *mesh = matrix->getRowFESpace()->getMesh();
associated_ = mesh->getPeriodicAssociations()[boundaryType];
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);
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
neighIndices_ = new DegreeOfFreedom[basFcts->getNumber()];
......@@ -131,7 +129,6 @@ namespace AMDiS {
int dim = rowFESpace->getMesh()->getDim();
if (dim > 1) {
int i, j;
DOFAdmin *admin = rowFESpace->getAdmin();
FixVec<int, WORLD> elFace(dim, NO_INIT);
......@@ -157,24 +154,25 @@ namespace AMDiS {
index = element->getVertexOfPosition(sideGeoIndex,
side,
vertex);
periodicDOFs[vertex] = (*associated_)[dofIndices[index]];
periodicDOFs[vertex] = (*associated)[dofIndices[index]];
}
Element *neigh = elInfo->getNeighbour(side);
basFcts->getLocalIndices(neigh, admin, neighIndices_);
int oppVertex = 0;
for(i = 0; i < dim + 1; i++) {
for (int i = 0; i < dim + 1; i++) {
// get vertex permutation
if(i == side) {
if (i == side) {
vertexPermutation[i] = 0;
} else {
DegreeOfFreedom periodicDOF =
periodicDOFs[element->getPositionOfVertex(side, i)];
for(j = 0; j < dim + 1; j++) {
if(neigh->getDOF(j, 0) == periodicDOF) break;
}
for (int j = 0; j < dim + 1; j++)
if (neigh->getDOF(j, 0) == periodicDOF)
break;
vertexPermutation[i] = j;
}
oppVertex += i - vertexPermutation[i];
......@@ -186,11 +184,9 @@ namespace AMDiS {
periodicDOFMapping_->getDOFPermutation(vertexPermutation);
// set associated dofs
for(i = 0; i < num; i++) {
if((*(basFcts->getCoords(i)))[side] == 0) {
(*associated_)[dofIndices[i]] = neighIndices_[dofPermutation[i]];
}
}
for (int i = 0; i < num; i++)
if ((*(basFcts->getCoords(i)))[side] == 0)
(*associated)[dofIndices[i]] = neighIndices_[dofPermutation[i]];
}
}
}
......@@ -201,8 +197,6 @@ namespace AMDiS {
{
FUNCNAME("PeriodicBC::exitMatrix()");
ERROR_EXIT("Not yet implemented, if you have an idea how to do it, than make it!!!\n");
TEST_EXIT(matrix)("no matrix\n");
if (matrix == masterMatrix_) {
......@@ -212,10 +206,12 @@ namespace AMDiS {
using namespace mtl;
associated_->print();
std::cout << "ASSOC = " << std::endl;
for (int i = 0; i < 10; i++)
std::cout << i << " = " << (*associated)[i] << std::endl;
// Compute reorder matrix (newRow and newCol yields transposed!!!)
matrix::traits::reorder<>::type R= matrix::reorder(*associated_);
matrix::traits::reorder<>::type R= matrix::reorder(*associated);
DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), B, D, E, TR;
A*= 0.5;
......@@ -243,12 +239,12 @@ namespace AMDiS {
DOFIterator<double> vecIt(vector, USED_DOFS);
Mesh *mesh = vector->getFESpace()->getMesh();
VertexVector *associated_ = mesh->getPeriodicAssociations()[boundaryType];
VertexVector *associated = mesh->getPeriodicAssociations()[boundaryType];
for (vecIt.reset(); !vecIt.end(); ++vecIt) {
dof = vecIt.getDOFIndex();
newDOF = (*associated_)[dof];
newDOF = (*associated)[dof];
if (dof < newDOF) {
entry = ((*vector)[dof] + (*vector)[newDOF]) * 0.5;
......
......@@ -34,7 +34,8 @@ namespace AMDiS {
class DimVecLess
{
public:
bool operator()(const DimVec<T> &v1, const DimVec<T> &v2) {
bool operator()(const DimVec<T> &v1, const DimVec<T> &v2)
{
int size = v1.getSize();
for (int i = 0; i < size; i++) {
if (v1[i] < v2[i])
......@@ -43,7 +44,7 @@ namespace AMDiS {
return false;
}
return false;
};
}
};
class PeriodicDOFMapping
......@@ -60,7 +61,9 @@ namespace AMDiS {
protected:
const BasisFunction *basFcts_;
std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> > dofPermutation_;
std::map<DimVec<double>, int, DimVecLess<double> > indexOfCoords_;
static std::vector<PeriodicDOFMapping*> mappings_;
......@@ -83,18 +86,18 @@ namespace AMDiS {
void initMatrix(DOFMatrix* matrix);
void fillBoundaryCondition(DOFMatrix *matrix,
ElInfo *elInfo,
void fillBoundaryCondition(DOFMatrix *matrix,
ElInfo *elInfo,
const DegreeOfFreedom *dofIndices,
const BoundaryType *localBound,
int nBasFcts);
const BoundaryType *localBound,
int nBasFcts);
void exitMatrix(DOFMatrix* matrix);
void exitVector(DOFVectorBase<double>* vector);
protected:
VertexVector *associated_;
VertexVector *associated;
PeriodicDOFMapping *periodicDOFMapping_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment