diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index eb5299a3201453a429256ec0f29c120a035c8087..b4e6f8f5a555b715fe261a708f6f10ef8f6caafb 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -116,6 +116,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+
// Remove all interpolation marks and, instead, set to each
// interpolation point its continous index starting from 0.
int i = 0;
@@ -126,7 +127,8 @@ namespace AMDiS {
}
// Traverse elements to create interpolation values.
- elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
+ counter = 0;
+ elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_ | Mesh::FILL_COORDS);
while (elInfo) {
if (!writeElem_ || writeElem_(elInfo))
addInterpData(elInfo);
@@ -268,102 +270,73 @@ namespace AMDiS {
{
FUNCNAME("DataCollector::addValueData()");
- const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
-
- /* vertex dofs */
- int dof_offset = localAdmin_->getNumberOfPreDOFs(VERTEX);
-
+ const BasisFunction *basisFcts = feSpace_->getBasisFcts();
+ const DegreeOfFreedom *localDOFs = basisFcts->getLocalIndices(elInfo->getElement(), localAdmin_, NULL);
+ const int nBasisFcts = basisFcts->getNumber();
+
+
for (int i = 0; i < mesh_->getGeo(VERTEX); i++) {
- (*interpPointInd_)[dof[i][dof_offset]] = -2; // mark as vertex
+ (*interpPointInd_)[localDOFs[i]] = -2; // mark as vertex
// get coords of this vertex
WorldVector<double> vertexCoords = elInfo->getCoord(i);
// search for coords at this dof
::std::list<WorldVector<double> >::iterator it =
- find((*dofCoords_)[dof[i][dof_offset]].begin(),
- (*dofCoords_)[dof[i][dof_offset]].end(),
+ find((*dofCoords_)[localDOFs[i]].begin(),
+ (*dofCoords_)[localDOFs[i]].end(),
vertexCoords);
// coords not yet in list?
- if (it == (*dofCoords_)[dof[i][dof_offset]].end()) {
+ if (it == (*dofCoords_)[localDOFs[i]].end()) {
// add new coords to list
- (*dofCoords_)[dof[i][dof_offset]].push_back(vertexCoords);
+ (*dofCoords_)[localDOFs[i]].push_back(vertexCoords);
}
}
+
- int nInterpPoints = 0;
- const BasisFunction *basisFcts = feSpace_->getBasisFcts();
-
- for (int i = 1; i <= dim_; i++) {
- int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
- int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
- dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-
- for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
- int node = node_offset + j;
-
- for (int k = 0; k < num_dofs; k++) {
- int dof_index = dof_offset + k;
-
-
- WorldVector<double> interpolCoords;
- elInfo->coordToWorld((*basisFcts->getCoords(mesh_->getGeo(VERTEX) + nInterpPoints)),
- &interpolCoords);
-
- nInterpPoints++;
-
- if ((*interpPointInd_)[dof[node][dof_index]] == -1) {
- // mark as interpolation point
- (*interpPointInd_)[dof[node][dof_index]] = -3;
-
- // search for interpolation point coordinates, and insert them to the
- // dof-entry, if not contained in the list
- ::std::list<WorldVector<double> >::iterator it =
- find((*interpPointCoords_)[dof[node][dof_index]].begin(),
- (*interpPointCoords_)[dof[node][dof_index]].end(),
+ for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts; i++) {
+ WorldVector<double> interpolCoords;
+ elInfo->coordToWorld(*basisFcts->getCoords(i), &interpolCoords);
+
+ if ((*interpPointInd_)[localDOFs[i]] == -1) {
+ // mark as interpolation point
+ (*interpPointInd_)[localDOFs[i]] = -3;
+
+ // search for interpolation point coordinates, and insert them to the
+ // dof-entry, if not contained in the list
+ ::std::list<WorldVector<double> >::iterator it =
+ find((*interpPointCoords_)[localDOFs[i]].begin(),
+ (*interpPointCoords_)[localDOFs[i]].end(),
interpolCoords);
- if (it == (*interpPointCoords_)[dof[node][dof_index]].end()) {
- (*interpPointCoords_)[dof[node][dof_index]].push_back(interpolCoords);
- }
-
-
- nInterpPoints_++;
- }
+
+ if (it == (*interpPointCoords_)[localDOFs[i]].end()) {
+ (*interpPointCoords_)[localDOFs[i]].push_back(interpolCoords);
+ nInterpPoints_++;
}
- }
+ }
}
-
+
return(0);
}
int DataCollector::addInterpData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addInterpData()");
-
- const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
-
+
::std::vector<int> elemInterpPoints;
elemInterpPoints.clear();
- for (int i = 1; i <= dim_; i++) {
- int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
- int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
- int dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-
- for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
- int node = node_offset + j;
+ const BasisFunction *basisFcts = feSpace_->getBasisFcts();
+ const DegreeOfFreedom *localDOFs = basisFcts->getLocalIndices(elInfo->getElement(), localAdmin_, NULL);
+ const int nBasisFcts = basisFcts->getNumber();
- for (int k = 0; k < num_dofs; k++) {
- int dof_index = dof_offset + k;
-
- elemInterpPoints.push_back((*interpPointInd_)[dof[node][dof_index]]);
- }
- }
+ for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts; i++) {
+ elemInterpPoints.push_back((*interpPointInd_)[localDOFs[i]]);
}
interpPoints_.push_back(elemInterpPoints);
-
+
return(0);
}
diff --git a/AMDiS/src/DataCollector.h b/AMDiS/src/DataCollector.h
index a7a23fb19f105c6b7fb8f1c7961491dfc38381ee..0ee307854714abf6be73f5930245a3d23c5a2e5e 100644
--- a/AMDiS/src/DataCollector.h
+++ b/AMDiS/src/DataCollector.h
@@ -294,6 +294,8 @@ namespace AMDiS {
* Pointer to a function which decides whether an element is considered.
*/
bool (*writeElem_)(ElInfo*);
+
+ int counter;
};
}
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index 39ba45978e9644cd7d722d0166d0d0386795bc80..6a367857f951fc2bbc47af459d3fd9637a146d07 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -12,8 +12,11 @@ namespace AMDiS {
{
ElementRegion_ED* red_;
- if (!elementData) return -1;
- red_=dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
+ if (!elementData)
+ return -1;
+
+ red_ = dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
+
if (red_)
return red_->getRegion();
@@ -22,8 +25,9 @@ namespace AMDiS {
void Element::setDOFPtrs()
{
- FUNCNAME("Element::setDOFPtrs");
+ FUNCNAME("Element::setDOFPtrs()");
TEST_EXIT(mesh)("no mesh!\n");
+
dof = mesh->createDOFPtrs();
}
@@ -33,27 +37,30 @@ namespace AMDiS {
index = mesh ? mesh->getNextElementIndex() : -1;
- child[0] = NULL;
- child[1] = NULL;
+ child[0] = NULL;
+ child[1] = NULL;
- dof = mesh ? mesh->createDOFPtrs() : NULL;
newCoord = NULL;
-
elementData = NULL;
+
+ if (mesh) {
+ setDOFPtrs();
+ } else {
+ mesh = NULL;
+ }
}
// call destructor through Mesh::freeElement !!!
Element::~Element()
{
- if(child[0]) DELETE child[0];
- if(child[1]) DELETE child[1];
-
- //if(elementData) DELETE elementData;
+ if (child[0])
+ DELETE child[0];
+ if (child[1])
+ DELETE child[1];
- if (newCoord)
- {
- DELETE newCoord;
- }
+ if (newCoord) {
+ DELETE newCoord;
+ }
}
/****************************************************************************/
@@ -65,7 +72,7 @@ namespace AMDiS {
/* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs */
#define CHANGE_DOFS_1(el) \
- ldof = el->dof[n0+i] + nd0; \
+ ldof = el->dof[n0 + i] + nd0; \
for (j = 0; j < nd; j++) { \
if ((k = ldof[j]) >= 0) { \
/* do it only once! (dofs are visited more than once) */ \
@@ -99,7 +106,7 @@ namespace AMDiS {
}
}
- if(mesh->getDim() > 1) {
+ if (mesh->getDim() > 1) {
if ((nd = admin->getNumberOfDOFs(EDGE))) {
nd0 = admin->getNumberOfPreDOFs(EDGE);
n0 = admin->getMesh()->getNode(EDGE);
@@ -109,7 +116,7 @@ namespace AMDiS {
}
}
- if (3==mesh->getDim()) {
+ if (mesh->getDim() == 3) {
if ((nd = admin->getNumberOfDOFs(FACE))) {
nd0 = admin->getNumberOfPreDOFs(FACE);
n0 = admin->getMesh()->getNode(FACE);
@@ -145,7 +152,7 @@ namespace AMDiS {
}
}
- if(mesh->getDim() > 1) {
+ if (mesh->getDim() > 1) {
if ((nd = admin->getNumberOfDOFs(EDGE))) {
nd0 = admin->getNumberOfPreDOFs(EDGE);
n0 = admin->getMesh()->getNode(EDGE);
@@ -155,7 +162,7 @@ namespace AMDiS {
}
}
- if (3==mesh->getDim()) {
+ if (mesh->getDim() == 3) {
if ((nd = admin->getNumberOfDOFs(FACE))) {
nd0 = admin->getNumberOfPreDOFs(FACE);
n0 = admin->getMesh()->getNode(FACE);
@@ -183,30 +190,30 @@ namespace AMDiS {
int Element::oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const
{
- int i, j, nv = 0, ov = 0;
+ int nv = 0, ov = 0;
int vertices = mesh->getGeo(VERTEX);
int dim = mesh->getDim();
- for (i = 0; i < vertices; i++)
- {
- if (nv < i-1) return(-1);
-
- for (j = 0; j < dim; j++)
- {
- if (dof[i] == pdof[j])
- {
- /****************************************************************************/
- /* i is a common vertex */
- /****************************************************************************/
- ov += i;
- nv++;
- break;
- }
- }
-
+ for (int i = 0; i < vertices; i++) {
+ if (nv < i-1)
+ return(-1);
+
+ for (int j = 0; j < dim; j++) {
+ if (dof[i] == pdof[j]) {
+ /****************************************************************************/
+ /* i is a common vertex */
+ /****************************************************************************/
+ ov += i;
+ nv++;
+ break;
+ }
}
- if (nv != mesh->getDim()) return(-1);
+
+ }
+
+ if (nv != mesh->getDim())
+ return(-1);
/****************************************************************************/
/* the opposite vertex is 3(6) - (sum of indices of common vertices) in */
/* 2d(3d) */
@@ -243,7 +250,7 @@ namespace AMDiS {
newCoord=NULL;
};
}
-
+
void Element::serialize(::std::ostream &out)
{
// write children
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index 2df1b7ff8c605424859313d9c0f437dfe7292cbd..7696a70ecfaa4a818026ca78503b9b9f939aec31 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -305,9 +305,9 @@ namespace AMDiS {
/** \brief
* Sets the pointer to the DOFs of the i-th node of Element
*/
- DegreeOfFreedom* setDOF(int i, DegreeOfFreedom* p) {
- dof[i] = p;
- return dof[i];
+ DegreeOfFreedom* setDOF(int pos, DegreeOfFreedom* p) {
+ dof[pos] = p;
+ return dof[pos];
};
/** \brief
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index dcf156e3523fa50d47c9abbdd9c27d38a16d8c87..886d615e56f46a9d011319e4c222d37e23cd9250 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -360,7 +360,9 @@ namespace AMDiS {
}
}
- Lagrange::Phi::~Phi() { DELETE [] vertices; }
+ Lagrange::Phi::~Phi() {
+ DELETE [] vertices;
+ }
Lagrange::GrdPhi::GrdPhi(Lagrange* owner_,
@@ -1003,16 +1005,17 @@ namespace AMDiS {
DegreeOfFreedom* result;
- if(indices) {
+ if (indices) {
result = indices;
} else {
- if(localVec) FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
+ if (localVec)
+ FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
localVecSize = nBasFcts;
result = localVec;
}
- for(pos=0, j=0; pos <= dim; pos++) {
+ for (pos = 0, j = 0; pos <= dim; pos++) {
posIndex = INDEX_OF_DIM(pos, dim);
n0 = admin->getNumberOfPreDOFs(posIndex);
node0 = admin->getMesh()->getNode(posIndex);
@@ -1032,33 +1035,6 @@ namespace AMDiS {
return result;
}
- // const double* Lagrange::getVec(const Element* el, const DOFVector<double> * dv,
- // double * d) const
- // {
- // static double* localVec = NULL;
- // const DOFAdmin* admin = dv->getFESpace()->getAdmin();
-
- // int i;
-
- // double* result;
-
- // if(d) {
- // result = d;
- // } else {
- // if(localVec) FREE_MEMORY(localVec, double, nBasFcts);
- // localVec = GET_MEMORY(double, nBasFcts);
- // result = localVec;
- // }
-
- // const DegreeOfFreedom *localIndices = getLocalIndices(el, admin, NULL);
-
- // for(i = 0; i < nBasFcts; i++) {
- // result[i] = (*dv)[localIndices[i]];
- // }
-
- // return result;
- // }
-
void Lagrange::l2ScpFctBas(Quadrature *q,
AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
DOFVector<WorldVector<double> >* fh)
@@ -1234,12 +1210,13 @@ namespace AMDiS {
RCNeighbourList* list,
int n, BasisFunction* basFct)
{
- FUNCNAME("Lagrange::refineInter2_2d");
+ FUNCNAME("Lagrange::refineInter2_2d()");
- if (n < 1) return;
+ if (n < 1)
+ return;
- Element *el;
- int node, n0;
+ Element *el;
+ int node, n0;
DegreeOfFreedom cdof;
const DegreeOfFreedom *pdof;
const DOFAdmin *admin = drv->getFESpace()->getAdmin();
@@ -1267,7 +1244,7 @@ namespace AMDiS {
cdof = el->getChild(0)->getDOF(node, n0);
(*drv)[cdof] =
- 0.375*(*drv)[pdof[0]] - 0.125*(*drv)[pdof[1]] + 0.75*(*drv)[pdof[5]];
+ 0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
/****************************************************************************/
/* node in the common edge of child[0] and child[1] */
@@ -1275,8 +1252,8 @@ namespace AMDiS {
cdof = el->getChild(0)->getDOF(node+1, n0);
(*drv)[cdof] =
- -0.125*((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25*(*drv)[pdof[5]]
- + 0.5*((*drv)[pdof[3]] + (*drv)[pdof[4]]);
+ -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
+ + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
/****************************************************************************/
/* midpoint of edge on child[1] at the refinement edge */
@@ -1284,22 +1261,20 @@ namespace AMDiS {
cdof = el->getChild(1)->getDOF(node+1, n0);
(*drv)[cdof] =
- -0.125*(*drv)[pdof[0]] + 0.375*(*drv)[pdof[1]] + 0.75*(*drv)[pdof[5]];
+ -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
- if (n > 1)
- {
- /****************************************************************************/
- /* adjust the value at the midpoint of the common edge of neigh's children */
- /****************************************************************************/
- el = list->getElement(1);
- pdof = basFct->getLocalIndices(el, admin, NULL);
-
- cdof = el->getChild(0)->getDOF(node+1, n0);
- (*drv)[cdof] =
- -0.125*((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25*(*drv)[pdof[5]]
- + 0.5*((*drv)[pdof[3]] + (*drv)[pdof[4]]);
- }
- return;
+ if (n > 1) {
+ /****************************************************************************/
+ /* adjust the value at the midpoint of the common edge of neigh's children */
+ /****************************************************************************/
+ el = list->getElement(1);
+ pdof = basFct->getLocalIndices(el, admin, NULL);
+
+ cdof = el->getChild(0)->getDOF(node+1, n0);
+ (*drv)[cdof] =
+ -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
+ + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
+ }
}
void Lagrange::refineInter2_3d(DOFIndexed<double> *drv,
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index c5d6183ed1fe4a36f91c91278e721635403361ed..51a38f009b0f5c5807d77540d34935a7059a6746 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -152,20 +152,6 @@ namespace AMDiS {
const DOFAdmin*,
DegreeOfFreedom*) const;
- /** \brief
- * Implements BasisFunction::getVec().
- */
- // const double *getVec(const Element* el,
- // const DOFVector<double> * dv,
- // double * d) const;
-
- /** \brief
- * Implements BasisFunction::getVec
- */
- // const WorldVector<double> *getVec(const Element* el,
- // const DOFVector<WorldVector<double> > * dv,
- // WorldVector<double> *d) const;
-
/** \brief
* Implements BasisFunction::l2ScpFctBas
*/
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 65b2c182f089b5da302e96d25e2223b76386aadb..a513b1278e01a746be1f707f88216ed4f2a2e332 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -133,16 +133,6 @@ namespace AMDiS {
element2->setElementData(ldp2);
}
- // DimVec<DegreeOfFreedom> periodicDOFsEl1(dim-1, NO_INIT);
- // DimVec<DegreeOfFreedom> periodicDOFsEl2(dim-1, NO_INIT);
-
- // for(j = 0; j < dim; j++) {
- // periodicDOFsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] =
- // melVertex[el2][vertexMapEl1[verticesEl1[j]]];
- // periodicDOFsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
- // melVertex[el1][vertexMapEl2[verticesEl2[j]]];
- // }
-
ldp1->addPeriodicInfo(mode,
boundaryType,
sideEl1,
@@ -228,8 +218,9 @@ namespace AMDiS {
for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
for (int k = 0; k < mesh->getGeo(VERTEX); k++) {
(*(mel + i))->setCoord(k, coords[melVertex[i][k]]);
+
const_cast<Element*>((*(mel+i))->getElement())->
- setDOF(k,dof[melVertex[i][k]]);
+ setDOF(k, dof[melVertex[i][k]]);
}
}
@@ -963,28 +954,18 @@ namespace AMDiS {
getElement()->
getDOF((k+l+1)%(dim+1)));
- for (j = i+1; j < mesh->getNumberOfLeaves(); j++)
- {
- //MacroElement *mel0 = mesh->getMacroElement(i);
- //Element *el = const_cast<Element*>(mel0->getElement());
- //m = el->oppVertex(dof);
- if ((m = mesh->getMacroElement(j)->getElement()->oppVertex(dof)) != -1)
- {
- mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
- mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
- mesh->getMacroElement(i)->setOppVertex(k, m);
- mesh->getMacroElement(j)->setOppVertex(m, k);
- break;
- }
+ for (j = i + 1; j < mesh->getNumberOfLeaves(); j++) {
+ if ((m = mesh->getMacroElement(j)->getElement()->oppVertex(dof)) != -1) {
+ mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
+ mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
+ mesh->getMacroElement(i)->setOppVertex(k, m);
+ mesh->getMacroElement(j)->setOppVertex(m, k);
+ break;
}
+ }
+
TEST_EXIT(j < mesh->getNumberOfLeaves())
("could not find neighbour %d of element %d\n", k, i);
- // if(j == mesh->getNumberOfLeaves()) {
- // ::std::cout << "neighbour " << k << " not found" << ::std::endl;
- // mesh->getMacroElement(i)->setBoundary(k, 1111);
- // mesh->getMacroElement(i)->setNeighbour(k, NULL);
- // mesh->getMacroElement(i)->setOppVertex(k, -1);
- // }
}
}
}
@@ -1002,13 +983,14 @@ namespace AMDiS {
void MacroReader::boundaryDOFs(Mesh *mesh)
{
- FUNCNAME("Mesh::boundaryDOFs");
- int i, lnode = mesh->getNode(EDGE);
- int k, lne = mesh->getNumberOfLeaves();
- int max_n_neigh = 0, n_neigh, ov;
+ FUNCNAME("Mesh::boundaryDOFs()");
+
+ int lnode = mesh->getNode(EDGE);
+ int k, lne = mesh->getNumberOfLeaves();
+ int max_n_neigh = 0, n_neigh, ov;
::std::deque<MacroElement*>::iterator mel;
const MacroElement* neigh;
- DegreeOfFreedom *dof;
+ DegreeOfFreedom *dof;
mesh->setNumberOfEdges(0);
mesh->setNumberOfFaces(0);
@@ -1017,52 +999,52 @@ namespace AMDiS {
switch(dim) {
case 2:
- for (mel = mesh->firstMacroElement(); mel!=mesh->endOfMacroElements(); mel++) {
-
+ for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
// check for periodic boundary
Element *el = const_cast<Element*>((*mel)->getElement());
ElementData *ed = el->getElementData(PERIODIC);
DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
- if(ed) {
+ if (ed) {
::std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos =
dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
::std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
- for(it = periodicInfos.begin(); it != end; ++it) {
- if(it->type != 0) {
+ for (it = periodicInfos.begin(); it != end; ++it) {
+ if (it->type != 0) {
periodic[it->elementSide] = true;
}
}
}
- for (i = 0; i < mesh->getGeo(NEIGH); i++) {
+ for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
if (!(*mel)->getNeighbour(i) ||
- ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex()))
- {
- mesh->incrementNumberOfEdges(1);;
- if (mesh->getNumberOfDOFs(EDGE)) {
- dof = el->setDOF(lnode+i, mesh->getDOF(EDGE));
-
- if ((*mel)->getNeighbour(i)) {
- Element *neigh =
- const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
- if (periodic[i]) {
- neigh->setDOF(lnode+(*mel)->getOppVertex(i), mesh->getDOF(EDGE));
- } else {
- neigh->setDOF(lnode+(*mel)->getOppVertex(i), dof);
- }
+ ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {
+
+ mesh->incrementNumberOfEdges(1);
+
+ if (mesh->getNumberOfDOFs(EDGE)) {
+ dof = el->setDOF(lnode + i, mesh->getDOF(EDGE));
+
+ if ((*mel)->getNeighbour(i)) {
+ Element *neigh = const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
+
+ if (periodic[i]) {
+ neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDOF(EDGE));
+ } else {
+ neigh->setDOF(lnode + (*mel)->getOppVertex(i), dof);
}
}
- }
+ }
+ }
}
}
break;
case 3:
lnode = mesh->getNode(FACE);
mel = mesh->firstMacroElement();
- for (i = 0; i < lne; i++) {
+ for (int i = 0; i < lne; i++) {
// check for periodic boundary
Element *el = const_cast<Element*>((*(mel+i))->getElement());
@@ -1152,9 +1134,9 @@ namespace AMDiS {
default: ERROR_EXIT("invalid dim\n");
}
- if (3==dim)
+ if (3 == dim) {
mesh->setMaxEdgeNeigh(::std::max(8, 2*max_n_neigh));
- else {
+ } else {
mesh->setMaxEdgeNeigh(dim-1);
}
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 1743d4e91fd7705765db6e15233ab487800e9de5..05a703cbd457b9f841d54bf87747f05c119e5aa1 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -363,161 +363,79 @@ namespace AMDiS {
void Mesh::dofCompress()
{
- FUNCNAME("Mesh::dofCompress");
- int size;
- Flag fill_flag;
-
- // int i;
- // for(i = 0; i < periodicAssociations[-2]->getSize(); i++) {
- // MSG("asso %d : dof %d -> dof %d\n",
- // periodicAssociations[-2],
- // i,
- // (*periodicAssociations[-2])[i]);
- // }
-
- for (iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++)
- {
- TEST_EXIT((compressAdmin = admin[iadmin]))
- ("no admin[%d] in mesh\n", iadmin);
-
- if ((size = compressAdmin->getSize()) < 1) continue;
- if (compressAdmin->getUsedDOFs() < 1) continue;
- if (compressAdmin->getHoleCount() < 1) continue;
-
- newDOF.resize(size);
+ FUNCNAME("Mesh::dofCompress()");
+ int size;
+ Flag fill_flag;
- compressAdmin->compress(newDOF);
+ for (iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
+ TEST_EXIT((compressAdmin = admin[iadmin]))
+ ("no admin[%d] in mesh\n", iadmin);
+
+ if ((size = compressAdmin->getSize()) < 1)
+ continue;
+ if (compressAdmin->getUsedDOFs() < 1)
+ continue;
+ if (compressAdmin->getHoleCount() < 1)
+ continue;
- if (preserveCoarseDOFs) {
- fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
- }
- else {
- fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
- }
-
- traverse( -1, fill_flag, newDOFFct1);
- traverse( -1, fill_flag, newDOFFct2);
-
- newDOF.resize(0);
+ newDOF.resize(size);
+
+ compressAdmin->compress(newDOF);
+
+ if (preserveCoarseDOFs) {
+ fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
+ } else {
+ fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
}
-
- // MSG("\n");
- // for(i = 0; i < periodicAssociations[-2]->getSize(); i++) {
- // MSG("asso %d : dof %d -> dof %d\n",
- // periodicAssociations[-2],
- // i,
- // (*periodicAssociations[-2])[i]);
- // }
-
- return;
+
+ traverse( -1, fill_flag, newDOFFct1);
+ traverse( -1, fill_flag, newDOFFct2);
+
+ newDOF.resize(0);
+ }
}
- // /****************************************************************************/
- // /* fill_boundary: */
- // /* fills boundary information for the vertices of the macro triangulation */
- // /* mel_vertex[k][i] = global index of vertex k of element i */
- // /****************************************************************************/
-
- // void Mesh::fillBoundary(int **mel_vertex)
- // {
- // FUNCNAME("Mesh::fillBoundary");
- // ::std::deque<MacroElement*>::const_iterator mel = firstMacroElement();
- // int i, j,k;
- // int lne = getNumberOfLeaves(), lnv = getNumberOfVertices();
- // int ned = getNumberOfEdges();
-
- // BoundaryType *bound = GET_MEMORY(BoundaryType, lnv+(3==dim)?ned:0);
-
- // for (i = 0; i < lnv; i++)
- // bound[i] = INTERIOR;
-
- // int facesPlusEdges =
- // (*mel)->getElement()->getGeo(FACE) +
- // (*mel)->getElement()->getGeo(EDGE);
-
- // for (i = 0; i < lne; i++) {
- // for (k = 0; k < getGeo(NEIGH); k++) {
- // //if ((*(mel+i))->getBoundary(k)) {
- // if ((*(mel+i))->getBoundary(k) >= DIRICHLET) {
- // for (j = 1; j < dim+1; j++)
- // bound[mel_vertex[i][(k+j)%(dim+1)]] = DIRICHLET;
- // }
- // else if ((*(mel+i))->getBoundary(k) <= NEUMANN) {
- // for (j = 1; j < dim+1; j++)
- // if (bound[mel_vertex[i][(k+j)%(dim+1)]] != DIRICHLET)
- // bound[mel_vertex[i][(k+j)%(dim+1)]] = NEUMANN;
- // }
- // //}
- // }
- // }
-
- // for (i = 0; i < lne; i++)
- // for (k = 0; k < getGeo(VERTEX); k++) {
- // (*(mel+i))->setBoundary(facesPlusEdges + k,
- // bound[mel_vertex[i][k]]);
- // }
- // FREE_MEMORY(bound, BoundaryType, lnv+(3==dim)?ned:0);
-
- // return;
- // };
DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
{
- FUNCNAME("Mesh::getDOF");
- const DOFAdmin *localAdmin;
- DegreeOfFreedom *dof;
- int i, j, n, n0, ndof;
+ FUNCNAME("Mesh::getDOF()");
TEST_EXIT(position >= CENTER && position <= FACE)
- ("unknown position %d\n",position);
-
- ndof = getNumberOfDOFs(position);
- if (ndof <= 0) return(NULL);
-
- dof = GET_MEMORY(DegreeOfFreedom, ndof);
+ ("unknown position %d\n", position);
- for (i = 0; i < getNumberOfDOFAdmin(); i++)
- {
- localAdmin = &getDOFAdmin(i);
- TEST_EXIT(localAdmin)("no admin[%d]\n", i);
-
- n = localAdmin->getNumberOfDOFs(position);
- n0 = localAdmin->getNumberOfPreDOFs(position);
+ int ndof = getNumberOfDOFs(position);
+ if (ndof <= 0)
+ return(NULL);
- TEST_EXIT(n+n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
+ DegreeOfFreedom *dof = GET_MEMORY(DegreeOfFreedom, ndof);
- for (j = 0; j < n; j++)
- {
- dof[n0+j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
- }
+ for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
+ const DOFAdmin *localAdmin = &getDOFAdmin(i);
+ TEST_EXIT(localAdmin)("no admin[%d]\n", i);
+
+ int n = localAdmin->getNumberOfDOFs(position);
+ int n0 = localAdmin->getNumberOfPreDOFs(position);
+
+ TEST_EXIT(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
+
+ for (int j = 0; j < n; j++) {
+ dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
}
+ }
return(dof);
}
- // const Boundary *Mesh::defaultBoundary(int bound)
- // {
- // static const Boundary dn[2] = {DIRICHLET, NEUMANN};
- // if (bound >= DIRICHLET)
- // return(dn);
- // else if (bound <= NEUMANN)
- // return(dn+1);
- // else
- // return(NULL);
- // }
-
DegreeOfFreedom **Mesh::createDOFPtrs()
{
- FUNCNAME("Mesh::createDOFPtrs");
- int i;
- DegreeOfFreedom **ptrs;
+ FUNCNAME("Mesh::createDOFPtrs()");
if (nNodeEl <= 0)
return(NULL);
- ptrs = GET_MEMORY(DegreeOfFreedom*, nNodeEl);
- for (i = 0; i < nNodeEl; i++)
+ DegreeOfFreedom **ptrs = GET_MEMORY(DegreeOfFreedom*, nNodeEl);
+ for (int i = 0; i < nNodeEl; i++)
ptrs[i] = NULL;
return(ptrs);
@@ -525,7 +443,7 @@ namespace AMDiS {
void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
{
- FUNCNAME("Mesh::freeDOFPtrs");
+ FUNCNAME("Mesh::freeDOFPtrs()");
TEST_EXIT(ptrs)("ptrs=NULL\n");
@@ -650,42 +568,23 @@ namespace AMDiS {
}
-
-
- //const Boundary* defaultBoundary(Mesh* m,int i) {return m->defaultBoundary(i);};
-
Element* Mesh::createNewElement(Element *parent)
{
FUNCNAME("Mesh::createNewElement()");
TEST_EXIT(elementPrototype)("no element prototype\n");
Element *el = parent ? parent->clone() : elementPrototype->clone();
-
- //el->setIndex(elementIndex++);
- //el->setDOFPtrs();
-
- // ElementData *elementData = NULL;
-
- // if(parent && parent->getElementData()) {
- // elementData = parent->getElementData()->clone();
- // } else if(elementDataPrototype) {
- // elementData = elementDataPrototype->clone();
- // }
- if(!parent && elementDataPrototype) {
+ if (!parent && elementDataPrototype) {
el->setElementData(elementDataPrototype->clone());
} else {
el->setElementData(NULL); // must be done in ElementData::refineElementData()
}
- // if(!elementData) {
- // WARNING("no element data!\n");
- // }
-
-
return el;
}
+
ElInfo* Mesh::createNewElInfo()
{
switch(dim) {
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 50920aa5cf72442fad399db498b4e95b3a3b032a..cddd6c8dffd725d45395aec29c7bddc1ae844d9e 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -263,11 +263,10 @@ namespace AMDiS {
if (feSpace_) {
WARNING("feSpace already created\n");
} else {
- if (initFlag.isSet(INIT_FE_SPACE) || (initFlag.isSet(INIT_SYSTEM)&&!adoptFlag.isSet(INIT_FE_SPACE))) {
+ if (initFlag.isSet(INIT_FE_SPACE) || (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
createFESpace();
}
- if (adoptProblem &&
- (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM)))
+ if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM)))
{
TEST_EXIT(!feSpace_)("feSpace already created");
feSpace_ = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace();
@@ -554,7 +553,7 @@ namespace AMDiS {
void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag)
{
- FUNCNAME("ProblemScal::buildAfterCoarsen");
+ FUNCNAME("ProblemScal::buildAfterCoarsen()");
clock_t first = clock();
diff --git a/AMDiS/src/VtkWriter.cc b/AMDiS/src/VtkWriter.cc
index f362f60a353435a33de3cadd7e19e914474cfcac..df47b2ab174ccb5f3260f6d443be7d3bd80c4e54 100644
--- a/AMDiS/src/VtkWriter.cc
+++ b/AMDiS/src/VtkWriter.cc
@@ -24,6 +24,9 @@ namespace AMDiS {
} else if ((dim_ == 2) && (degree_ == 3)) {
nVertices += (*dc_)[0]->getNumberInterpPoints();
nElements *= 9;
+ } else if ((dim_ == 2) && (degree_ == 4)) {
+ nVertices += (*dc_)[0]->getNumberInterpPoints();
+ nElements *= 16;
}
::std::ofstream file;
@@ -70,12 +73,8 @@ namespace AMDiS {
file << " </DataArray>" << ::std::endl;
file << " <DataArray type=\"Int32\" Name=\"connectivity\">" << ::std::endl;
- if ((dim_ == 2) && (degree_ == 2)) {
- writeConnectivity_dim2_degree2(file);
- } else {
- writeConnectivity(file);
- }
-
+ writeConnectivity(file);
+
file << " </DataArray>" << ::std::endl;
file << " </Cells>" << ::std::endl;
file << " <PointData>" << ::std::endl;
@@ -171,17 +170,25 @@ namespace AMDiS {
void VtkWriter::writeConnectivity(::std::ofstream &file)
{
- ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
- ::std::list<ElementInfo>::iterator elementIt;
- int vertices = (*dc_)[0]->getMesh()->getGeo(VERTEX);
-
- for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
- // for all vertices
- for (int i = 0; i < vertices; i++) {
- file << " " << elementIt->vertexInfo[i]->outputIndex;
- }
- file << ::std::endl;
- }
+ if ((dim_ == 2) && (degree_ == 2)) {
+ writeConnectivity_dim2_degree2(file);
+ } else if ((dim_ == 2) && (degree_ == 3)) {
+ writeConnectivity_dim2_degree3(file);
+ } else if ((dim_ == 2) && (degree_ == 4)) {
+ writeConnectivity_dim2_degree4(file);
+ } else {
+ ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
+ ::std::list<ElementInfo>::iterator elementIt;
+ int vertices = (*dc_)[0]->getMesh()->getGeo(VERTEX);
+
+ for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
+ // for all vertices
+ for (int i = 0; i < vertices; i++) {
+ file << " " << elementIt->vertexInfo[i]->outputIndex;
+ }
+ file << ::std::endl;
+ }
+ }
}
@@ -219,6 +226,143 @@ namespace AMDiS {
}
+ void VtkWriter::writeConnectivity_dim2_degree3(::std::ofstream &file)
+ {
+ ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
+ ::std::list<ElementInfo>::iterator elementIt;
+
+ ::std::vector< ::std::vector<int> > *interpPoints = (*dc_)[0]->getInterpPoints();
+ ::std::vector< ::std::vector<int> >::iterator pointIt;
+
+ int nVertices = (*dc_)[0]->getNumberVertices();
+
+ for (pointIt = interpPoints->begin(), elementIt = elements->begin();
+ pointIt != interpPoints->end();
+ ++pointIt, ++elementIt) {
+
+ file << " " << elementIt->vertexInfo[0]->outputIndex
+ << " " << (*pointIt)[3] + nVertices
+ << " " << (*pointIt)[4] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[4] + nVertices
+ << " " << (*pointIt)[5] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[3] + nVertices
+ << " " << (*pointIt)[4] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[2] + nVertices
+ << " " << (*pointIt)[3] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << elementIt->vertexInfo[1]->outputIndex
+ << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[5] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[6] + nVertices
+ << " " << (*pointIt)[5] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[2] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << elementIt->vertexInfo[2]->outputIndex
+ << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[2] + nVertices << ::std::endl;
+
+
+ }
+
+ }
+
+
+ void VtkWriter::writeConnectivity_dim2_degree4(::std::ofstream &file)
+ {
+ ::std::list<ElementInfo> *elements = (*dc_)[0]->getElementInfos();
+ ::std::list<ElementInfo>::iterator elementIt;
+
+ ::std::vector< ::std::vector<int> > *interpPoints = (*dc_)[0]->getInterpPoints();
+ ::std::vector< ::std::vector<int> >::iterator pointIt;
+
+ int nVertices = (*dc_)[0]->getNumberVertices();
+
+ for (pointIt = interpPoints->begin(), elementIt = elements->begin();
+ pointIt != interpPoints->end();
+ ++pointIt, ++elementIt) {
+
+ file << " " << elementIt->vertexInfo[0]->outputIndex
+ << " " << (*pointIt)[5] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[5] + nVertices
+ << " " << (*pointIt)[9] + nVertices
+ << " " << (*pointIt)[6] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[6] + nVertices
+ << " " << (*pointIt)[7] + nVertices
+ << " " << (*pointIt)[9] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[7] + nVertices
+ << " " << (*pointIt)[9] + nVertices
+ << " " << (*pointIt)[10] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[7] + nVertices
+ << " " << (*pointIt)[8] + nVertices
+ << " " << (*pointIt)[10] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[8] + nVertices
+ << " " << (*pointIt)[10] + nVertices << ::std::endl;
+
+ file << " " << elementIt->vertexInfo[1]->outputIndex
+ << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[8] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[4] + nVertices
+ << " " << (*pointIt)[5] + nVertices
+ << " " << (*pointIt)[9] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[4] + nVertices
+ << " " << (*pointIt)[9] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[9] + nVertices
+ << " " << (*pointIt)[10] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[10] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[0] + nVertices
+ << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[10] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[3] + nVertices
+ << " " << (*pointIt)[4] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[2] + nVertices
+ << " " << (*pointIt)[3] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << (*pointIt)[1] + nVertices
+ << " " << (*pointIt)[2] + nVertices
+ << " " << (*pointIt)[11] + nVertices << ::std::endl;
+
+ file << " " << elementIt->vertexInfo[2]->outputIndex
+ << " " << (*pointIt)[2] + nVertices
+ << " " << (*pointIt)[3] + nVertices << ::std::endl;
+
+ }
+ }
+
int VtkWriter::updateAnimationFile(::std::string valueFilename,
::std::vector< ::std::string > *paraViewAnimationFrames,
const char *animationFilename)
diff --git a/AMDiS/src/VtkWriter.h b/AMDiS/src/VtkWriter.h
index b8aa1db122ff05227e39a12319a0ace1ee15d8e6..3de2e8806384b4ef377be0c776a00ad9f56304f0 100644
--- a/AMDiS/src/VtkWriter.h
+++ b/AMDiS/src/VtkWriter.h
@@ -80,6 +80,18 @@ namespace AMDiS {
void writeConnectivity_dim2_degree2(::std::ofstream &file);
+ /** \brief
+ *
+ */
+ void writeConnectivity_dim2_degree3(::std::ofstream &file);
+
+
+ /** \brief
+ *
+ */
+ void writeConnectivity_dim2_degree4(::std::ofstream &file);
+
+
/** \brief
*
*/