more debug printing

dune/tectonic/utils/debugutils.hh

@@ -14,6 +14,42 @@
 
 //using namespace std;
 
+    template <typename VectorType>
+    auto average(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType res(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                res += vec[i];
+            }
+            res *= 1.0/len;
+        }
+
+        return res;
+    }
+
+    template <typename VectorType>
+    auto minmax(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType min(0.0), max(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                if (vec[i] > max) {
+                    max = vec[i];
+                }
+                if (vec[i] < min) {
+                    min = vec[i];
+                }
+            }
+        }
+        return std::array{min, max};
+    }
 
     template <int s>
     void print(const Dune::BitSetVector<s>& x, std::string message){

src/foam/foam.cc

@@ -194,7 +194,7 @@ int main(int argc, char *argv[]) {
     for (size_t i=0; i<bodyCount; i++) {
         nVertices[i] = contactNetwork.body(i)->nVertices();
     }
-
+    print(nVertices, "#dofs: ");
 
     MyProgramState programState(nVertices);
 
@@ -381,7 +381,7 @@ int main(int argc, char *argv[]) {
     // ------------------------
     // assemble global friction
     // ------------------------
-    print(programState.weightedNormalStress, "weightedNormalStress");
+    //print(programState.weightedNormalStress, "weightedNormalStress");
 
     contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
 
@@ -491,6 +491,7 @@ int main(int argc, char *argv[]) {
       current.rate_->extractVelocity(programState.v);
       current.rate_->extractAcceleration(programState.a);
       current.state_->extractAlpha(programState.alpha);
+
       globalFriction.updateAlpha(programState.alpha);
 
       /*print(programState.u, "current u:");