Hi Michael, I’m glad that you’re on the right track now.
> If you could give some advice or refer to some references on this topic, I > would appreciate that very mush. Unfortunately, this is where pretty much I hit the limit of what I can say on this topic as it would conflict with the expectations of my employer. There’s a little by-line on one approach in this paper https://www.sciencedirect.com/science/article/abs/pii/S0997753814001508 <https://www.sciencedirect.com/science/article/abs/pii/S0997753814001508> But as for other references, I’m hoping that someone else can give you some guidance there. (There’s at least one fairly common book in nonlinear mechanics that discusses the "linearisation of external virtual work”. Maybe a search in this direction might get you somewhere.) Another option would be to try the auto-differentiation tools to do all of the heavy lifting for you. With few restrictions, it would be compatible with any general traction load. At the very least and AD implementation it could help you validate any linearisation that you do by hand. You can take a look at steps 70 and 71 to get some idea as to what it might be able to do for you. Using the approach in step-71 is probably easiest to start off in your case. You could even limit it to just linearising the contribution that comes from this boundary / traction contribution. Best, Jean-Paul > On 4. Aug 2021, at 23:22, Michael Li <[email protected]> wrote: > > Hi Jean-Paul, > > Yes, I set up the boundary check when taking care of the traction. It worked > as expected using Physics::Transformations::nansons_formula(face_normal, F). > > To make a quick test, I did not linearize the load stiffness but just used > the load at last Newton-Raphson step, so it is like a mixture of N-R and > Picard iteration. The technique is not correct and there is some discrepancy > compared with the reference. But the result is better than not considering > the face area change at all. So in some sense, it tells me that the surface > transform is correct. I’m working on how to linearize an arbitrary surface > traction. I have found some references dealing with the normal(pressure) load > but not arbitrary surface load (the force is not normal to the face in my > case). If you could give some advice or refer to some references on this > topic, I would appreciate that very mush. > > You’re totally right, it converges only if a very small load step is used. > The residual grows up fast and the result is meaningless if the load > increment is large. > > Thanks for your great help! > > Michael > > > > From: Jean-Paul Pelteret <mailto:[email protected]> > Sent: Wednesday, August 4, 2021 12:19 AM > To: [email protected] <mailto:[email protected]> > Subject: Re: [deal.II] Computing the solution gradient at the quadrature > point on a face > > Hi Michael, > > The one thing that stands out about the snippet of code that you shared is > that there is no check to > if (cell->face(face)->at_boundary() == true && <ON BOUNDARY OF INTEREST>)... > Do you still have something like that? > > Otherwise, after a cursory glance, everything appears to be alright with > this. It mimics what you’d do to get the deformation gradient at a QP in a > cell. > > To debug further, what happens when you apply a very small load? Maybe its > not scaled appropriately for the stiffness of the material. Also, if you’ve > now added a deformation dependent part to the load (i.e. the area > transformation) then you might need to consistently linearise this in order > to attain convergence for large loads. > > Best, > Jean-Paul > > > On 2. Aug 2021, at 18:03, Michael Li <[email protected] > <mailto:[email protected]>> wrote: > > Hi Jean-Paul, > > Thank you for your hints. > > I initialized a local variable solution_grads_u_face_total with > fe_face_values_ref[u_fe].get_function_gradients to get the gradient of > displacement at the quadrature point on the face (see codes below). But the > gradient acquired is not sensible and the resultant deformation is totally > wrong. > > void assemble_neumann_contribution_one_cell(…) > { > ……………. > const FEValuesExtractors::Vector &u_fe = data.solid->u_fe; > > std::vector<Tensor<2, dim,NumberType> > > solution_grads_u_face_total(data.solid->qf_face.size()); > > scratch.fe_face_values_ref.reinit(cell, face); > > > scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total, > > solution_grads_u_face_total); > for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; > ++face) > { > for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; ++f_q_point) > { > const Tensor<2,dim,NumberType> &grad_u = > solution_grads_u_face_total[f_q_point]; > const Tensor<2,dim,NumberType> F = > Physics::Elasticity::Kinematics::F(grad_u); > …………. > } > ……….. > } > ……….. > } > > Can you tell me what’s wrong with the code above? > > Thank you! > Michael > > > From: Jean-Paul Pelteret <mailto:[email protected]> > Sent: Sunday, August 1, 2021 1:46 PM > To: [email protected] <mailto:[email protected]> > Subject: Re: [deal.II] Computing the solution gradient at the quadrature > point on a face > > Hi Michael, > > The problem here is that “scratch.solution_grads_u_total” is initialised with > the cell FEValues, and so has n_cell_quadrature_points entries. You need to > pass > "scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total, > scratch.solution_grads_u_total);” something that has size > n_face_quadrature_points. That’s what this error is informing you. > > > > As a side note, Physics::Transformations::nansons_formula() > <https://dealii.org/developer/doxygen/deal.II/namespacePhysics_1_1Transformations.html#aa82925b742c3708f625a48a7abc440bc> > can compute the transformation of the normal vector + area change for you, > should you like the traction direction to follow the moving boundary surface > as well. > > Best, > Jean-Paul > > > > On 1. Aug 2021, at 20:05, Michael Lee <[email protected] > <mailto:[email protected]>> wrote: > > fe_face_values.get_function_gradients returns a rank 2 tensor as the solution > is the displacement vector. > > Let me be more specific. In the code gallery " Quasi-Static Finite-Strain > Compressible Elasticity (The deal.II Library: The 'Quasi-Static Finite-Strain > Compressible Elasticity' code gallery program (dealii.org) > <https://www.dealii.org/current/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html> > ). It does not consider the follower force as described in the > documentation: " > // We specify the traction in reference configuration. > // For this problem, a defined total vertical force is applied > // in the reference configuration. > // The direction of the applied traction is assumed not to > // evolve with the deformation of the domain." > > So I made a modification to consider the effect of the deformed area on the > traction as follows according to the formula da/dA = J sqrt(N C^-1 N) . > > void > assemble_neumann_contribution_one_cell(const typename > DoFHandler<dim>::active_cell_iterator &cell, > ScratchData_ASM &scratch, > PerTaskData_ASM &data) > { > // Aliases for data referenced from the Solid class > const unsigned int &n_q_points_f = data.solid->n_q_points_f; > const unsigned int &dofs_per_cell = data.solid->dofs_per_cell; > const Time &time = data.solid->time; > const FESystem<dim> &fe = data.solid->fe; > const unsigned int &u_dof = data.solid->u_dof; > const FEValuesExtractors::Vector &u_fe = data.solid->u_fe; > > // Next we assemble the Neumann contribution. We first check to see it > the > // cell face exists on a boundary on which a traction is applied and add > // the contribution if this is the case. > for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; > ++face) > if (cell->face(face)->at_boundary() == true > && cell->face(face)->boundary_id() == 11) > { > scratch.fe_face_values_ref.reinit(cell, face); > > scratch.fe_face_values_ref[u_fe].get_function_gradients(scratch.solution_total, > > scratch.solution_grads_u_total); > > for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; > ++f_q_point) > { > // Note that the contributions to the right hand side vector > we > // compute here only exist in the displacement components of > the > // vector. > const double time_ramp = (time.current() / time.end()); > const double magnitude = 0.5*1.95*1000*0.32*0.32*time_ramp; > // (Total force) / (RHS surface area) > Tensor<1,dim> dir; > dir[0] = 1.0; > Tensor<1,dim> face_normal = > scratch.fe_face_values_ref.normal_vector(f_q_point); > > const Tensor<2,dim,NumberType> &grad_u = > scratch.solution_grads_u_total[f_q_point]; > const Tensor<2,dim,NumberType> F = > Physics::Elasticity::Kinematics::F(grad_u); > const SymmetricTensor<2,dim,NumberType> C = > Physics::Elasticity::Kinematics::C(F); > const SymmetricTensor<2,dim,NumberType> C_inv = invert(C); > const Tensor<1, dim> traction = magnitude*dir; > > for (unsigned int i = 0; i < dofs_per_cell; ++i) > { > const unsigned int i_group = > fe.system_to_base_index(i).first.first; > > if (i_group == u_dof) > { > const unsigned int component_i = > fe.system_to_component_index(i).first; > const double Ni = > scratch.fe_face_values_ref.shape_value(i, > f_q_point); > const double JxW = > scratch.fe_face_values_ref.JxW(f_q_point); > > data.cell_rhs(i) += (Ni * traction[component_i]) > * sqrt(face_normal * C_inv * > face_normal) > * JxW; > } > } > } > } > } > > > The errors I got are as follows: > If I set the parameters as: set Polynomial degree = 2 set Quadrature > order = 2 > An error occurred in line <645> of file > </home/michael/dealii-9.2.0/source/fe/fe_values.cc <http://fe_values.cc/>> in > function > void dealii::FEValuesViews::internal::do_function_derivatives(const > dealii::ArrayView<Number>&, const dealii::Table<2, dealii::Tensor<order, > spacedim> >&, const std::vector<typename dealii::FEValuesViews::Vector<dim, > spacedim>::ShapeFunctionData>&, std::vector<typename > dealii::ProductType<Number, dealii::Tensor<(order + 1), spacedim> >::type>&) > [with int order = 1; int dim = 3; int spacedim = 3; Number = double; typename > dealii::FEValuesViews::Vector<dim, spacedim>::ShapeFunctionData = > dealii::FEValuesViews::Vector<3, 3>::ShapeFunctionData; typename > dealii::ProductType<Number, dealii::Tensor<(order + 1), spacedim> >::type = > dealii::Tensor<2, 3, double>] > The violated condition was: > static_cast<typename ::dealii::internal::argument_type<void( typename > std::common_type<decltype(derivatives.size()), > decltype(n_quadrature_points)>::type)>::type>(derivatives.size()) == > static_cast<typename ::dealii::internal::argument_type<void( typename > std::common_type<decltype(derivatives.size()), > decltype(n_quadrature_points)>::type)>::type>(n_quadrature_points) > Additional information: > Dimension 8 not equal to 4. > > If I use set Polynomial degree = 2 set Quadrature order = 1 > An error occurred in line <335> of file > </home/michael/dealii-9.2.0/source/lac/sparse_direct.cc > <http://sparse_direct.cc/>> in function > void dealii::SparseDirectUMFPACK::factorize(const Matrix&) [with Matrix = > dealii::SparseMatrix<double>] > The violated condition was: > status == UMFPACK_OK > Additional information: > UMFPACK routine umfpack_dl_numeric returned error status 1. > > On Sunday, August 1, 2021 at 12:08:37 AM UTC-6 [email protected] > <http://gmail.com/> wrote: > The return type should be `Tensor<2,dim,NumberType>` shouldn't it? > > Hope this helps. If not, could describe how the code does not work? It it not > compiling? Is the result wrong? > > PM > > On Sunday, 1 August 2021 at 04:09:45 UTC+2 [email protected] > <http://gmail.com/> wrote: > Dear All, > > I want to do a surface integration which needs to evaluate the gradients of > the solution on the face. The following code does not work > > fe_face_values.get_function_gradients(solution, solution_grads); > for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; > ++f_q_point) > { > const Tensor<2,dim,NumberType> &grad_u = > solution_grads[f_q_point]; > ... > } > > Can anybody give a hint on that? > > Thanks, > Michael > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/3b018015-206e-41f7-9ff3-88c909927af7n%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/3b018015-206e-41f7-9ff3-88c909927af7n%40googlegroups.com?utm_medium=email&utm_source=footer>. > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/E4775A97-1A0A-48FF-9CE6-54ADD1448E53%40gmail.com > > <https://groups.google.com/d/msgid/dealii/E4775A97-1A0A-48FF-9CE6-54ADD1448E53%40gmail.com?utm_medium=email&utm_source=footer>. > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/2B729A00-167E-41F5-9363-7B4ADD1A8468%40hxcore.ol > > <https://groups.google.com/d/msgid/dealii/2B729A00-167E-41F5-9363-7B4ADD1A8468%40hxcore.ol?utm_medium=email&utm_source=footer>. > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/C5A123B1-E09B-4916-ADE4-D1DB8142B720%40gmail.com > > <https://groups.google.com/d/msgid/dealii/C5A123B1-E09B-4916-ADE4-D1DB8142B720%40gmail.com?utm_medium=email&utm_source=footer>. > > > -- > The deal.II project is located at http://www.dealii.org/ > <http://www.dealii.org/> > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > <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 [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/25249C1B-9F64-4922-92A6-6A1D5AF3B52D%40hxcore.ol > > <https://groups.google.com/d/msgid/dealii/25249C1B-9F64-4922-92A6-6A1D5AF3B52D%40hxcore.ol?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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/27C2A35B-D9ED-414B-98E0-722747C2A914%40gmail.com.
