Dear Qin, If I understand this statement
> I have a vector named traction[dof_handler.n_dofs()], which gives the > various traction on each dof correctly, you have a pre-assembled global traction vector (of size n_dofs) that you wish to add to the system. Why is it that you want to “assemble” this vector again? Can you not just call system_rhs += traction_vector? If the Neumann boundaries overlapped with the Dirichlet boundaries, then you might be able to first call constraints.set_zero(traction_vector) so zero out the contributions from the overlapping DoFs. If that’s not a viable solution then one problem in the code that you’ve shown is that your assembly itself is incorrect. You need to extract the traction vector at the face quadrature point and assemble with that. Currently you’re performing assembly at each face quadrature point using the value of some element of the total traction vector as located at a support point. Remember that if you were to integrate some local traction vectors over the Neumann boundaries then you should end up with the global traction vector that you have precomputed. You can get the quadrature point values of the local traction vector from the global vector using FEValues::get_function_values() <https://www.dealii.org/9.0.0/doxygen/deal.II/classFEValuesBase.html#a3f84ee49f0fbb32bfe01b2f4a473bb47> or a similar call. I hope that this helps. Best, Jean-Paul > On 07 Oct 2018, at 17:29, Lucky Qin <1296944...@qq.com> wrote: > > Hello everyone, > > I would like to apply a traction boundary condition on the surface of a > cantilever beam. Suppose I have a vector named > traction[dof_handler.n_dofs()], which gives the various traction on each dof. > I implemented the boundary condition as follows: > > > for ( unsigned int i = 0; i < dofs_per_cell; ++i ) > { > const unsigned int component_i = > fe.system_to_component_index( i ).first; > for ( unsigned int face = 0; face < > GeometryInfo<dim>::faces_per_cell; ++face ) > { > if ( cell->face( face )->at_boundary() == true && > cell->face( face )->boundary_id() != 0 ) //boundary_id() == 0 means the > Dirichlet boundary condition > { > fe_face_values.reinit( cell, face ); > > std::vector<unsigned int> local_face_dof_indices( > dofs_per_face ); > cell->face( face )->get_dof_indices( > local_face_dof_indices ); > > > assert( dofs_per_face == > fe_face_values.n_quadrature_points * dim ); > > > for ( unsigned int q = 0; q < > fe_face_values.n_quadrature_points; ++q ) > { > cell_rhs( i ) += traction( > local_face_dof_indices[q * dim + component_i] ) * fe_face_values.shape_value( > i, q ) * fe_face_values.JxW( q ); > } > } > } > } > > At present, I set the traction to be nonzero in only one of the surfaces of a > cantilever beam to make quantitative comparison with the previous literature. > However, the results do not seem to be correct. > Another thing I would like to know is what this traction may mean in this > code. Is this pressure or force? > I would like to ask for some help from you about these two questions. > > Best regards， > Qin > > -- > 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 dealii+unsubscr...@googlegroups.com > <mailto:dealii+unsubscr...@googlegroups.com>. > For more options, visit https://groups.google.com/d/optout > <https://groups.google.com/d/optout>. -- 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. For more options, visit https://groups.google.com/d/optout.
