Hi Hamed, I was in the same boat a few weeks ago. After reading through some of the other posts on this list (which Jean-Paul linked), I wrote the following method that you and maybe others will find useful.
It takes in a vector containing the list of components of a field that need point constraints (rigidBodyModeComponents) and adds a constraint where needed: https://github.com/prisms-center/phaseField/blob/next/src/matrixfree/boundaryConditions.cc (lines 43-75) Cheers! Steve On Thursday, October 20, 2016 at 7:38:19 PM UTC-4, Hamed Babaei wrote: > > Hi friends, > > For an elastic problem, I am going to apply zero boundary displacements > for three specific points on the center of -x, -y and -z planes of a cubic > domain. > I have already done this but for the boundary surface not a boundary point > (the same as incremental_boundary_displacement in step-18). > The following is what I wrote to do so which doesn't work properly. In > fact, it fixes the whole -x, -y and -z surfaces, not just for the three > points on them that I intended. > > template <int dim> > class BoundaryCondition : public Function<dim> > { > public: > BoundaryCondition (const int boundary_id); > virtual void vector_value (const Point<dim> &p, > Vector<double> &values) const; > virtual void vector_value_list (const std::vector<Point<dim> > > &points, > std::vector<Vector<double> > > &value_list) const; > private: > const int boundary_id; > > }; > template <int dim> > BoundaryCondition<dim>::BoundaryCondition (const int boundery_id) > : > Function<dim> (dim), > boundary_id(boundary_id) > > {} > template <int dim> > inline > void > BoundaryCondition<dim>::vector_value (const Point<dim> &p, > Vector<double> &values) const > { > Assert (values.size() == dim, > ExcDimensionMismatch (values.size(), dim)); > > Point<dim> point_x; > point_x(1) = 5; > point_x(2) = 5; > > Point<dim> point_y; > point_y(0) = 5; > point_y(2) = 5; > > Point<dim> point_z; > point_z(0) = 5; > point_z(1) = 5; > > > if (boundary_id ==0 && ((p-point_x).norm_square() > <(0.5e-9)*(0.5e-9))) > values(0) = 0; > else if (boundary_id ==2 && ((p-point_y).norm_square() < > (0.5e-9)*(0.5e-9))) > values(1)= 0; > else if (boundary_id ==4 && ((p-point_z).norm_square() < > (0.5e-9)*(0.5e-9))) > values(2)= 0; > > } > template <int dim> > void > BoundaryCondition<dim>::vector_value_list (const > std::vector<Point<dim> > &points, > std::vector<Vector<double> > > &value_list) const > { > const unsigned int n_points = points.size(); > Assert (value_list.size() == n_points, > ExcDimensionMismatch (value_list.size(), n_points)); > for (unsigned int p=0; p<n_points; ++p) > BoundaryCondition<dim>::vector_value (points[p], > value_list[p]); > } > > .....and in the constraint I have > added VectorTools::interpolate_boundary_values for every boundary plane, > for example for the -x plane it looks like : > > > VectorTools::interpolate_boundary_values(dof_handler, > boundary_id, > > BoundaryCondition<dim> (boundary_id), > constraints, > > fe.component_mask(x_displacement)); > > > I don't know why it doesn't recognize the if condition " if > (boundary_id ==0 && ((p-point_x).norm_square() <(0.5e-9)*(0.5e-9)))" !!! > > I was wondering if you know where I am making mistake, or if there is any > step in which this boundary condition has applied. > > Thanks, > Hamed > -- 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.
