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> 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> 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] 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] 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/ 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/3b018015-206e-41f7-9ff3-88c909927af7n%40googlegroups.com.
