Hi Alex,
I just got a chance to have Trilinos installed but the same error occurs
regarding *EnergyFunctional*. I got the following notice when trying to
test *step-31*:
Error! This tutorial requires a deal.II library that was configured with
the following options:
DEAL_II_WITH_TRILINOS = ON
However, the deal.II library found at /usr/local was configured with these
options
DEAL_II_WITH_TRILINOS = OFF
which conflict with the requirements.
-- Configuring incomplete, errors occurred!
Is there a way to let deal.II know Trillinos is installed? Do I have to
reinstall deal.ii?
Or, can I compile the code without using Trillinos with some modification?
Thanks!
Michael
On Friday, June 11, 2021 at 4:04:20 AM UTC-6 alexanderc...@gmail.com wrote:
> Hi Michael,
>
> I think that is the first point in my code that something from the
> automatic differentiation module is referred to, which is a part of
> trilinos. Was your version of deal.ii compiled with trilinos, and if so,
> was your version of trilinos compiled with sacado?
>
> To answer the question I posed, I did some more tests and found that
> refining in the width and height of the beam has no effect on the
> convergence of the shear force. I went down to just a single division in z
> and 2 in y (in other words there are just two faces when looking down the
> end of the beam), and could get the same level of agreement with beam
> theory with one order of magnitude fewer degrees of freedom. I guess I am
> still mildly surprised how refined I had to made the x direction, but
> perhaps if I also used adaptive refinement things would further improve.
>
> Best,
> Alex
>
> On Thursday, June 10, 2021 at 6:20:58 p.m. UTC+2 lian...@gmail.com wrote:
>
>> 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+un...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/4cade24b-48aa-432b-a590-a75535c0f2f2n%40googlegroups.com
>>
>> <https://groups.google.com/d/msgid/dealii/4cade24b-48aa-432b-a590-a75535c0f2f2n%40googlegroups.com?utm_medium=email&utm_source=footer>
>> .
>>
>>
>>
>
--
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/76a3d008-3e08-40c6-8887-04f3e3c9df8an%40googlegroups.com.
