sonstiges

extensions/GradientCalculations.h

+/** \file GradientCalculations.h */
+
+#ifndef GRADIENT_CALCULATIONS_H
+#define GRADIENT_CALCULATIONS_H
+
+#include "AMDiS.h"
+
+using namespace AMDiS;
+
+namespace experimental {
+
+  inline void getRecoveryGradientSym(const DOFVector<double> *u,
+				     DOFVector<WorldVector<double> > *result, 
+				     bool asym = true)
+  {
+    FUNCNAME("getRecoveryGradientSym()");
+    const FiniteElemSpace *feSpace = u->getFeSpace();
+    int dim = feSpace->getMesh()->getDim();
+
+    if (!result)
+      result= new DOFVector<WorldVector<double> >(feSpace, "grad");
+
+    WorldVector<double> grd;
+    nullify(grd);
+
+    result->set(grd);
+
+    DOFVector<double> volume(feSpace, "volume");
+    volume.set(0.0);
+
+    DOFVector<WorldVector<double> > nrm(feSpace, "normal");
+
+    const BasisFunction *basFcts = feSpace->getBasisFcts();
+    int nBasisFcts = basFcts->getNumber();
+    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
+
+    Mesh *mesh = feSpace->getMesh();
+    Flag flags(Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_COORDS |
+	      Mesh::FILL_GRD_LAMBDA | Mesh::FILL_BOUND);
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, flags);
+
+    mtl::dense_vector<double> localUh(nBasisFcts);
+    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
+    std::map<DegreeOfFreedom, bool> visited;
+
+    while (elInfo) {
+      double det = elInfo->getDet();
+      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+      u->getLocalVector(elInfo->getElement(), localUh);
+      basFcts->evalGrdUh(bary, grdLambda, localUh, grd);
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+
+      for (int i = 0; i < (dim + 1); i++) {
+	  (*result)[localIndices[i]] += grd * det;
+	  volume[localIndices[i]] += det;
+
+	  if (elInfo->getBoundary(i)) {
+	    WorldVector<double> n, rgrd;
+	    elInfo->getNormal(i, n);
+
+	    if (asym)
+	      rgrd = grd;
+	    else
+	      rgrd = grd - 2.0*n*(grd*n); // reflect gradient at boundary
+
+	    for (int j = 0; j < (dim + 1); j++)
+	      if (j != i) {
+		(*result)[localIndices[j]] += rgrd * det; // elements with boundary face
+		volume[localIndices[j]] += det;
+
+		nrm[localIndices[j]] += n;        // normal for other elements sharing boundary nodes
+		visited[localIndices[j]] = true;
+	      }
+	  }
+      }
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    elInfo = stack.traverseFirst(mesh, -1, flags);
+    while (elInfo) {
+      double det = elInfo->getDet();
+      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+
+      for (int i = 0; i < (dim + 1); i++) {
+	if (visited[localIndices[i]]) {
+	  int j = 0;
+	  for (; j < (dim + 1); j++)
+	    if (j != i && elInfo->getBoundary(j))
+	      break;
+
+	  if (j == (dim + 1)) {
+	    WorldVector<double> rgrd, n;
+	    u->getLocalVector(elInfo->getElement(), localUh);
+	    basFcts->evalGrdUh(bary, grdLambda, localUh, grd);
+
+	    if (asym)
+	      rgrd= grd;
+	    else {
+	      n = 1.0/norm(nrm[localIndices[i]]) * nrm[localIndices[i]];
+	      rgrd = grd - 2.0*n*(grd*n);
+	    }
+
+	    (*result)[localIndices[i]] += rgrd * det;    // elements with boundary node but no face
+	    volume[localIndices[i]] += det;
+	  }
+	}
+      }
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    DOFIterator<double> volIt(&volume, USED_DOFS);
+    DOFIterator<WorldVector<double> > grdIt(result, USED_DOFS);
+
+    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
+      if (*volIt != 0.0)
+	*grdIt *= 1.0 / (*volIt);
+  }
+
+} // end namespace experimental
+
+#endif // GRADIENT_CALCULATIONS_H
\ No newline at end of file

extensions/MetaTools.h

@@ -19,10 +19,11 @@
 // }
 
 namespace tools {
-    
-  /// binomial coefficients
-  /// _______________________________________________________________________________
   
+  /**
+    * \defgroup binomial_coefficients binomial coefficients
+    * @{
+    */
   namespace details {
 
     template <int N, int K, bool N_pos = true, bool K_pos = true>
@@ -51,54 +52,74 @@ namespace tools {
     };
   }
   
-  // interface to calc binomial coefficients (N, K) at compile-time for N,K >= 0
+  /// \brief
+  /// interface to calc binomial coefficients (N, K) at compile-time for N,K >= 0
+  ///
   template <int N, int K>
   struct Binomial
   {
     BOOST_STATIC_CONSTANT(int, value = (details::Binomial<N, K, (N>0), (K>0)>::value) );
   };
+  /**@}*/
 
-  /// some mathematical basics
-  /// ____________________________________________________________________________
+// ______________________________________________________________________________________________________
+  
+  /**
+    * \defgroup math_basics some mathematical basics
+    * @{
+    */
   
-  // check divisibility
+  /// \brief
+  /// check divisibility
+  ///
   template <int u, int v>
   struct Divisible
   {
     BOOST_STATIC_CONSTANT(bool, value = u % v == 0 ? true : false );
   };
     
-  // check divide by zero condition
+  /// \brief
+  /// check divide by zero condition
+  ///
   template <int u>
   struct Divisible<u, 0>
   {
     BOOST_STATIC_CONSTANT(bool, value = false ); // better: static assert
   };
   
-  // check number is even or not
+  /// \brief
+  /// check number is even or not
+  ///
   template <int u>
   struct IsEven
   {
     BOOST_STATIC_CONSTANT(bool, value = (Divisible<u, 2>::value) );
   };
     
-  // check number is odd or not
+  /// \brief
+  /// check number is odd or not
+  ///
   template <int u>
   struct IsOdd
   {
     BOOST_STATIC_CONSTANT(bool, value = (-Divisible<u, 2>::value) );
   };
   
+  /// \brief
+  /// absolute value calculation
+  ///
   template <int N>
   struct Abs
   {
     BOOST_STATIC_CONSTANT(bool, value = (N >= 0) ? N : -N );
   };
-  
-  // greatest common devisor and least common multiple
-  // ________________________________________________________________________
 
-  // calculate gcd
+// ______________________________________________________________________________________________________
+//  greatest common devisor and least common multiple
+
+  /// \brief
+  /// calculate gcd of template parameters u and v
+  ///
   template <int u, int v>
   struct Gcd
   {
@@ -117,18 +138,21 @@ namespace tools {
     BOOST_STATIC_CONSTANT(int, value = -1 );
   };
   
-  // calculate lcm
+  /// \brief
+  /// calculate the least common multiple of template parameters u and v
+  ///
   template <int u, int v>
   struct Lcm
   {
     BOOST_STATIC_CONSTANT(int, value = (Abs<u * v>::value / Gcd<u, v>::value) );
   };
-  
-  
-  // powers of values
-  // _______________________________________________________________________________
-  
-  // to calculate the power
+
+// ______________________________________________________________________________________________________
+// powers of values
+ 
+  /// \brief
+  /// calculate the power of template parameters a<sup>p</sup>
+  ///
   template <int a, int p>
   struct Pow
   {
@@ -152,11 +176,13 @@ namespace tools {
   {
     BOOST_STATIC_CONSTANT(int, value = (Pow<n, 2>::value) );
   };
-  
-  // roots of values
-  // _______________________________________________________________________________
-  
-  // template to compute sqrt(N) via iteration
+
+// ______________________________________________________________________________________________________
+// roots of values
+ 
+  /// \brief
+  /// template to compute sqrt(a) via iteration
+  ///
   template <int a, int I=1>
   struct Sqrt {
       // instantiate next step or result type as branch
@@ -179,7 +205,9 @@ namespace tools {
     BOOST_STATIC_CONSTANT(int, value = 0 );
   };
   
-  // template to compute sqrt(N) via iteration
+  /// \brief
+  /// template to compute the Nth root of a via iteration
+  ///
   template <int a, int N, int I=1>
   struct Root {
     // instantiate next step or result type as branch
@@ -201,10 +229,10 @@ namespace tools {
   struct Root<0, N, I> {
     BOOST_STATIC_CONSTANT(int, value = 0 );
   };
-  
-  
-  /// monomials and polynomials
-  /// _______________________________________________________________________________
+  /**@}*/
+
+// ______________________________________________________________________________________________________
+// monomials and polynomials
 
   template<int N>
   struct Monomial {
@@ -256,10 +284,10 @@ namespace tools {
   private:
     Coefficients& coefficients;
   };
-  
-  /// Bernstein polynomials and Bezier-Surfaces
-  /// ____________________________________________________________________________
-  
+
+// ______________________________________________________________________________________________________
+// Bernstein polynomials and Bezier-Surfaces
+
   template<int I, int order>
   struct BernsteinPolynomial 
   {
@@ -309,9 +337,9 @@ namespace tools {
       return BernsteinPolynomial<I, order>::value(x) * BernsteinPolynomial<J, order>::value(x) * BernsteinPolynomial<K, order>::value(z);
     }
   };
-  
-  // Bezier-Polynomial
-  //________________________________________________________________________
+
+// ______________________________________________________________________________________________________
+// Bezier-Polynomial
   
   namespace details {
     
@@ -350,8 +378,8 @@ namespace tools {
     Coefficients& coefficients;
   };
 
-  // Bezier-Polynomial 2d
-  //________________________________________________________________________
+// ______________________________________________________________________________________________________
+
   
   namespace details {
     
@@ -375,7 +403,9 @@ namespace tools {
     };
   }
     
-  // interface for bezier-surfaces of order (N,N)
+  /// \brief
+  /// interface for bezier-surfaces of order (N,N)
+  ///
   template<typename Coefficients>
   struct BezierPolynomial2
   {
@@ -396,10 +426,9 @@ namespace tools {
   private:
     Coefficients& coefficients; // p in points is element of result_type
   };
-  
-  
-  /// for loops: for<i, n [,{incr,decr}]>, for<mpl::range_c>, for<mpl::vector_c>
-  ///____________________________________________________________________________
+
+// ______________________________________________________________________________________________________
+// for loops: for<i, n [,{incr,decr}]>, for<mpl::range_c>, for<mpl::vector_c>
 
   template<int I> struct incr { BOOST_STATIC_CONSTANT(int, value = I+1); };
   template<int I> struct decr { BOOST_STATIC_CONSTANT(int, value = I-1); };
@@ -433,15 +462,17 @@ namespace tools {
     static void loop() { FOR_LOOP<I, N, decr, finished, F>::loop(); }
   };
 
-  // interface for classical for loop: for j=I..N
-  // the index will be incremented, if I<N and decremented otherwise
+  /// \brief
+  /// interface for classical for loop: for j=I..N
+  /// the index will be incremented, if I&lt;N and decremented otherwise
+  ///
   template<int I, int N, template<int> class F> 
   struct FOR 
   { 
     static void loop() { FOR_LOOP_CHECK<I, N, I < N, I == N, F>::loop(); } 
   };
 
-  // 'for-loop' over sequence of values
+  // for-loop over sequence of values
   template<typename Seq, long I, long N, template<int> class F>
   struct FOR_SUBSEQ
   {
@@ -459,17 +490,18 @@ namespace tools {
     static void loop() {} 
   };
 
-  // interface for loop over sequence of values
-  // for i in {seq[0]..seq[end]}
+  /// \brief
+  /// interface for loop over sequence of values
+  /// for i in {seq[0]..seq[end]}
+  ///
   template<typename Seq, template<int> class F> 
   struct FOR_SEQ 
   { 
     static void loop() { FOR_SUBSEQ<Seq, 0, boost::mpl::size<Seq>::value, F>::loop(); } 
   };
 
-
-  /// for-each loops
-  ///____________________________________________________________________________
+//____________________________________________________________________________
+// for-each loops
 
   template<int I> struct for_each
   {

extensions/VtuReader.hh

@@ -16,7 +16,6 @@ namespace AMDiS {
   namespace VtuReader {
     
     /** \brief
-     * Copies the values of a VTU file to a DOF vector.
      *
      * The DOF vector must have been created by a corresponding mesh file. The
      * function now reads the corresponding VTU file and
@@ -50,7 +49,6 @@ namespace AMDiS {
       TEST_EXIT(dofVectors.size() > 0)("no DOF vectors specified\n");
       TEST_EXIT(componentNames.size() > 0)("no componentName specified\n");
 
-
       if(!boost::filesystem::exists(filename))
 	throw(std::runtime_error(filename + " does not exist!"));
 
@@ -80,9 +78,9 @@ namespace AMDiS {
 	    int nComponents = -1;
 	    if (DataArray.attribute("NumberOfComponents"))
 	      nComponents = DataArray.attribute("NumberOfComponents").as_int();
-	    if (nComponents != -1 && num_rows(test) > nComponents)
+	    if (nComponents != -1 && vector_operations::num_rows(test) > nComponents)
 	      throw(std::runtime_error("Can not store values in DOFVector with given value type. Too many components!"));
-	    string2valueList(values, valueList[i], num_rows(test), nComponents);
+	    string2valueList(values, valueList[i], vector_operations::num_rows(test), nComponents);
 	    break;
 	  }
 	}