diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 4842a445d6c09b89d4cf36f51f2665ea8704d57d..aca2b01ff42281bd1ec60ffbcee5755489e47990 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -251,40 +251,43 @@ int main(int argc, char *argv[]) { MatrixType massMatrix; VectorType gravityFunctional; { - // Assemble mass matrix and compute density - MassAssembler<GridType, P1Basis::LocalFiniteElement, - P1Basis::LocalFiniteElement> const localMass; - OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis) - .assemble(localMass, massMatrix); + double const gravity = 9.81; + double const density = parset.get<double>("body.density"); - double volume = 1.0; - for (int i = 0; i < dim; ++i) - volume *= (upperRight[i] - lowerLeft[i]); + { // Assemble mass matrix + MassAssembler<GridType, P1Basis::LocalFiniteElement, + P1Basis::LocalFiniteElement> const localMass; + OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis) + .assemble(localMass, massMatrix); + massMatrix *= density; + } - double area = 1.0; - for (int i = 0; i < dim; ++i) - if (i != 1) - area *= (upperRight[i] - lowerLeft[i]); + { // Compute normal stress + double volume = 1.0; + for (int i = 0; i < dim; ++i) + volume *= (upperRight[i] - lowerLeft[i]); - double const gravity = 9.81; - double const density = parset.get<double>("body.density"); + double area = 1.0; + for (int i = 0; i < dim; ++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; + } - // 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 - SmallVector weightedGravitationalDirection(0); - weightedGravitationalDirection[1] = -density * gravity; - ConstantFunction<SmallVector, SmallVector> const gravityFunction( - weightedGravitationalDirection); - L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler( - gravityFunction); - FunctionalAssembler<P1Basis>(p1Basis) - .assemble(gravityFunctionalAssembler, gravityFunctional, true); + { // Compute gravitational body force + SmallVector weightedGravitationalDirection(0); + weightedGravitationalDirection[1] = -density * gravity; + ConstantFunction<SmallVector, SmallVector> const gravityFunction( + weightedGravitationalDirection); + L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler( + gravityFunction); + FunctionalAssembler<P1Basis>(p1Basis) + .assemble(gravityFunctionalAssembler, gravityFunctional, true); + } } // Assemble elastic force on the body