Commit f615ae6e authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Nothing really important

parent ada4194a
...@@ -38,7 +38,7 @@ namespace AMDiS { ...@@ -38,7 +38,7 @@ namespace AMDiS {
{ {
DELETE nDOF; DELETE nDOF;
for (int i = 0; i < grdTmpVec1.size(); i++) { for (int i = 0; i < static_cast<int>(grdTmpVec1.size()); i++) {
DELETE grdTmpVec1[i]; DELETE grdTmpVec1[i];
DELETE grdTmpVec2[i]; DELETE grdTmpVec2[i];
} }
......
...@@ -80,7 +80,6 @@ namespace AMDiS { ...@@ -80,7 +80,6 @@ namespace AMDiS {
if (localBCs.size() > 0) { if (localBCs.size() > 0) {
// get boundaries of all DOFs // get boundaries of all DOFs
const BoundaryType *localBound = basisFcts->getBound(elInfo, NULL); const BoundaryType *localBound = basisFcts->getBound(elInfo, NULL);
// get dof indices // get dof indices
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices); basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
...@@ -88,11 +87,11 @@ namespace AMDiS { ...@@ -88,11 +87,11 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
if (!(*it).second->isDirichlet()) { if (!(*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts); (*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts);
} }
} }
} }
// apply dirichlet boundary conditions // apply dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) { for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) { if ((*it).second) {
......
...@@ -2050,58 +2050,59 @@ namespace AMDiS { ...@@ -2050,58 +2050,59 @@ namespace AMDiS {
int MacroReader::basicDOFCheckFct(ElInfo* elinfo) int MacroReader::basicDOFCheckFct(ElInfo* elinfo)
{ {
FUNCNAME("MacroReader::basicDOFCheckFct"); FUNCNAME("MacroReader::basicDOFCheckFct()");
Mesh* mesh = Mesh::traversePtr; Mesh* mesh = Mesh::traversePtr;
Element* el = elinfo->getElement(); Element* el = elinfo->getElement();
const DOFAdmin& admin = mesh->getDOFAdmin(mesh->iadmin);
const DOFAdmin& adm = mesh->getDOFAdmin(mesh->iadmin); const Element *neig;
const Element *neig; const DegreeOfFreedom *dof;
const DegreeOfFreedom *dof;
int i, j, jdof, ndof, i0, j0, ov;
if (0 == mesh->dof_used.size()) if (0 == mesh->dof_used.size())
return 0; return 0;
if ((ndof = adm.getNumberOfDOFs(VERTEX))) { int ndof = admin.getNumberOfDOFs(VERTEX);
j0 = adm.getNumberOfPreDOFs(VERTEX); if (ndof) {
TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(VERTEX)) int j0 = admin.getNumberOfPreDOFs(VERTEX);
("adm.getNumberOfPreDOFs(VERTEX) %d + nDOF %d > mesh->nDOF %d\n", TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(VERTEX))
j0, ndof, mesh->getNumberOfDOFs(VERTEX)); ("admin.getNumberOfPreDOFs(VERTEX) %d + nDOF %d > mesh->nDOF %d\n",
i0 = mesh->getNode(VERTEX); j0, ndof, mesh->getNumberOfDOFs(VERTEX));
for (i = 0; i < mesh->getGeo(VERTEX); i++) int i0 = mesh->getNode(VERTEX);
{ for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
if ((dof = el->getDOF(i0+i)) == NULL) if ((dof = el->getDOF(i0+i)) == NULL) {
ERROR("no vertex dof %d on element %d\n", i, el->getIndex()); ERROR("no vertex dof %d on element %d\n", i, el->getIndex());
else } else {
for (j = 0; j < ndof; j++) for (int j = 0; j < ndof; j++) {
{ int jdof = dof[j0 + j];
jdof = dof[j0 + j]; TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size())) ("vertex dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
("vertex dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size()); mesh->dof_used[jdof]++;
mesh->dof_used[jdof]++;
}
} }
/* neighbour vertex dofs have been checked in check_fct() */ }
} }
/* neighbour vertex dofs have been checked in check_fct() */
}
if (mesh->getDim() > 1) { if (mesh->getDim() > 1) {
if ((ndof = adm.getNumberOfDOFs(EDGE))) { ndof = admin.getNumberOfDOFs(EDGE);
j0 = adm.getNumberOfPreDOFs(EDGE); if (ndof) {
int j0 = admin.getNumberOfPreDOFs(EDGE);
TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(EDGE)) TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(EDGE))
("adm.getNumberOfPreDOFs(EDGE) %d + nDOF %d > mesh->nDOF %d\n", ("admin.getNumberOfPreDOFs(EDGE) %d + nDOF %d > mesh->nDOF %d\n",
j0, ndof, mesh->getNumberOfDOFs(EDGE)); j0, ndof, mesh->getNumberOfDOFs(EDGE));
i0 = mesh->getNode(EDGE); int i0 = mesh->getNode(EDGE);
for (i = 0; i < mesh->getGeo(EDGE); i++) { for (int i = 0; i < mesh->getGeo(EDGE); i++) {
if ((dof = el->getDOF(i0 + i)) == NULL) { dof = el->getDOF(i0 + i);
if (dof == NULL) {
ERROR("no edge dof %d on element %d\n", i, el->getIndex()); ERROR("no edge dof %d on element %d\n", i, el->getIndex());
} else { } else {
for (j = 0; j < ndof; j++) { for (int j = 0; j < ndof; j++) {
jdof = dof[j0 + j]; int jdof = dof[j0 + j];
TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size())) TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
("edge dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size()); ("edge dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
mesh->dof_used[jdof]++; mesh->dof_used[jdof]++;
...@@ -2109,9 +2110,10 @@ namespace AMDiS { ...@@ -2109,9 +2110,10 @@ namespace AMDiS {
} }
if (el->getFirstChild() == NULL) { if (el->getFirstChild() == NULL) {
if(mesh->getDim() == 2) { if (mesh->getDim() == 2) {
if ((neig = elinfo->getNeighbour(i))) { neig = elinfo->getNeighbour(i);
ov = elinfo->getOppVertex(i); if (neig) {
int ov = elinfo->getOppVertex(i);
TEST(neig->getDOF(i0 + ov) == dof) TEST(neig->getDOF(i0 + ov) == dof)
("el %d edge %d dof %8X: wrong dof %8X in neighbour %d edge %d\n", ("el %d edge %d dof %8X: wrong dof %8X in neighbour %d edge %d\n",
...@@ -2119,15 +2121,14 @@ namespace AMDiS { ...@@ -2119,15 +2121,14 @@ namespace AMDiS {
neig->getIndex(), ov); neig->getIndex(), ov);
} }
} else { // dim == 3 } else { // dim == 3
int in, k, found; for (int in = 0; in < mesh->getGeo(NEIGH); in++) {
for (in = 0; in < mesh->getGeo(NEIGH); in++) {
if ((in != el->getVertexOfEdge(i,0)) && if ((in != el->getVertexOfEdge(i,0)) &&
(in != el->getVertexOfEdge(i,1)) && (in != el->getVertexOfEdge(i,1)) &&
(neig = elinfo->getNeighbour(in))) { (neig = elinfo->getNeighbour(in))) {
found = 0; int found = 0;
for (k = 0; k < mesh->getGeo(EDGE); k++) for (int k = 0; k < mesh->getGeo(EDGE); k++) {
if (neig->getDOF(i0 + k) == dof) found++; if (neig->getDOF(i0 + k) == dof) found++;
}
TEST(found==1)("el %d edge %d dof found=%d in neighbour %d\n", TEST(found==1)("el %d edge %d dof found=%d in neighbour %d\n",
el->getIndex(), i, found, neig->getIndex()); el->getIndex(), i, found, neig->getIndex());
} }
...@@ -2138,17 +2139,18 @@ namespace AMDiS { ...@@ -2138,17 +2139,18 @@ namespace AMDiS {
} }
} }
if(mesh->getDim()==3) { if (mesh->getDim() == 3) {
if ((ndof = adm.getNumberOfDOFs(FACE))) { ndof = admin.getNumberOfDOFs(FACE);
j0 = adm.getNumberOfPreDOFs(FACE); if (ndof) {
int j0 = admin.getNumberOfPreDOFs(FACE);
TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(FACE)) TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(FACE))
("admin->n0_dof[FACE] %d + nDOF %d > mesh->nDOF %d\n", ("admin->n0_dof[FACE] %d + nDOF %d > mesh->nDOF %d\n",
j0, ndof, mesh->getNumberOfDOFs(FACE)); j0, ndof, mesh->getNumberOfDOFs(FACE));
i0 = mesh->getNode(FACE); int i0 = mesh->getNode(FACE);
for (i = 0; i < mesh->getGeo(FACE); i++) { for (int i = 0; i < mesh->getGeo(FACE); i++) {
TEST(dof = el->getDOF(i0 + i))("no face dof %d ???\n", i); TEST(dof = el->getDOF(i0 + i))("no face dof %d ???\n", i);
for (j = 0; j < ndof; j++) { for (int j = 0; j < ndof; j++) {
jdof = dof[j0 + j]; int jdof = dof[j0 + j];
TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size())) TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
("face dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size()); ("face dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
mesh->dof_used[jdof]++; mesh->dof_used[jdof]++;
...@@ -2156,7 +2158,8 @@ namespace AMDiS { ...@@ -2156,7 +2158,8 @@ namespace AMDiS {
if (el->getChild(0) == NULL) { if (el->getChild(0) == NULL) {
if ((neig = elinfo->getNeighbour(i))) { if ((neig = elinfo->getNeighbour(i))) {
ov = elinfo->getOppVertex(i); int ov = elinfo->getOppVertex(i);
TEST(neig->getDOF(i0 + ov) == dof) TEST(neig->getDOF(i0 + ov) == dof)
("el %d face %d dof %8X: wrong dof %8X in neighbour %d face %d\n", ("el %d face %d dof %8X: wrong dof %8X in neighbour %d face %d\n",
el->getIndex(), i, dof, neig->getDOF(i0 + ov), neig->getIndex(), el->getIndex(), i, dof, neig->getDOF(i0 + ov), neig->getIndex(),
...@@ -2167,20 +2170,21 @@ namespace AMDiS { ...@@ -2167,20 +2170,21 @@ namespace AMDiS {
} }
} }
if ((ndof = adm.getNumberOfDOFs(CENTER))) { ndof = admin.getNumberOfDOFs(CENTER);
i0 = mesh->getNode(CENTER); if (ndof) {
TEST(dof = el->getDOF(i0))("no center dof???\n"); int i0 = mesh->getNode(CENTER);
j0 = adm.getNumberOfPreDOFs(CENTER); TEST(dof = el->getDOF(i0))("no center dof???\n");
TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(CENTER)) int j0 = admin.getNumberOfPreDOFs(CENTER);
("adm.getNumberOfPreDOFs(CENTER) %d + nDOF %d > mesh->nDOF %d\n", TEST_EXIT(j0 + ndof <= mesh->getNumberOfDOFs(CENTER))
j0, ndof, mesh->getNumberOfDOFs(CENTER)); ("admin.getNumberOfPreDOFs(CENTER) %d + nDOF %d > mesh->nDOF %d\n",
for (j = 0; j < ndof; j++) { j0, ndof, mesh->getNumberOfDOFs(CENTER));
jdof = dof[j0 + j]; for (int j = 0; j < ndof; j++) {
TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size())) int jdof = dof[j0 + j];
("center dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size()); TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
mesh->dof_used[jdof]++; ("center dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
} mesh->dof_used[jdof]++;
} }
}
return 0; return 0;
} }
......
...@@ -1069,21 +1069,21 @@ namespace AMDiS { ...@@ -1069,21 +1069,21 @@ namespace AMDiS {
bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) { bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
::std::vector<DegreeOfFreedom> associatedToDOF1; ::std::vector<DegreeOfFreedom> associatedToDOF1;
int i, size;
::std::map<BoundaryType, VertexVector*>::iterator it; ::std::map<BoundaryType, VertexVector*>::iterator it;
::std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end(); ::std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
DegreeOfFreedom dof, assDOF; DegreeOfFreedom dof, assDOF;
associatedToDOF1.push_back(dof1); associatedToDOF1.push_back(dof1);
for(it = periodicAssociations.begin(); it != end; ++it) { for (it = periodicAssociations.begin(); it != end; ++it) {
size = static_cast<int>(associatedToDOF1.size()); int size = static_cast<int>(associatedToDOF1.size());
for(i = 0; i < size; i++) { for (int i = 0; i < size; i++) {
dof = associatedToDOF1[i]; dof = associatedToDOF1[i];
assDOF = (*(it->second))[dof]; assDOF = (*(it->second))[dof];
if(assDOF == dof2) { if (assDOF == dof2) {
return true; return true;
} else { } else {
if(assDOF != dof) associatedToDOF1.push_back(assDOF); if (assDOF != dof)
associatedToDOF1.push_back(assDOF);
} }
} }
} }
......
...@@ -617,18 +617,6 @@ namespace AMDiS { ...@@ -617,18 +617,6 @@ namespace AMDiS {
*/ */
inline bool isInitialized() { return initialized; }; inline bool isInitialized() { return initialized; };
// inline void addPeriodicBC(BoundaryType type) {
// periodicBoundaryTypes.insert(type);
// };
// inline bool isPeriodicBC(BoundaryType type) {
// return (periodicBoundaryTypes.find(type) != periodicBoundaryTypes.end());
// };
// inline ::std::map<BoundaryType, PeriodicBC*>& getPeriodicBCMap() {
// return periodicBoundaryConditions;
// };
inline ::std::map<BoundaryType, VertexVector*>& getPeriodicAssociations() { inline ::std::map<BoundaryType, VertexVector*>& getPeriodicAssociations() {
return periodicAssociations; return periodicAssociations;
}; };
...@@ -817,11 +805,6 @@ namespace AMDiS { ...@@ -817,11 +805,6 @@ namespace AMDiS {
*/ */
bool preserveCoarseDOFs; bool preserveCoarseDOFs;
// /** \brief
// * List of all Meshes. Can be accessed via Mesh::begin() and Mesh::end()
// */
// static ::std::list<Mesh*> meshes;
/** \brief /** \brief
* Number of all DOFs on a single element * Number of all DOFs on a single element
*/ */
......
...@@ -109,18 +109,14 @@ namespace AMDiS { ...@@ -109,18 +109,14 @@ namespace AMDiS {
if (!masterMatrix_) { if (!masterMatrix_) {
masterMatrix_ = matrix; masterMatrix_ = matrix;
Mesh *mesh = matrix->getRowFESpace()->getMesh(); Mesh *mesh = matrix->getRowFESpace()->getMesh();
associated_ = mesh->getPeriodicAssociations()[boundaryType]; associated_ = mesh->getPeriodicAssociations()[boundaryType];
TEST_EXIT_DBG(associated_)("no associations for periodic boundary condition %d\n", TEST_EXIT(associated_)
boundaryType); ("no associations for periodic boundary condition %d\n", boundaryType);
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int num = basFcts->getNumber();
neighIndices_ = GET_MEMORY(DegreeOfFreedom, num); neighIndices_ = GET_MEMORY(DegreeOfFreedom,
rowFESpace->getBasisFcts()->getNumber());
} }
} }
...@@ -131,30 +127,22 @@ namespace AMDiS { ...@@ -131,30 +127,22 @@ namespace AMDiS {
int nBasFcts) int nBasFcts)
{ {
if (matrix == masterMatrix_) { if (matrix == masterMatrix_) {
int dim = rowFESpace->getMesh()->getDim(); int dim = rowFESpace->getMesh()->getDim();
if (dim > 1) { if (dim > 1) {
DOFAdmin *admin = rowFESpace->getAdmin(); DOFAdmin *admin = rowFESpace->getAdmin();
FixVec<int, WORLD> elFace(dim, NO_INIT); FixVec<int, WORLD> elFace(dim, NO_INIT);
FixVec<int, WORLD> neighFace(dim, NO_INIT); FixVec<int, WORLD> neighFace(dim, NO_INIT);
DimVec<int> vertexPermutation(dim, NO_INIT); DimVec<int> vertexPermutation(dim, NO_INIT);
const BasisFunction *basFcts = rowFESpace->getBasisFcts(); const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int num = basFcts->getNumber(); int num = basFcts->getNumber();
Element *element = elInfo->getElement(); Element *element = elInfo->getElement();
DimVec<DegreeOfFreedom> periodicDOFs(dim-1, NO_INIT); DimVec<DegreeOfFreedom> periodicDOFs(dim-1, NO_INIT);
int vertex, index, side; int vertex, index, side;
GeoIndex sideGeoIndex = INDEX_OF_DIM(dim-1, dim); GeoIndex sideGeoIndex = INDEX_OF_DIM(dim-1, dim);
for (side = 0; side < dim + 1; side++) { for (side = 0; side < dim + 1; side++) {
if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) { if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
for (vertex = 0; vertex < dim; vertex++) { for (vertex = 0; vertex < dim; vertex++) {
index = element->getVertexOfPosition(sideGeoIndex, index = element->getVertexOfPosition(sideGeoIndex,
side, side,
...@@ -163,10 +151,9 @@ namespace AMDiS { ...@@ -163,10 +151,9 @@ namespace AMDiS {
} }
Element *neigh = elInfo->getNeighbour(side); Element *neigh = elInfo->getNeighbour(side);
basFcts->getLocalIndices(neigh, admin, neighIndices_); basFcts->getLocalIndices(neigh, admin, neighIndices_);
int oppVertex = 0; int oppVertex = 0;
for (int i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
// get vertex permutation // get vertex permutation
if (i == side) { if (i == side) {
...@@ -184,12 +171,13 @@ namespace AMDiS { ...@@ -184,12 +171,13 @@ namespace AMDiS {
} }
oppVertex += i - vertexPermutation[i]; oppVertex += i - vertexPermutation[i];
} }
vertexPermutation[side] = oppVertex; vertexPermutation[side] = oppVertex;
// get DOF permutation // get DOF permutation
const DegreeOfFreedom *dofPermutation = const DegreeOfFreedom *dofPermutation =
periodicDOFMapping_->getDOFPermutation(vertexPermutation); periodicDOFMapping_->getDOFPermutation(vertexPermutation);
// set associated dofs // set associated dofs
for (int i = 0; i < num; i++) { for (int i = 0; i < num; i++) {
if ((*(basFcts->getCoords(i)))[side] == 0) { if ((*(basFcts->getCoords(i)))[side] == 0) {
......
...@@ -603,15 +603,15 @@ namespace AMDiS { ...@@ -603,15 +603,15 @@ namespace AMDiS {
Mesh::FILL_DET | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH); Mesh::FILL_NEIGH);
// for all elements ... // for all elements ...
while (elInfo) { while (elInfo) {
if (systemMatrix_->getBoundaryManager()) if (systemMatrix_->getBoundaryManager())
systemMatrix_->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix_); systemMatrix_->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix_);
if (rhs_->getBoundaryManager()) if (rhs_->getBoundaryManager())
rhs_->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs_); rhs_->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs_);
if (solution_->getBoundaryManager()) if (solution_->getBoundaryManager())
solution_->getBoundaryManager()->fillBoundaryConditions(elInfo, solution_); solution_->getBoundaryManager()->fillBoundaryConditions(elInfo, solution_);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
......
...@@ -42,6 +42,8 @@ namespace AMDiS { ...@@ -42,6 +42,8 @@ namespace AMDiS {
fineSolutions.clear(); fineSolutions.clear();
coarseSolutions.clear(); coarseSolutions.clear();
} }
virtual ~ParaRealProblemBase() {};
void storeSolution(T *vec) void storeSolution(T *vec)
{ {
......
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