More structure

src/one-body-sample.cc

@@ -251,12 +251,18 @@ int main(int argc, char *argv[]) {
     MatrixType massMatrix;
     VectorType gravityFunctional;
     {
-      // Assemble mass matrix and compute density
+      double const gravity = 9.81;
+      double const density = parset.get<double>("body.density");
+
+      { // Assemble mass matrix
         MassAssembler<GridType, P1Basis::LocalFiniteElement,
                       P1Basis::LocalFiniteElement> const localMass;
         OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
             .assemble(localMass, massMatrix);
+        massMatrix *= density;
+      }
 
+      { // Compute normal stress
         double volume = 1.0;
         for (int i = 0; i < dim; ++i)
           volume *= (upperRight[i] - lowerLeft[i]);
@@ -266,17 +272,13 @@ int main(int argc, char *argv[]) {
           if (i != 1)
             area *= (upperRight[i] - lowerLeft[i]);
 
-      double const gravity = 9.81;
-      double const density = parset.get<double>("body.density");
-
         // volume * gravity * density / area    = normal stress
         // V      * g       * rho     / A       = sigma_n
         // m^d    * N/kg    * kg/m^d  / m^(d-1) = N/m^(d-1)
         normalStress = volume * gravity * density / area;
+      }
 
-      massMatrix *= density;
-
-      // Compute gravitational body force
+      { // Compute gravitational body force
         SmallVector weightedGravitationalDirection(0);
         weightedGravitationalDirection[1] = -density * gravity;
         ConstantFunction<SmallVector, SmallVector> const gravityFunction(
@@ -286,6 +288,7 @@ int main(int argc, char *argv[]) {
         FunctionalAssembler<P1Basis>(p1Basis)
             .assemble(gravityFunctionalAssembler, gravityFunctional, true);
       }
+    }
 
     // Assemble elastic force on the body
     MatrixType stiffnessMatrix;