Prescribe density rather than normal stress

src/one-body-sample.cc

@@ -179,8 +179,6 @@ int main(int argc, char *argv[]) {
       OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
           .assemble(localMass, massMatrix);
 
-      // We have volume*gravity*density/area = normalstress (V*g*rho/A =
-      // sigma_n)
       double volume = 1.0;
       for (int i = 0; i < dim; ++i)
         volume *= (upperRight[i] - lowerLeft[i]);
@@ -191,11 +189,12 @@ int main(int argc, char *argv[]) {
           area *= (upperRight[i] - lowerLeft[i]);
 
       double const gravity = 9.81;
-      normalStress = parset.get<double>("boundary.friction.normalstress");
+      double const density = parset.get<double>("body.density");
 
-      // rho    = sigma     * A       / V   / g
-      // kg/m^d = N/m^(d-1) * m^(d-1) / m^d / (N/kg)
-      double const density = normalStress * area / volume / gravity;
+      // 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;
       if (parset.get<bool>("enableTimer"))

src/one-body-sample.parset

@@ -21,6 +21,7 @@ refinements = 4
 [body]
 E = 5e7
 nu = 0.3 # The closer we get to 0.5, the more wiggly everything gets
+density = 5000
 
 [solver]
 tolerance = 1e-10
@@ -49,7 +50,6 @@ acceptFactor = 1.0
 requiredResidual = 1e-12
 
 [boundary.friction]
-normalstress = 50000
 # "mu_0 is the nominal coefficient of friction that has values near 0.6"
 # -- James H. Dieterich and Brian Kilgore: Implications of fault
 #    constitutive properties for earthquake prediction