diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index da3020b3c1731e7e80dde5fd30cf1a0b98cc3741..ef9035948a370f6c7cf66b94a4bf727b9bd0df78 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -178,26 +178,6 @@ int main(int argc, char *argv[]) { functions.get("velocityDirichletCondition"), &neumannFunction = functions.get("neumannCondition"); - // Set up normal stress, mass matrix, and gravity functional - double normalStress; - { - double volume = 1.0; - for (size_t i = 0; i < dims; ++i) - volume *= (upperRight[i] - lowerLeft[i]); - - double area = 1.0; - for (size_t i = 0; i < dims; ++i) - if (i != 1) - area *= (upperRight[i] - lowerLeft[i]); - - // 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; - } - FrictionData const frictionData(parset.sub("boundary.friction"), - normalStress); - using MyAssembler = MyAssembler<GridView, dims>; using Matrix = MyAssembler::Matrix; using LocalMatrix = Matrix::block_type; @@ -235,6 +215,25 @@ int main(int argc, char *argv[]) { Vector ell(fineVertexCount); computeExternalForces(0.0, ell); + double normalStress; + { + double volume = 1.0; + for (size_t i = 0; i < dims; ++i) + volume *= (upperRight[i] - lowerLeft[i]); + + double area = 1.0; + for (size_t i = 0; i < dims; ++i) + if (i != 1) + area *= (upperRight[i] - lowerLeft[i]); + + // 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; + } + FrictionData const frictionData(parset.sub("boundary.friction"), + normalStress); + // {{{ Initial conditions ScalarVector alpha_initial(fineVertexCount); alpha_initial =