From 3abcc7558115bafee981d29e69709ca010768f2e Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 29 Sep 2013 14:58:31 +0000 Subject: [PATCH] Add a few comments. Fix a couple of oversights where the same function is used for both cell integrals and boundary integrals (though the evaluated function values are, in fact, only used for the boundary integrals). git-svn-id: https://svn.dealii.org/trunk@31011 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 98 +++++++++++++---------------- 1 file changed, 45 insertions(+), 53 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index b58953756a..2d5a3dd279 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -260,19 +260,19 @@ namespace Step42 stress_strain_tensor_linearized += stress_strain_tensor_kappa; } - //

Equation data: right hand side, boundary values, obstacles

+ //

Equation data: boundary forces, boundary values, obstacles

// // The following should be relatively standard. We need classes for - // the right hand side forcing term (which we here choose to be zero) + // the boundary forcing term (which we here choose to be zero) // and boundary values on those part of the boundary that are not part // of the contact surface (also chosen to be zero here). namespace EquationData { template - class RightHandSide : public Function + class BoundaryForce : public Function { public: - RightHandSide (); + BoundaryForce (); virtual double value (const Point &p, @@ -284,7 +284,7 @@ namespace Step42 }; template - RightHandSide::RightHandSide () + BoundaryForce::BoundaryForce () : Function(dim) {} @@ -292,7 +292,7 @@ namespace Step42 template double - RightHandSide::value (const Point &, + BoundaryForce::value (const Point &, const unsigned int) const { return 0.; @@ -300,11 +300,11 @@ namespace Step42 template void - RightHandSide::vector_value (const Point &p, + BoundaryForce::vector_value (const Point &p, Vector &values) const { for (unsigned int c = 0; c < this->n_components; ++c) - values(c) = RightHandSide::value(p, c); + values(c) = BoundaryForce::value(p, c); } @@ -356,7 +356,10 @@ namespace Step42 // at position $x=y=0.5, z=z_{\text{surface}}+0.59$ and radius $r=0.6$, // where $z_{\text{surface}}$ is the vertical position of the (flat) // surface of the deformable body. The function's value - // returns the location of the obstacle for a given $x,y$ value. + // returns the location of the obstacle for a given $x,y$ value if the + // point actually lies below the sphere, or a large positive value that + // can't possibly interfere with the deformation if it lies outside + // the "shadow" of the sphere. template class SphereObstacle : public Function { @@ -385,28 +388,28 @@ namespace Step42 template - double - SphereObstacle::value ( - const Point &p, const unsigned int component) const - { - if (component == 0) - return p(0); - else if (component == 1) - return p(1); - else if (component == 2) - { - double return_value = 1000.0; - if ((p(0) - 0.5) * (p(0) - 0.5) + (p(1) - 0.5) * (p(1) - 0.5) - < 0.36) - return_value = (-std::sqrt( - 0.36 - (p(0) - 0.5) * (p(0) - 0.5) + double + SphereObstacle::value ( + const Point &p, const unsigned int component) const + { + if (component == 0) + return p(0); + else if (component == 1) + return p(1); + else if (component == 2) + { + if ((p(0) - 0.5) * (p(0) - 0.5) + (p(1) - 0.5) * (p(1) - 0.5) + < 0.36) + return (-std::sqrt( + 0.36 - (p(0) - 0.5) * (p(0) - 0.5) - (p(1) - 0.5) * (p(1) - 0.5)) + z_surface + 0.59); - return return_value; - } + else + return 1000; + } - Assert(false, ExcNotImplemented()); - return 1e9; // an unreasonable value; ignored in debug mode because of the preceding Assert - } + Assert(false, ExcNotImplemented()); + return 1e9; // an unreasonable value; ignored in debug mode because of the preceding Assert + } template @@ -1005,11 +1008,9 @@ namespace Step42 const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values(n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face(n_face_q_points, - Vector(dim)); + const EquationData::BoundaryForce boundary_force; + std::vector > boundary_force_values(n_face_q_points, + Vector(dim)); FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); Vector cell_rhs(dofs_per_cell); @@ -1029,9 +1030,6 @@ namespace Step42 cell_matrix = 0; cell_rhs = 0; - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); - std::vector > strain_tensor(n_q_points); fe_values[displacement].get_function_symmetric_gradients(u, strain_tensor); @@ -1084,15 +1082,14 @@ namespace Step42 { fe_values_face.reinit(cell, face); - right_hand_side.vector_value_list( - fe_values_face.get_quadrature_points(), - right_hand_side_values_face); + boundary_force.vector_value_list(fe_values_face.get_quadrature_points(), + boundary_force_values); for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) { Tensor<1, dim> rhs_values; - rhs_values[2] = right_hand_side_values[q_point][2]; + rhs_values[2] = boundary_force_values[q_point][2]; for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values @@ -1106,7 +1103,7 @@ namespace Step42 local_dof_indices, system_matrix_newton, system_rhs_newton, true); - }; + } system_matrix_newton.compress(VectorOperation::add); system_rhs_newton.compress(VectorOperation::add); @@ -1133,11 +1130,9 @@ namespace Step42 const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); - const EquationData::RightHandSide right_hand_side; - std::vector > right_hand_side_values(n_q_points, - Vector(dim)); - std::vector > right_hand_side_values_face(n_face_q_points, - Vector(dim)); + const EquationData::BoundaryForce boundary_force; + std::vector > boundary_force_values(n_face_q_points, + Vector(dim)); Vector cell_rhs(dofs_per_cell); @@ -1160,9 +1155,6 @@ namespace Step42 fe_values.reinit(cell); cell_rhs = 0; - right_hand_side.vector_value_list(fe_values.get_quadrature_points(), - right_hand_side_values); - std::vector > strain_tensor(n_q_points); fe_values[displacement].get_function_symmetric_gradients(current_solution, strain_tensor); @@ -1205,14 +1197,14 @@ namespace Step42 { fe_values_face.reinit(cell, face); - right_hand_side.vector_value_list(fe_values_face.get_quadrature_points(), - right_hand_side_values_face); + boundary_force.vector_value_list(fe_values_face.get_quadrature_points(), + boundary_force_values); for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) { Tensor<1, dim> rhs_values; - rhs_values[2] = right_hand_side_values[q_point][2]; + rhs_values[2] = boundary_force_values[q_point][2]; for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values * fe_values_face.JxW(q_point)); -- 2.39.5