FileWriter in parallel corrected

AMDiS/src/CouplingProblemStat.h

@@ -69,7 +69,7 @@ namespace AMDiS {
       virtual void initialize(Flag initFlag,
 			      ProblemStatSeq *adoptProblem = nullptr,
 			      Flag adoptFlag = INIT_NOTHING)
-      {
+      { FUNCNAME("CouplingProblemStat::initialize()");
 
 	// create one refinement-/coarseningmanager for all problems
 	if (refinementManager != nullptr && coarseningManager != nullptr) { 
@@ -118,7 +118,6 @@ namespace AMDiS {
 	      meshesForProblems[i].insert(meshByName[meshName]);
 
 	    problems[i]->setComponentMesh(j, meshByName[meshName]);
-		    MSG((problems[i]->getName() + "->setComponentMesh(" + boost::lexical_cast<std::string>(j) + ") = " + meshName + "\n").c_str());
 	      
 	    // feSpace
 	    int degree = 1;
@@ -172,7 +171,7 @@ namespace AMDiS {
 
       void createRefCoarseManager() 
       {
-	FUNCNAME("ProblemStat::createRefCoarseManager()");
+	FUNCNAME("CouplingProblemStat::createRefCoarseManager()");
     
 	int dim = 0;
 	Parameters::get(name + "->dim", dim);

AMDiS/src/FixVec.h

@@ -660,6 +660,12 @@ namespace AMDiS {
   {
     typedef typename ProductType<T1,T2>::type type;
   };
+  
+  
+  
+  template<typename T,GeoIndex d>
+  inline void resize(FixVec<T,d>& vec, size_t newSize)
+  { }
 }
 
 

AMDiS/src/GenericOperatorTerm.h

@@ -271,7 +271,7 @@ struct GenericFirstOrderTerm_b : public FirstOrderTerm
     const int nPoints = static_cast<int>(Lb.size());
 
     for (int iq = 0; iq < nPoints; iq++)
-      lb(grdLambda, term(iq), Lb[iq], fac);
+      lb(grdLambda, term(iq), Lb[iq], 1.0);
   }
 };
 

AMDiS/src/MatrixVector.h

@@ -531,7 +531,20 @@ namespace AMDiS {
     z[0] = x[1] * y[2] - x[2] * y[1];
     z[1] = x[2] * y[0] - x[0] * y[2];
     z[2] = x[0] * y[1] - x[1] * y[0];
+  }  
+  
+  template<typename T>
+  inline size_t num_rows(const Vector<T>& vec)
+  {
+    return static_cast<size_t>(vec.getSize());
+  }
+  
+  template<typename T>
+  inline void resize(Vector<T>& vec, size_t newSize)
+  {
+    vec.resize(newSize);
   }
+  
 }
 
 #endif // AMDIS_MATRIXVECTOR_H

AMDiS/src/RobinBC.cc

@@ -171,7 +171,6 @@ namespace AMDiS {
       robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
       for (int i = 0; i < dim + 1; i++) 
 	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
-      WARNING("Sign of alpha changed in RobinBC!\n");
     }
 
   }
@@ -188,10 +187,11 @@ namespace AMDiS {
 
     int dim = elInfo->getMesh()->getDim();
 
-    if (neumannOperators)
+    if (neumannOperators) {
       for (int i = 0; i < dim + 1; i++)
 	if (elInfo->getBoundary(i) == boundaryType)
 	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
+    }
   }
 
 

AMDiS/src/Triangle.cc

@@ -299,7 +299,7 @@ namespace AMDiS {
 	child[1]->getSubBoundary(nextBound1, subBound);
 	child[0]->getSubBoundary(nextBound0, subBound);
       } else {
-	child[0]->getSubBoundary(nextBound0, subBound);
+	child[0]->getSubBoundary(nextBound0, subBound); // TODO: check this!
 	child[1]->getSubBoundary(nextBound1, subBound);
       }
     } else {

AMDiS/src/expressions/AtomicExpression.h

@@ -38,6 +38,7 @@ namespace AMDiS {
   
 namespace result_of {
       
+  /// initialize a coordinate vector at quadrature points
   struct Coords : public LazyOperatorTermBase
   {
     typedef WorldVector<double> value_type;
@@ -69,8 +70,6 @@ namespace result_of {
     }
 
     inline value_type operator()(const int& iq) const { return x[iq]; }
-    
-    inline value_type eval(const int& iq) const { return x[iq]; }
     inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
   };
   
@@ -106,12 +105,157 @@ namespace result_of {
     }
 
     inline double operator()(const int& iq) const { return x[iq][I]; }
-    
-    inline double eval(const int& iq) const { return x[iq][I]; }
     inline double derivative(const int& iq, int identifier) const { return 0.0; }
   };
   
   
+  /// initialize a normal vector at quadrature points
+  struct Normals : public LazyOperatorTermBase
+  {
+    typedef WorldVector<double> value_type;
+    mutable WorldVector<double> normal;
+    int side;
+
+    Normals(int side_) : side(side_) {}
+
+    template<typename List>
+    void insertFeSpaces(List& feSpaces) const {}
+    
+    int getDegree() const
+    {
+      return 1;
+    }
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* elInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      elInfo->getNormal(side, normal);
+    }
+
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      smallElInfo->getNormal(side, normal);
+    }
+
+    inline value_type operator()(const int& iq) const { return normal; }
+    inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
+  };
+  
+  
+  struct Normal : public LazyOperatorTermBase
+  {
+    typedef double value_type;
+    mutable WorldVector<double> normal;
+    int side;
+    int I;
+
+    Normal(int side_, int I_) : side(side_), I(I_) {}
+
+    template<typename List>
+    void insertFeSpaces(List& feSpaces) const {}
+    
+    int getDegree() const
+    {
+      return 1;
+    }
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* elInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      elInfo->getNormal(side, normal);
+    }
+
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      smallElInfo->getNormal(side, normal);
+    }
+
+    inline value_type operator()(const int& iq) const { return normal[I]; }
+    inline value_type derivative(const int& iq, int identifier) const { 0.0; }
+  };
+  
+  
+  
+  /// initialize a normal vector at quadrature points
+  struct ElementNormals : public LazyOperatorTermBase
+  {
+    typedef WorldVector<double> value_type;
+    mutable WorldVector<double> elementNormal;
+
+    ElementNormals() {}
+
+    template<typename List>
+    void insertFeSpaces(List& feSpaces) const {}
+    
+    int getDegree() const
+    {
+      return 1;
+    }
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* elInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      elInfo->getElementNormal(elementNormal);
+    }
+
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      smallElInfo->getElementNormal(elementNormal);
+    }
+
+    inline value_type operator()(const int& iq) const { return elementNormal; }
+    inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
+  };
+  
+  
+  struct ElementNormal : public LazyOperatorTermBase
+  {
+    typedef double value_type;
+    mutable WorldVector<double> elementNormal;
+    int I;
+
+    ElementNormal(int I_) : I(I_) {}
+
+    template<typename List>
+    void insertFeSpaces(List& feSpaces) const {}
+    
+    int getDegree() const
+    {
+      return 1;
+    }
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* elInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      elInfo->getElementNormal(elementNormal);
+    }
+
+
+    template<typename OT>
+    void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      smallElInfo->getElementNormal(elementNormal);
+    }
+
+    inline value_type operator()(const int& iq) const { return elementNormal[I]; }
+    inline value_type derivative(const int& iq, int identifier) const { 0.0; }
+  };
+  
+  
   template<typename T>
   struct Const : public LazyOperatorTermBase
   {
@@ -136,8 +280,6 @@ namespace result_of {
 	    SubAssembler* subAssembler, Quadrature *quad) {}
 
     inline value_type operator()(const int& iq) const { return value; }
-    
-    inline value_type eval(const int& iq) const { return value; }
     inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
   };
   
