RE: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal

[Michael Li]

Alex,   Thank you for sharing your codes. I have some compiling errors relating to “EnergyFunctional’: solve_ring_nonlinear.cpp:520:43: error: ‘EnergyFunctional’ in namespace ‘dealii::Differentiation::AD’ does not name a template type   520 |     using ADHelper = Differentiation::AD::EnergyFunctional<       |                                           ^~~~~~~~~~~~~~~~ Can you help me with that?   For the large error, I noticed you’re using linear element. I encountered the same large error when comparing it with Abaqus with FE_Q(1). But the error came down with less grids when I used higher finite element FE_Q<dim>(2). I remember the deflection of beam is a cubic function of coordinate. You may try that see if it improves.   Best, Michael   From: Alex Cumberworth Sent: Wednesday, June 9, 2021 9:54 AM To: deal.II User Group Subject: [deal.II] Re: Integrated material and spatial traction forces on boundary not equal   Hello,   I have attached the most recent version of my code here. I have tried to make setting boundary conditions in the parameter file more convenient for myself; you can set boundary domains, boundary conditions that use these domains, and stages if the displacement is large. However, there are not many comments, so you may want to just remove this part for your own purposes.   However, I have been a bit surprised in my comparisons at how fine a mesh is required to achieve convergence with the beam theory result. I am now using a beam that is 1 x 2 x 20, and using the subdivided rectangle helper function, I set the number of subdivisions to be 5 and 10 for the width and height, respectively. I then varied the number of subdivision in the length between 10 and 1000. The beam theory result is that the shear force has a magnitude of 0.001 for a displacement on the right side of 0.1. Even at 1000 subdivisions, the FEM result is 0.00113 (from 0.00129 at 500). The system has 200 000 degrees of freedom, and the result is still off by 13%. Is it expected that even to calculate a shear force in this simple problem that such a large number of degrees of freedom are required?     Best, Alex   On Monday, June 7, 2021 at 3:39:44 p.m. UTC+2 lian...@gmail.com wrote: Hi Alex, I'm learning deal.ii and trying do the similar verification. If it is possible for you to share the code with me?   Thank you! Michael On Tuesday, May 11, 2021 at 4:46:55 AM UTC-6 alexanderc...@gmail.com wrote: Hello,   As a test to validate my code, I am solving the equations for geometrically nonlinear elasticity (the Saint Venant-Kirchhoff model) for a beam with a small displacement boundary condition on the right end such that I can compare with Euler-Bernoulli beam theory. I can compare both the displacement and the shear force between the FEM solution and the beam theory solution. In my FEM integration, I output the normal and shear forces for both sides of the beam in both the material and spatial reference. The left and right sides are balanced, as expected, but the spatial and material forces are not quite equal.   Shouldn't it be the case that spatial and material force is the same? Here are the outputted forces for the right side   Right boundary material normal force: 0.0694169 Right boundary spatial normal force: 0.0724468 Right boundary material shear force: 0.152057 Right boundary spatial shear force: 0.152864   Further, beam theory gives a shear force with a magnitude of exactly 0.2. If I make the displacement smaller the FEM and beam theory shear forces do not converge. Is it expected for them to converge?   Below is the deformed system with the stress vectors on the faces included. The black grid is the deformed FEM solution, while the solid red is the beam theory solution. If there is an issue, I would guess it would be in the integration. In converting the material normal vector to the spatial reference, I first only applied the inverse transpose of the deformation gradient, and did not multiply by the determinant until calculating the force vector. I did this so that I can get the unit normal spatial vectors to add up and later average so that I have an average normal vector for the whole boundary face to calculate the normal and shear force vectors. I have pasted the function in below:   template <int dim> void SolveRing<dim>::integrate_over_boundaries() {     QGauss<dim - 1> quadrature_formula(fe.degree + 1);     FEFaceValues<dim> fe_face_values(             fe,             quadrature_formula,             update_values | update_gradients | update_quadrature_points |                     update_JxW_values | update_normal_vectors);     std::vector<Tensor<1, dim, double>> material_force(2);     std::vector<Tensor<1, dim, double>> spatial_force(2);     std::vector<Tensor<1, dim, double>> ave_material_normal(2);     std::vector<Tensor<1, dim, double>> ave_spatial_normal(2);     const FEValuesExtractors::Vector displacements(0);     for (const auto& cell: dof_handler.active_cell_iterators()) {         for (const auto face_i: GeometryInfo<dim>::face_indices()) {             const unsigned int boundary_id {cell->face(face_i)->boundary_id()};             if (not(boundary_id == 1 or boundary_id == 2)) {                 continue;             }             fe_face_values.reinit(cell, face_i);             std::vector<Tensor<1, dim, double>> normal_vectors {                     fe_face_values.get_normal_vectors()};             std::vector<Tensor<2, dim, double>> solution_gradients(                     fe_face_values.n_quadrature_points);             fe_face_values[displacements].get_function_gradients(                     present_solution, solution_gradients);             for (const auto q_i: fe_face_values.quadrature_point_indices()) {                 const Tensor<2, dim, double> grad_u {solution_gradients[q_i]};                 const Tensor<1, dim, double> material_normal {                         normal_vectors[q_i]};                 const Tensor<2, dim, double> grad_u_T {transpose(grad_u)};                 const Tensor<2, dim, double> green_lagrange_strain {                         0.5 * (grad_u + grad_u_T + grad_u_T * grad_u)};                 const Tensor<2, dim, double> piola_kirchhoff {                         lambda * trace(green_lagrange_strain) *                                 unit_symmetric_tensor<dim>() +                         2 * mu * green_lagrange_strain};                 ave_material_normal[boundary_id - 1] += material_normal;                 material_force[boundary_id - 1] += piola_kirchhoff *                                                    material_normal *                                                    fe_face_values.JxW(q_i);                 const Tensor<2, dim, double> deformation_grad {                         grad_u + unit_symmetric_tensor<dim>()};                 const double deformation_grad_det {                         determinant(deformation_grad)};                 const Tensor<2, dim, double> cauchy {                         deformation_grad * piola_kirchhoff *                         transpose(deformation_grad) / deformation_grad_det};                 Tensor<1, dim, double> spatial_normal {                         transpose(invert(deformation_grad)) * material_normal};                 spatial_force[boundary_id - 1] += cauchy * spatial_normal *                                                   fe_face_values.JxW(q_i) *                                                   deformation_grad_det;                 spatial_normal /= spatial_normal.norm();                 ave_spatial_normal[boundary_id - 1] += spatial_normal;             }         }     }     ave_material_normal[0] /= ave_material_normal[0].norm();     ave_spatial_normal[0] /= ave_spatial_normal[0].norm();     ave_material_normal[1] /= ave_material_normal[1].norm();     ave_spatial_normal[1] /= ave_spatial_normal[1].norm();     const Tensor<1, dim, double> left_material_normal_force = {             (material_force[0] * ave_material_normal[0]) *             ave_material_normal[0]};     const Tensor<1, dim, double> left_material_shear_force = {             material_force[0] - left_material_normal_force};     const Tensor<1, dim, double> left_spatial_normal_force = {             (spatial_force[0] * ave_spatial_normal[0]) * ave_spatial_normal[0]};     const Tensor<1, dim, double> left_spatial_shear_force = {             spatial_force[0] - left_spatial_normal_force};     const Tensor<1, dim, double> right_material_normal_force = {             (material_force[1] * ave_material_normal[1]) *             ave_material_normal[1]};     const Tensor<1, dim, double> right_material_shear_force = {             material_force[1] - right_material_normal_force};     const Tensor<1, dim, double> right_spatial_normal_force = {             (spatial_force[1] * ave_spatial_normal[1]) * ave_spatial_normal[1]};     const Tensor<1, dim, double> right_spatial_shear_force = {             spatial_force[1] - right_spatial_normal_force};     cout << "Left boundary material normal force: "          << left_material_normal_force.norm() << std::endl;     cout << "Right boundary material normal force: "          << right_material_normal_force.norm() << std::endl;     cout << "Left boundary material shear force: "          << left_material_shear_force.norm() << std::endl;     cout << "Right boundary material shear force: "          << right_material_shear_force.norm() << std::endl;     cout << "Left boundary spatial normal force: "          << left_spatial_normal_force.norm() << std::endl;     cout << "Right boundary spatial normal force: "          << right_spatial_normal_force.norm() << std::endl;     cout << "Left boundary spatial shear force: "          << left_spatial_shear_force.norm() << std::endl;     cout << "Right boundary spatial shear force: "          << right_spatial_shear_force.norm() << std::endl;   I will also paste in the core of my assembly loop:   fe_values.reinit(cell); cell_matrix = 0; cell_rhs = 0; std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); cell->get_dof_indices(local_dof_indices); const unsigned int n_independent_variables = local_dof_indices.size(); ADHelper ad_helper(n_independent_variables); ad_helper.register_dof_values(evaluation_point, local_dof_indices); const std::vector<ADNumberType>& dof_values_ad =         ad_helper.get_sensitive_dof_values(); ADNumberType energy_ad = ADNumberType(0.0); std::vector<Tensor<2, dim, ADNumberType>> old_solution_gradients(         fe_values.n_quadrature_points); fe_values[displacements].get_function_gradients_from_local_dof_values(         dof_values_ad, old_solution_gradients); for (const unsigned int q_index: fe_values.quadrature_point_indices()) {     const Tensor<2, dim, ADNumberType> grad_u {             old_solution_gradients[q_index]};     const Tensor<2, dim, ADNumberType> grad_u_T {transpose(grad_u)};     const Tensor<2, dim, ADNumberType> green_lagrange_strain_tensor {             0.5 * (grad_u + grad_u_T + grad_u_T * grad_u)};     ADNumberType t1 = lambda / 2 *                       std::pow(trace(green_lagrange_strain_tensor), 2);     ADNumberType t2 = mu * double_contract<0, 0, 1, 1>(                                    green_lagrange_strain_tensor,                                    green_lagrange_strain_tensor);     ADNumberType pi {t1 + t2};     energy_ad += pi * fe_values.JxW(q_index); } ad_helper.register_energy_functional(energy_ad); present_energy += ad_helper.compute_energy(); ad_helper.compute_residual(cell_rhs); cell_rhs *= -1.0; // RHS = - residual   Any thoughts would be greatly appreciated,   Alex -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group. To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/K-lMxbtZUdQ/unsubscribe. To unsubscribe from this group and all its topics, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4cade24b-48aa-432b-a590-a75535c0f2f2n%40googlegroups.com.   -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/24E2EF49-EA72-47BB-B0B3-7889B534CBD3%40hxcore.ol.