Dear Wolfgang, It should not be zero. In addition, the local quadrature point is for the cell.
Thanks, Benhour On Thursday, December 8, 2016 at 8:16:03 PM UTC-6, [email protected] wrote: > > Dear All, > I still have problem with implying boundary condition on the curved > domain. consider we produce quarter of a ring in 2D with the function > "quarter_hyper_shell" and want to impose Neumann boundary condition on the > outer surface. In fact we have a force normal to this surface. The other > boundary conditions for straight lines are axisymmetry boundary condition. > What I did is given as follow: > > 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); > for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point) > { > const double time_ramp = (time.current() / time.end()); > const double eta_orderM = > (local_quadrature_point_history[f_q_point].get_eta_boundary_value_mechanic()) > * time_ramp; *// The multiplication of Neumann value (for example heat > flux in heat equation or traction force normal to the surface in mechanic) > and time_ramp*. > static const Tensor<1, dim> normals = > scratch.fe_face_values_ref.normal_vector(f_q_point); * // Extracting the > normal vectors of the selected face* > const Tensor<1, dim> neumann_value = eta_orderM * normals; > > 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 JxW = scratch.fe_face_values_ref.JxW(f_q_point); > data.cell_rhs(i) -= neumann_value[component_i] * > scratch.fe_face_values_ref.shape_value(i, f_q_point) * JxW; > } > } > } > } > The problem is that with changing the sign of the data.cell_rhs(i), the > direction of the force did not change. In fact, The ring is under > compression, however, in the result, it seems that it is under tension. It > would be very kind of you if you tell me which part I am making mistake. I > know this website is not for code reviewing, but I really appreciate your > time and help. > > Bests, > Benhour > -- 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]. For more options, visit https://groups.google.com/d/optout.
