merge of two FileWriter versions

AMDiS/CMakeLists.txt

@@ -150,7 +150,6 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/est/SimpleResidualEstimator.cc
 	      ${SOURCE_DIR}/io/ArhReader.cc
 	      ${SOURCE_DIR}/io/ArhWriter.cc
-	      ${SOURCE_DIR}/io/DataCollector.cc
 	      ${SOURCE_DIR}/io/DofWriter.cc
 	      ${SOURCE_DIR}/io/ElementFileWriter.cc
 	      ${SOURCE_DIR}/io/FileWriter.cc

AMDiS/src/AMDiS.h

@@ -124,6 +124,10 @@
 #include "io/ValueReader.h"
 #include "io/ValueWriter.h"
 #include "io/VtkWriter.h"
+#include "io/VtkVectorWriter.h"
+
+#include "io/DataCollector.hh"
+#include "io/VtkVectorWriter.hh"
 
 #include "nonlin/ProblemNonLin.h"
 #include "nonlin/NonLinSolver.h"

AMDiS/src/CouplingProblemStat.h

@@ -110,7 +110,7 @@ namespace AMDiS {
                           refSet);
           if (refSet < 0)
             refSet = 0;
-          
+
           if (meshForRefinementSet.find(refSet) == meshForRefinementSet.end()) {
             Mesh *newMesh = new Mesh(meshName, dim);
             meshForRefinementSet[refSet] = newMesh;
@@ -126,7 +126,7 @@ namespace AMDiS {
           int degree = 1;
           Parameters::get(problems[i]->getName() + "->polynomial degree[" + 
                           boost::lexical_cast<string>(j) + "]", degree);
-    
+
           if (feSpaceMap[pair<Mesh*, int>(meshForRefinementSet[refSet], degree)] == NULL) {
             stringstream s;
             s << problems[i]->getName() << "->feSpace[" << j << "]";

AMDiS/src/Initfile.cc

@@ -145,10 +145,14 @@ namespace AMDiS {
         posVarBegin = posVar;
       }
       std::string varName = inputSwap.substr(posVarBegin + 1 , posVarEnd - posVarBegin - 1);
-      std::string varParam = "";
-      int error_code = checkedGet(varName, varParam);
+
       // if varname is found in parameter list then replace variable by value
       // otherwise throw tagNotFound exception
+      std::string varParam;
+      int error_code = checkedGet(varName, varParam);
+      if (error_code != 0)
+	throw TagNotFoundBreak("required tag '" + varName + "' for variable substitution not found");
+      
       std::string replaceName = inputSwap.substr(posVar , posVarEnd - posVar + (posVarBegin - posVar));
       inputSwap.replace(inputSwap.find(replaceName), replaceName.length(), varParam);
 

AMDiS/src/Initfile.h

@@ -446,8 +446,11 @@ namespace AMDiS {
 	convert(valStr, value);
       } else if(error_code == TAG_NOT_FOUND_BREAK)
 	throw TagNotFoundBreak("required tag '" + tag + "' not found");
-      else if (error_code == TAG_NOT_FOUND && debugInfo == 2)
-	std::cout << "there is no tag '" + tag + "'" << std::endl;
+      else if (error_code == TAG_NOT_FOUND) {
+	if (debugInfo == 2)
+	  std::cout << "there is no tag '" + tag + "'" << std::endl;
+      } else
+	throw std::runtime_error("unknown error_code (" + boost::lexical_cast<std::string>(error_code) + ") in checkedGet(...) returned for tag '" + tag + "'");
       
       if (debugInfo == 2) {
 	std::cout << "Parameter '" << tag << "'"
@@ -460,7 +463,6 @@ namespace AMDiS {
     static InitEntry get(const std::string tag) 
     {
       using namespace InitfileInternal;
-      singlett->getMsgInfo();
       InitEntry result;
       
       std::string valStr;
@@ -472,6 +474,8 @@ namespace AMDiS {
 	throw TagNotFoundBreak("required tag '" + tag + "' not found");
       else if (error_code == TAG_NOT_FOUND)
 	throw TagNotFound("there is no tag '" + tag + "'"); // exception must be thrown, because an empty object would be return otherwise
+      else
+	throw std::runtime_error("unknown error_code in checkedGet(...) returned for tag '" + tag + "'");
 	
       return result;
     }
@@ -592,7 +596,6 @@ protected:
 	  return TAG_NOT_FOUND_BREAK;
       }
       valStr = it->second;
-
       return 0;
     }
 

AMDiS/src/MatrixVector.h

@@ -432,8 +432,8 @@ namespace AMDiS {
   }
 
   /// y = a * x + y.
-  template<typename T>
-  inline const Vector<T>& axpy(const T& a,
+  template<typename T, typename S>
+  inline const Vector<T>& axpy(const S& a,
 			       const Vector<T> &x,
 			       Vector<T> &y)
   {

AMDiS/src/SubAssembler.hh

@@ -74,6 +74,7 @@ namespace AMDiS {
 
     cachedValuesAtQPs[vec]->valid = true;
     vecAtQPs = values;
+
   }
 
 

AMDiS/src/io/DataCollector.cc

-//
-// Software License for AMDiS
-//
-// Copyright (c) 2010 Dresden University of Technology 
-// All rights reserved.
-// Authors: Simon Vey, Thomas Witkowski et al.
-//
-// This file is part of AMDiS
-//
-// See also license.opensource.txt in the distribution.
-
-
-#include "DataCollector.h"
-#include "Traverse.h"
-#include "DOFVector.h"
-#include "SurfaceRegion_ED.h"
-#include "ElementRegion_ED.h"
-#include "Projection.h"
-
-namespace AMDiS {
-
-  DataCollector::DataCollector(const FiniteElemSpace *fe,
-			       DOFVector<double> *val,
-			       int l,
-			       Flag flag,
-			       bool (*writeElemFct)(ElInfo*))
-    : values(val),
-      level(l),
-      traverseFlag(flag),
-      feSpace(fe),
-      nVertices(0),
-      nElements(0),
-      nInterpPoints(0),
-      nConnection(0),
-      elementDataCollected(false),
-      valueDataCollected(false),
-      periodicDataCollected(false),
-      writeElem(writeElemFct)
-  {
-    FUNCNAME("DataCollector::DataCollector()");
-
-    TEST_EXIT(feSpace)("No finite elem space defined!\n");
-   
-    // === get mesh  ===    
-    mesh = feSpace->getMesh();
-
-    // === get admin ===    
-    localAdmin = const_cast<DOFAdmin*>(feSpace->getAdmin());
-    // === create vertex info vector ===
-    vertexInfos = new DOFVector<std::list<VertexInfo> >(feSpace, "vertex infos");
-
-    interpPointInd = new DOFVector<int>(feSpace, "interpolation point indices");
-    interpPointCoords = new DOFVector< std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
-    dofCoord = new DOFVector< std::list<WorldVector<double> > >(feSpace, "dof coords");
-
-    dim = mesh->getDim();    
-    nPreDofs = localAdmin->getNumberOfPreDofs(VERTEX);
-  } 
-
-
-  DataCollector::~DataCollector()
-  {
-    delete vertexInfos;
-    delete interpPointInd;
-    delete interpPointCoords;
-    delete dofCoord;
-  }
-
-
-  void DataCollector::fillAllData()
-  {
-    if (!elementDataCollected)
-      startCollectingElementData();
-
-    if (!periodicDataCollected)
-      startCollectingPeriodicData();
-
-    if (!valueDataCollected)
-      startCollectingValueData();
-  }
-
-
-  void DataCollector::startCollectingElementData()
-  {
-    FUNCNAME("DataCollector::startCollectingElementData()");
-
-    Flag flag = traverseFlag;
-    flag |= 
-      Mesh::FILL_NEIGH      | 
-      Mesh::FILL_COORDS     |
-      Mesh::FILL_OPP_COORDS |
-      Mesh::FILL_BOUND;
-
-    elements.resize(0, dim);
-
-    TraverseStack stack;
-
-    // Traverse elements to create continuous element indices
-    ElInfo *elInfo = stack.traverseFirst(mesh, level, flag);
-    while (elInfo) {
-      if (!writeElem || writeElem(elInfo))
-	outputIndices[elInfo->getElement()->getIndex()] = nElements++;
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    // Traverse elements to create element information
-    elInfo = stack.traverseFirst(mesh, level, flag);
-
-    while (elInfo) {
-      if (!writeElem || writeElem(elInfo))
-	addElementData(elInfo);
-
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    elementDataCollected = true;
-  }
-
-
-  void DataCollector::startCollectingValueData()
-  {
-    FUNCNAME("DataCollector::startCollectingValueData()");
-
-    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
-    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
-      (*intPointIt) = -1;
-
-    interpPoints.clear();
-
-    basisFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
-    nBasisFcts = basisFcts->getNumber();
-    localDOFs = new DegreeOfFreedom[nBasisFcts];
-  
-    TraverseStack stack;
-
-    // Traverse elements to add value information and to mark
-    // interpolation points.
-    ElInfo *elInfo = stack.traverseFirst(mesh, level, 
-					 traverseFlag | Mesh::FILL_COORDS);
-    while (elInfo) {
-      if (!writeElem || writeElem(elInfo))
-	addValueData(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    // If there are interpolation points, add them to the corresponding 
-    // data array.
-    if (nInterpPoints > 0) {
-      // Remove all interpolation marks and, instead, set to each
-      // interpolation point its continous index starting from 0.
-      int i = 0;
-      for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
-	if (*intPointIt == -3)
-	  *intPointIt = i++;
-      
-      // Traverse elements to create interpolation values.
-      elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
-      while (elInfo) {
-	if (!writeElem || writeElem(elInfo))
-	  addInterpData(elInfo);
-	elInfo = stack.traverseNext(elInfo);
-      }
-    }
-   
-    delete [] localDOFs;
-    valueDataCollected = true;
-  }
-
-
-  void DataCollector::startCollectingPeriodicData()
-  {    
-    FUNCNAME("DataCollector::startCollectingPeriodicData()");
-
-    periodicConnections.clear();
-    
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
-    while (elInfo) {
-      if (!writeElem || writeElem(elInfo)) {
-	LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
-	  (elInfo->getElement()->
-	   getElementData()->
-	   getElementData(PERIODIC));
-	
-	if (ldp)
-	  nConnection += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
-
-	periodicConnections.push_back(DimVec<bool>(dim, DEFAULT_VALUE, false));
-      }
-      elInfo = stack.traverseNext(elInfo);
-    }   
-
-    nConnection /= 2;
-
-    periodicInfos.clear();
-
-    Flag flag = traverseFlag;
-    flag |= 
-      Mesh::FILL_COORDS    |
-      Mesh::FILL_OPP_COORDS|
-      Mesh::FILL_NEIGH     | 
-      Mesh::FILL_BOUND;
-
-    elInfo = stack.traverseFirst(mesh, level, flag);
-    while (elInfo) {
-      if (!writeElem || writeElem(elInfo))
-	addPeriodicData(elInfo);
-      elInfo = stack.traverseNext(elInfo);
-    }   
-
-    periodicDataCollected = true;
-  }
-
-
-  void DataCollector::addElementData(ElInfo* elInfo)
-  {
-    FUNCNAME("DataCollector::addElementData()");
-
-    const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
-
-    // create ElementInfo
-    ElementInfo elementInfo(dim);
-    
-    // read element region
-    ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
-    
-    if (ed)
-      elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
-    else 
-      elementInfo.elementRegion = -1;    
-   
-    // read surface regions to element info
-    ed = elInfo->getElement()->getElementData(SURFACE_REGION);
-    elementInfo.surfaceRegions.set(-1);
-    while (ed) {
-      SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
-      elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
-      ed = ed->getDecorated(SURFACE_REGION);
-    }
-
-    // for all vertices
-    for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
-      // get coords of this vertex
-      WorldVector<double> &vertexCoords = elInfo->getCoord(i);
-      
-      // get dof index of this vertex
-      DegreeOfFreedom vertexDof = dof[i][nPreDofs];
-
-      // search for coords at this dof
-      std::list<VertexInfo>::iterator it =
-	  find((*vertexInfos)[vertexDof].begin(), (*vertexInfos)[vertexDof].end(),
-	       vertexCoords);
-      
-      // coords not yet in list?
-      if (it == (*vertexInfos)[vertexDof].end()) {
-	// create new vertex info and increment number of vertices
-	VertexInfo newVertexInfo = {vertexCoords, nVertices};
-
-	// add new vertex info to list
-	(*vertexInfos)[vertexDof].push_front(newVertexInfo);
-
-	// set iterator to new vertex info
-	it = (*vertexInfos)[vertexDof].begin();	
-
-	nVertices++;
-      }
-      
-      // fill element info
-      elementInfo.vertexInfo[i] = it;
-      elementInfo.boundary[i] = elInfo->getBoundary(i);
-      elementInfo.projection[i] = elInfo->getProjection(i);
-      elementInfo.neighbour[i] = 
-	elInfo->getNeighbour(i) ?
-	outputIndices[elInfo->getNeighbour(i)->getIndex()] :
-	-1;
-    }
-    
-    if (dim == 3)
-      elementInfo.type = elInfo->getType();
-
-    // remember element info
-    elements.push_back(elementInfo);
-  }
-
-
-  void DataCollector::addValueData(ElInfo *elInfo)
-  {
-    FUNCNAME("DataCollector::addValueData()");
-    
-    WorldVector<double> vertexCoords;
-    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
-
-    // First, traverse all DOFs at the vertices of the element, determine
-    // their coordinates and add them to the corresponding entry in dofCoords.
-    for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
-      DegreeOfFreedom dofi = localDOFs[i];
-
-      (*interpPointInd)[dofi] = -2; // mark as vertex
-      
-      // get coords of this vertex
-      vertexCoords = elInfo->getCoord(i);
-
-      // search for coords at this dof
-      std::list<WorldVector<double> >::iterator it =
-	  find((*dofCoord)[dofi].begin(),
-	       (*dofCoord)[dofi].end(),
-	       vertexCoords);
-      
-      // coords not yet in list?
-      if (it == (*dofCoord)[dofi].end()) {
-	// add new coords to list
-	(*dofCoord)[dofi].push_back(vertexCoords);
-      }
-    }
-   
-    // Then, traverse all interpolation DOFs of the element, determine
-    // their coordinates and add them to the corresponding entry in 
-    // interpPointCoords.
-    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++) {
-      DegreeOfFreedom dofi = localDOFs[i];
-
-      elInfo->coordToWorld(*basisFcts->getCoords(i), vertexCoords);
-      
-      if ((*interpPointInd)[dofi] == -1) {
-	// mark as interpolation point
-	(*interpPointInd)[dofi] = -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)[dofi].begin(), 
-		     (*interpPointCoords)[dofi].end(),
-		     vertexCoords);
-	
-	if (it == (*interpPointCoords)[dofi].end()) {
-	  (*interpPointCoords)[dofi].push_back(vertexCoords); 
-	  nInterpPoints++;
-	}
-      }      
-    }  
-  }
-
-
-  void DataCollector::addInterpData(ElInfo *elInfo)
-  {
-    FUNCNAME("DataCollector::addInterpData()");
-    
-    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
-
-    std::vector<int> elemInterpPoints(0);
-    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++)
-      elemInterpPoints.push_back((*interpPointInd)[localDOFs[i]]);
-    
-    interpPoints.push_back(elemInterpPoints); 
-  }
-
-
-  void DataCollector::addPeriodicData(ElInfo *elInfo) 
-  {
-    FUNCNAME("DataCollector::addPeriodicData");
-
-    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
-      (elInfo->getElement()->
-       getElementData()->
-       getElementData(PERIODIC));
-    
-    if (ldp) {
-      std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
-      
-      for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
-	   it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
-	   ++it) {
-
-	int outputIndex = outputIndices[elInfo->getElement()->getIndex()];
-	int neighIndex  = outputIndices[elInfo->
-					 getNeighbour(it->elementSide)->
-					 getIndex()];
-	
-	if (!periodicConnections[outputIndex][it->elementSide]) {
-	  PeriodicInfo periodicInfo;
-	  
-	  periodicInfo.mode = it->periodicMode;
-	  periodicInfo.type = it->type;
-	  periodicInfo.outputIndex = outputIndex;
-	  periodicInfo.neighIndex = neighIndex;
-	  periodicInfo.vertexMap.clear();
-	  
-	  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),
-								   it->elementSide,
-								   i);
-	    int dof1 = elInfo->getElement()->getDof(index1, nPreDofs);
-	    
-	    for (int j = 0; j < dim; j++) {
-	      int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim),
-								     elInfo->getOppVertex(it->elementSide),
-								     j);
-	      int dof2 = elInfo->getNeighbour(it->elementSide)->getDof(index2, nPreDofs);
-	      
-	      if ((dof1 == dof2) || (mesh->associated(dof1, dof2))) {
-		periodicInfo.vertexMap[index1] = index2;
-		break;
-	      } 
-	    }
-	  }
-	  
-	  periodicInfos.push_back(periodicInfo);
-	}
-      }
-    }
-  }
-
-
-  std::list<ElementInfo>* DataCollector::getElementInfos()
-  {        
-    if (!elementDataCollected)
-      startCollectingElementData();
-      
-    return &elements;
-  }
-
-
-  DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
-  {
-    if (!elementDataCollected)
-      startCollectingElementData();
-
-    return vertexInfos;
-  }
-
-
-  std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
-  {
-    if (!periodicDataCollected)
-      startCollectingPeriodicData();
-
-    return &periodicInfos;
-  }
-
-
-  int DataCollector::getNumberVertices()
-  {
-    if (!elementDataCollected)
-      startCollectingElementData();
-
-    return nVertices;
-  }
-
-
-  int DataCollector::getNumberElements()
-  {
-    if (!elementDataCollected)
-      startCollectingElementData();
-
-    return nElements;
-  }
-
-
-  int DataCollector::getNumberInterpPoints()
-  {
-    if (!valueDataCollected)
-      startCollectingValueData();
-
-    return nInterpPoints;
-  }
-
-  
-  int DataCollector::getNumberConnections()
-  {
-    if (!periodicDataCollected)
-      startCollectingPeriodicData();
-
-    return nConnection;
-  }
-
-
-  DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
-  {
-    if (!valueDataCollected)
-      startCollectingValueData();
-
-    return dofCoord;
-  }
-
-
-  DOFVector<int>* DataCollector::getInterpPointInd()
-  {
-    if (!valueDataCollected)
-      startCollectingValueData();
-
-    return interpPointInd;
-  }