Dear all,
I am implementing a very simple linear elastic problem, and I've succeeded
in using Dirichlet conditions. However, reading the tutorial I am unable to
implement a Neumann condition. In my simple example I have something like
this:
^^^^^^^^^^ f
+----------+
| |
| |
| |
| |
+----------+
//////////// clamped
I don't understand how I can implement the integral on the top boundary, at
the moment I've written this code (which does not compile):
// Force (neumann bc), no force for now
ConstantFunction<2> right_hand_side(std::vector<double>({ 0.0, 0.0 }));
// Force on boundary (neumann bc)
ConstantFunction<2> right_face_side(std::vector<double>({ 0.0, 111.0
}));
// For each cell
for (auto cell : dof_handler.active_cell_iterators())
{
// Zero local matrix and vector
cell_matrix = 0;
cell_rhs = 0;
// Compute FE values for this cell
fe_values.reinit(cell);
// Compute lambda and mu (in case we can use custom functions)
lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
mu.value_list(fe_values.get_quadrature_points(), mu_values);
// Neumann force vector
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
// ...assemble lhs matrix and rhs vector
// Neumann (boundary forces, 2 is top boundary)
for (unsigned int face_number = 0; face_number < GeometryInfo<dim>::
faces_per_cell; ++face_number)
{
if (cell->face(face_number)->at_boundary() && (cell->face
(face_number)->boundary_id() == 2))
{
const unsigned int n_face_q_points =
face_quadrature_formula.size();
fe_face_values.reinit(cell, face_number);
auto qps = fe_face_values.get_quadrature_points();
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
const double neumann_value =
right_face_side.value(fe_face_values.quadrature_point(q_point));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) += (neumann_value *
fe_face_values.shape_value(i,
q_point) *
fe_face_values.JxW(q_point));
}
}
}
} // Boundary forces
I have also tried to use the normal vector, but weirdly enough, I don't get
a tension, but a shear with the following code:
const double neumann_value = (right_face_side.value(qps[
q_point], 0) * fe_face_values.normal_vector(q_point)[0])
+ (right_face_side.value(qps[
q_point], 1) * fe_face_values.normal_vector(q_point)[1]);
How can I implement this?
Thanks!
--
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.