@@ -167,10 +309,15 @@ namespace result_of {
 }
 
 result_of::Coords X() { return result_of::Coords(); }
-
 template<int I>
 result_of::Coord<I> X() { return result_of::Coord<I>(); }
 
+result_of::Normals N(int side) { return result_of::Normals(side); }
+result_of::Normal N(int side, int I) { return result_of::Normal(side, I); }
+
+result_of::ElementNormals M() { return result_of::ElementNormals(); }
+result_of::ElementNormal M(int I) { return result_of::ElementNormal(I); }
+
 template<typename T>
 result_of::Const<T> constant(const T& value) { return result_of::Const<T>(value); }
 

AMDiS/src/expressions/LazyExpression.h

@@ -572,7 +572,7 @@ template<typename Term>
 inline typename boost::enable_if<
   typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
   result_of::Abs<Term> >::type
-abs(const Term& t) { return result_of::Abs<Term>(t); }
+abs_(const Term& t) { return result_of::Abs<Term>(t); } // TODO: Funktionsnamen ohne Unterstrich
 
 
 // signum of a term

AMDiS/src/expressions/ValueOf.h

@@ -82,7 +82,7 @@ namespace result_of {
   };
   
   
-  template<typename T, template<class> class Vector>
+  template<template<class> class Vector, typename T>
   struct ValueOf<Vector<DOFVector<T>*> > : public LazyOperatorTermBase
   {
     typedef Vector<T> value_type;
@@ -108,12 +108,12 @@ namespace result_of {
     inline void initElement(OT& ot, const ElInfo* elInfo,
 	    SubAssembler* subAssembler, Quadrature *quad)
     {
-      Vector<mtl::dense_vector<T> > helper(num_rows(vecDV));
+      Vector<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV));
       for (size_t i = 0; i < num_rows(vecDV); i++)
 	ot.getVectorAtQPs(vecDV[i], elInfo, subAssembler, quad, helper[i]);
       vec.change_dim(num_rows(helper[0]));
       for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
-	value_type tmp(num_rows(vecDV));
+	value_type tmp; resize(tmp, num_rows(vecDV));
 	for (size_t i = 0; i < num_rows(vecDV); i++)
 	  tmp[i] = helper[i][iq];
 	vec[iq] = tmp;
@@ -125,12 +125,12 @@ namespace result_of {
     inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
 	    SubAssembler* subAssembler, Quadrature *quad)
     {      
-      Vector<mtl::dense_vector<T> > helper(num_rows(vecDV));
+      Vector<mtl::dense_vector<T> > helper; resize(helper, num_rows(vecDV));
       for (size_t i = 0; i < num_rows(vecDV); i++)
 	ot.getVectorAtQPs(vecDV[i], smallElInfo, largeElInfo, subAssembler, quad, helper[i]);
       vec.change_dim(num_rows(helper[0]));
       for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
-	value_type tmp(num_rows(vecDV));
+	value_type tmp; resize(tmp, num_rows(vecDV));
 	for (size_t i = 0; i < num_rows(vecDV); i++)
 	  tmp[i] = helper[i][iq];
 	vec[iq] = tmp;
@@ -193,6 +193,51 @@ namespace result_of {
       return static_cast<value_type>(identifier == I);
     }
   };
+  
+  template<typename Vector>
+  struct ComponentOf : public LazyOperatorTermBase {};
+  
+  template<template<class> class Vector, typename T>
+  struct ComponentOf<DOFVector<Vector<T> > > : public LazyOperatorTermBase
+  {
+    typedef T value_type;
+
+    DOFVector<Vector<T> >* vecDV;
+    mutable mtl::dense_vector<Vector<T> > vec;
+    int I;
+
+    ComponentOf(DOFVector<Vector<T> >& vector, int I_) : vecDV(&vector), I(I_) {}
+    ComponentOf(DOFVector<Vector<T> >* vector, int I_) : vecDV(vector), I(I_) {}
+
+    template<typename List>
+    inline void insertFeSpaces(List& feSpaces) const
+    {
+      feSpaces.insert(vecDV->getFeSpace());
+    }
+    
+    inline int getDegree() const
+    {
+      return vecDV->getFeSpace()->getBasisFcts()->getDegree();
+    }
+
+    template<typename OT>
+    inline void initElement(OT& ot, const ElInfo* elInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+    }
+
+
+    template<typename OT>
+    inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler, Quadrature *quad)
+    {
+      ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+    }
+
+    inline value_type operator()(const int& iq) const { return vec[iq][I]; }
+    inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
+  };
 
 
 } // end namespace result_of
