Commit 6640e0b5 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* VTK Output possible for Lagrange elements of second order

parent 8f323098
...@@ -96,14 +96,14 @@ namespace AMDiS { ...@@ -96,14 +96,14 @@ namespace AMDiS {
* compares two BasisFunction objects. * compares two BasisFunction objects.
*/ */
virtual bool operator==(const BasisFunction& a) const { virtual bool operator==(const BasisFunction& a) const {
return a.getName()==name; return a.getName() == name;
}; };
/** \brief /** \brief
* Returns !(*this == b) * Returns !(*this == b)
*/ */
inline bool operator!=(const BasisFunction& b) const { inline bool operator!=(const BasisFunction& b) const {
return !operator==(b); return !(operator == (b));
}; };
/** \brief /** \brief
......
...@@ -303,18 +303,21 @@ namespace AMDiS { ...@@ -303,18 +303,21 @@ namespace AMDiS {
} }
const int DOFAdmin::getNumberOfPreDOFs(int i) const { const int DOFAdmin::getNumberOfPreDOFs(int i) const {
TEST_EXIT((0<=i)&&(4>i))(""); TEST_EXIT((0 <= i) && (4 > i))("");
return nr0DOF[i]; return nr0DOF[i];
} }
void DOFAdmin::setNumberOfDOFs(int i,int v) { void DOFAdmin::setNumberOfDOFs(int i,int v) {
TEST_EXIT((0<=i)&&(4>i))(""); TEST_EXIT((0 <= i) && (4 > i))("");
nrDOF[i]=v;
nrDOF[i] = v;
} }
void DOFAdmin::setNumberOfPreDOFs(int i, int v) { void DOFAdmin::setNumberOfPreDOFs(int i, int v) {
TEST_EXIT((0<=i)&&(4>i))(""); TEST_EXIT((0 <= i) && (4 > i))("");
nr0DOF[i]=v;
nr0DOF[i] = v;
} }
DOFAdmin::~DOFAdmin() DOFAdmin::~DOFAdmin()
......
...@@ -13,7 +13,7 @@ namespace AMDiS { ...@@ -13,7 +13,7 @@ namespace AMDiS {
bool (*writeElem)(ElInfo*)) bool (*writeElem)(ElInfo*))
: writeElem_(writeElem) : writeElem_(writeElem)
{ {
FUNCNAME("DataCollector::DataCollector"); FUNCNAME("DataCollector::DataCollector()");
TEST_EXIT(feSpace)("no feSpace\n"); TEST_EXIT(feSpace)("no feSpace\n");
TEST_EXIT(values)("no value Vector\n"); TEST_EXIT(values)("no value Vector\n");
...@@ -27,6 +27,7 @@ namespace AMDiS { ...@@ -27,6 +27,7 @@ namespace AMDiS {
// === create vertex info vector === // === create vertex info vector ===
vertexInfos_ = NEW DOFVector< ::std::list<VertexInfo> >(feSpace, "vertex infos"); vertexInfos_ = NEW DOFVector< ::std::list<VertexInfo> >(feSpace, "vertex infos");
interpPointInd_ = NEW DOFVector<int>(feSpace, "interpolation point indices"); interpPointInd_ = NEW DOFVector<int>(feSpace, "interpolation point indices");
interpPointCoords_ = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
dofCoords_ = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "dof coords"); dofCoords_ = NEW DOFVector< ::std::list<WorldVector<double> > >(feSpace, "dof coords");
dim_ = mesh_->getDim(); dim_ = mesh_->getDim();
...@@ -51,12 +52,13 @@ namespace AMDiS { ...@@ -51,12 +52,13 @@ namespace AMDiS {
{ {
DELETE vertexInfos_; DELETE vertexInfos_;
DELETE interpPointInd_; DELETE interpPointInd_;
DELETE interpPointCoords_;
DELETE dofCoords_; DELETE dofCoords_;
} }
int DataCollector::startCollectingElementData() int DataCollector::startCollectingElementData()
{ {
FUNCNAME("DataCollector::startCollectingElementData"); FUNCNAME("DataCollector::startCollectingElementData()");
Flag flag = traverseFlag_; Flag flag = traverseFlag_;
flag |= flag |=
...@@ -71,16 +73,16 @@ namespace AMDiS { ...@@ -71,16 +73,16 @@ namespace AMDiS {
// Traverse elements to create continuous element indices // Traverse elements to create continuous element indices
ElInfo *elInfo = stack.traverseFirst(mesh_, level_, flag); ElInfo *elInfo = stack.traverseFirst(mesh_, level_, flag);
while(elInfo) { while (elInfo) {
if(!writeElem_ || writeElem_(elInfo)) if (!writeElem_ || writeElem_(elInfo))
outputIndices_[elInfo->getElement()->getIndex()] = nElements_++; outputIndices_[elInfo->getElement()->getIndex()] = nElements_++;
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
// Traverse elements to create element information // Traverse elements to create element information
elInfo = stack.traverseFirst(mesh_, level_, flag); elInfo = stack.traverseFirst(mesh_, level_, flag);
while(elInfo) { while (elInfo) {
if(!writeElem_ || writeElem_(elInfo)) if (!writeElem_ || writeElem_(elInfo))
addElementData(elInfo); addElementData(elInfo);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
...@@ -92,10 +94,10 @@ namespace AMDiS { ...@@ -92,10 +94,10 @@ namespace AMDiS {
int DataCollector::startCollectingValueData() int DataCollector::startCollectingValueData()
{ {
FUNCNAME("DataCollector::startCollectingValueData"); FUNCNAME("DataCollector::startCollectingValueData()");
DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS); DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS);
for(intPointIt.reset(); !intPointIt.end(); ++intPointIt) { for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
(*intPointIt) = -1; (*intPointIt) = -1;
} }
...@@ -109,7 +111,7 @@ namespace AMDiS { ...@@ -109,7 +111,7 @@ namespace AMDiS {
level_, level_,
traverseFlag_ | Mesh::FILL_COORDS); traverseFlag_ | Mesh::FILL_COORDS);
while(elInfo) { while(elInfo) {
if(!writeElem_ || writeElem_(elInfo)) if (!writeElem_ || writeElem_(elInfo))
addValueData(elInfo); addValueData(elInfo);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
...@@ -117,16 +119,16 @@ namespace AMDiS { ...@@ -117,16 +119,16 @@ namespace AMDiS {
// Remove all interpolation marks and, instead, set to each // Remove all interpolation marks and, instead, set to each
// interpolation point its continous index starting from 0. // interpolation point its continous index starting from 0.
int i = 0; int i = 0;
for(intPointIt.reset(); !intPointIt.end(); ++intPointIt) { for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
if(*intPointIt == -3) { if (*intPointIt == -3) {
*intPointIt = i++; *intPointIt = i++;
} }
} }
// Traverse elements to create interpolation values. // Traverse elements to create interpolation values.
elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_); elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
while(elInfo) { while (elInfo) {
if(!writeElem_ || writeElem_(elInfo)) if (!writeElem_ || writeElem_(elInfo))
addInterpData(elInfo); addInterpData(elInfo);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
...@@ -138,22 +140,21 @@ namespace AMDiS { ...@@ -138,22 +140,21 @@ namespace AMDiS {
int DataCollector::startCollectingPeriodicData() int DataCollector::startCollectingPeriodicData()
{ {
FUNCNAME("DataCollector::startCollectingPeriodicData"); FUNCNAME("DataCollector::startCollectingPeriodicData()");
periodicConnections_.clear(); periodicConnections_.clear();
TraverseStack stack; TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_); ElInfo *elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
while(elInfo) { while (elInfo) {
if(!writeElem_ || writeElem_(elInfo)) { if (!writeElem_ || writeElem_(elInfo)) {
LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*> LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
(elInfo->getElement()-> (elInfo->getElement()->
getElementData()-> getElementData()->
getElementData(PERIODIC)); getElementData(PERIODIC));
if(ldp) { if (ldp) {
nConnections_ += nConnections_ += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
} }
periodicConnections_.push_back(DimVec<bool>(dim_, DEFAULT_VALUE, false)); periodicConnections_.push_back(DimVec<bool>(dim_, DEFAULT_VALUE, false));
...@@ -173,8 +174,8 @@ namespace AMDiS { ...@@ -173,8 +174,8 @@ namespace AMDiS {
Mesh::FILL_BOUND; Mesh::FILL_BOUND;
elInfo = stack.traverseFirst(mesh_, level_, flag); elInfo = stack.traverseFirst(mesh_, level_, flag);
while(elInfo) { while (elInfo) {
if(!writeElem_ || writeElem_(elInfo)) if (!writeElem_ || writeElem_(elInfo))
addPeriodicData(elInfo); addPeriodicData(elInfo);
elInfo = stack.traverseNext(elInfo); elInfo = stack.traverseNext(elInfo);
} }
...@@ -186,9 +187,8 @@ namespace AMDiS { ...@@ -186,9 +187,8 @@ namespace AMDiS {
int DataCollector::addElementData(ElInfo* elInfo) int DataCollector::addElementData(ElInfo* elInfo)
{ {
FUNCNAME("DataCollector::addElementData"); FUNCNAME("DataCollector::addElementData()");
int i;
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
DegreeOfFreedom vertexDOF; DegreeOfFreedom vertexDOF;
WorldVector<double> vertexCoords; WorldVector<double> vertexCoords;
...@@ -197,10 +197,9 @@ namespace AMDiS { ...@@ -197,10 +197,9 @@ namespace AMDiS {
ElementInfo elementInfo(dim_); ElementInfo elementInfo(dim_);
// read element region // read element region
ElementData *ed = ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
elInfo->getElement()->getElementData(ELEMENT_REGION);
if(ed) { if (ed) {
elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion(); elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
} else { } else {
elementInfo.elementRegion = -1; elementInfo.elementRegion = -1;
...@@ -211,14 +210,14 @@ namespace AMDiS { ...@@ -211,14 +210,14 @@ namespace AMDiS {
elementInfo.surfaceRegions.set(-1); elementInfo.surfaceRegions.set(-1);
while(ed) { while (ed) {
SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed); SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion(); elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
ed = ed->getDecorated(SURFACE_REGION); ed = ed->getDecorated(SURFACE_REGION);
} }
// for all vertices // for all vertices
for(i = 0; i < dim_ + 1; i++) { for (int i = 0; i < dim_ + 1; i++) {
// get coords of this vertex // get coords of this vertex
vertexCoords = elInfo->getCoord(i); vertexCoords = elInfo->getCoord(i);
...@@ -232,7 +231,7 @@ namespace AMDiS { ...@@ -232,7 +231,7 @@ namespace AMDiS {
vertexCoords); vertexCoords);
// coords not yet in list? // coords not yet in list?
if(it == (*vertexInfos_)[vertexDOF].end()) { if (it == (*vertexInfos_)[vertexDOF].end()) {
// create new vertex info and increment number of vertices // create new vertex info and increment number of vertices
VertexInfo newVertexInfo = {vertexCoords, nVertices_}; VertexInfo newVertexInfo = {vertexCoords, nVertices_};
...@@ -255,7 +254,7 @@ namespace AMDiS { ...@@ -255,7 +254,7 @@ namespace AMDiS {
-1; -1;
} }
if(dim_ == 3) { if (dim_ == 3) {
elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType()); elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
} }
...@@ -267,16 +266,14 @@ namespace AMDiS { ...@@ -267,16 +266,14 @@ namespace AMDiS {
int DataCollector::addValueData(ElInfo *elInfo) int DataCollector::addValueData(ElInfo *elInfo)
{ {
FUNCNAME("DataCollector::addValueData"); FUNCNAME("DataCollector::addValueData()");
int i, j, k;
int node_offset, dof_offset, num_dofs, node, dof_index;
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
/* vertex dofs */ /* vertex dofs */
dof_offset = localAdmin_->getNumberOfPreDOFs(VERTEX); int dof_offset = localAdmin_->getNumberOfPreDOFs(VERTEX);
for (i = 0; i < mesh_->getGeo(VERTEX); i++) { for (int i = 0; i < mesh_->getGeo(VERTEX); i++) {
(*interpPointInd_)[dof[i][dof_offset]] = -2; // mark as vertex (*interpPointInd_)[dof[i][dof_offset]] = -2; // mark as vertex
// get coords of this vertex // get coords of this vertex
...@@ -289,53 +286,76 @@ namespace AMDiS { ...@@ -289,53 +286,76 @@ namespace AMDiS {
vertexCoords); vertexCoords);
// coords not yet in list? // coords not yet in list?
if(it == (*dofCoords_)[dof[i][dof_offset]].end()) { if (it == (*dofCoords_)[dof[i][dof_offset]].end()) {
// add new coords to list // add new coords to list
(*dofCoords_)[dof[i][dof_offset]].push_back(vertexCoords); (*dofCoords_)[dof[i][dof_offset]].push_back(vertexCoords);
} }
} }
for(i = 1; i <= dim_; i++) { int nInterpPoints = 0;
num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_)); const BasisFunction *basisFcts = feSpace_->getBasisFcts();
node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_)); for (int i = 1; i <= dim_; i++) {
for (j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) { int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
node = node_offset + j; int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
for(k = 0; k < num_dofs; k++) { dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
dof_index = dof_offset + k;
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) { if ((*interpPointInd_)[dof[node][dof_index]] == -1) {
// mark as interp. point // mark as interpolation point
(*interpPointInd_)[dof[node][dof_index]] = -3; (*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(),
interpolCoords);
if (it == (*interpPointCoords_)[dof[node][dof_index]].end()) {
(*interpPointCoords_)[dof[node][dof_index]].push_back(interpolCoords);
}
nInterpPoints_++; nInterpPoints_++;
} }
} }
} }
} }
return(0); return(0);
} }
int DataCollector::addInterpData(ElInfo *elInfo) int DataCollector::addInterpData(ElInfo *elInfo)
{ {
FUNCNAME("DataCollector::addInterpData"); FUNCNAME("DataCollector::addInterpData()");
int i, j, k;
int node_offset, dof_offset, num_dofs, node, dof_index;
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF(); const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
::std::list<int> elemInterpPoints; ::std::vector<int> elemInterpPoints;
elemInterpPoints.clear(); elemInterpPoints.clear();
for(i = 1; i <= dim_; i++) { for (int i = 1; i <= dim_; i++) {
num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_)); int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_)); int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_)); int dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
for (j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
node = node_offset + j; for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
for(k = 0; k < num_dofs; k++) { int node = node_offset + j;
dof_index = dof_offset + k;
for (int k = 0; k < num_dofs; k++) {
int dof_index = dof_offset + k;
elemInterpPoints.push_back((*interpPointInd_)[dof[node][dof_index]]); elemInterpPoints.push_back((*interpPointInd_)[dof[node][dof_index]]);
} }
...@@ -355,57 +375,54 @@ namespace AMDiS { ...@@ -355,57 +375,54 @@ namespace AMDiS {
getElementData()-> getElementData()->
getElementData(PERIODIC)); getElementData(PERIODIC));
if(ldp) { if (ldp) {
::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it; ::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
for(it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin(); for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end(); it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
++it) ++it) {
{
int outputIndex = outputIndices_[elInfo->getElement()->getIndex()]; int outputIndex = outputIndices_[elInfo->getElement()->getIndex()];
int neighIndex = outputIndices_[elInfo-> int neighIndex = outputIndices_[elInfo->
getNeighbour(it->elementSide)-> getNeighbour(it->elementSide)->
getIndex()]; getIndex()];
if (!periodicConnections_[outputIndex][it->elementSide]) {
PeriodicInfo periodicInfo;
if(!periodicConnections_[outputIndex][it->elementSide]) { periodicInfo.mode = it->periodicMode;
PeriodicInfo periodicInfo; periodicInfo.type = it->type;
periodicInfo.outputIndex = outputIndex;
periodicInfo.mode = it->periodicMode; periodicInfo.neighIndex = neighIndex;
periodicInfo.type = it->type; periodicInfo.vertexMap.clear();
periodicInfo.outputIndex = outputIndex;
periodicInfo.neighIndex = neighIndex; periodicConnections_[outputIndex][it->elementSide] = true;
periodicInfo.vertexMap.clear(); periodicConnections_
[neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
periodicConnections_[outputIndex][it->elementSide] = true;
periodicConnections_
[neighIndex][elInfo->getOppVertex(it->elementSide)] = true; for (int i = 0; i < dim_; i++) {
int index1 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
int index1, index2, dof1, dof2, i, j; it->elementSide,
i);
for(i = 0; i < dim_; i++) { int dof1 = elInfo->getElement()->getDOF(index1, nPreDofs_);
index1 =
elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_), for (int j = 0; j < dim_; j++) {
it->elementSide, int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
i); elInfo->getOppVertex(it->elementSide),
dof1 = elInfo->getElement()->getDOF(index1, nPreDofs_); j);
int dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs_);
for(j = 0; j < dim_; j++) {
index2 = if ((dof1 == dof2) || (mesh_->associated(dof1, dof2))) {
elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_), periodicInfo.vertexMap[index1] = index2;
elInfo->getOppVertex(it->elementSide), break;
j); }
dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs_);
if((dof1 == dof2) || (mesh_->associated(dof1, dof2))) {
periodicInfo.vertexMap[index1] = index2;
break;
}
}
} }
periodicInfos_.push_back(periodicInfo);
} }
periodicInfos_.push_back(periodicInfo);
} }
}
} }
return(0); return(0);
...@@ -511,7 +528,16 @@ namespace AMDiS { ...@@ -511,7 +528,16 @@ namespace AMDiS {
return interpPointInd_; return interpPointInd_;
} }
::std