@@ -204,6 +249,15 @@ result_of::ValueOf<DOFVector<T> > valueOf(DOFVector<T>& vector) { return result_
 template<typename T>
 result_of::ValueOf<DOFVector<T> > valueOf(DOFVector<T>* vector) { return result_of::ValueOf<DOFVector<T> >(vector); }
 
+template<template<class> class Vector, typename T>
+result_of::ValueOf<Vector<DOFVector<T>*> > valueOf(Vector<DOFVector<T>*> &vector) { return result_of::ValueOf<Vector<DOFVector<T>*> >(vector); }
+
+template<template<class> class Vector, typename T>
+result_of::ComponentOf<DOFVector<Vector<T> > > componentOf(DOFVector<Vector<T> >& vector, int I) { return result_of::ComponentOf<DOFVector<Vector<T> > >(vector, I); }
+
+template<template<class> class Vector, typename T>
+result_of::ComponentOf<DOFVector<Vector<T> > > componentOf(DOFVector<Vector<T> >* vector, int I) { return result_of::ComponentOf<DOFVector<Vector<T> > >(vector, I); }
+
 } // end namespace AMDiS
 
 #endif // AMDIS_VALUE_OF_H
\ No newline at end of file

AMDiS/src/io/FileWriter.cc

@@ -127,9 +127,11 @@ namespace AMDiS
 	path vtu_path = fn;
 	path vtu_filename = vtu_path.filename();
 	vtu_path.remove_filename() /= vtu_filename.stem();
-	create_directory(vtu_path); 
-	vtu_path /= vtu_filename;	  
-	fn = vtu_path.string();
+	try {
+	  create_directory(vtu_path); 
+	  vtu_path /= vtu_filename;	  
+	  fn = vtu_path.string();
+	} catch (...) {}
       }
 
   #if HAVE_PARALLEL_DOMAIN_AMDIS

AMDiS/src/io/FileWriter.hh

@@ -130,6 +130,7 @@ namespace AMDiS
       paraViewPrecision = 0;
       writeParaViewVectorFormat = 0;
       writeAs3dVector = false;
+      createParaViewSubDir = false;
       writeParaViewAnimation = 0;
       writePeriodicFormat = 0;
       writePngFormat = 0;

AMDiS/src/io/detail/VtkWriter.cc

@@ -250,7 +250,9 @@ namespace AMDiS  { namespace io {
 	      << "    </PCells>\n";
 	  file << "    <PPointData>\n";
 	  
-	  for (size_t i = 0; i < (writeAsVector ? 1 : componentNames.size()); i++) { 
+	  int nValues = static_cast<int>(componentNames.size());
+	  nValues = writeAsVector ? std::min(nValues,1) : nValues;	
+	  for (size_t i = 0; i < nValues; i++) {
 	    if(highPrecision && format != ::AMDiS::io::VtkWriter::ASCII)
 	      file << "      <PDataArray type=\"Float64\" Name=\"";
 	    else