In case that link dies, here's the snippet of code:
// Set constraints to pin the solution if there are no Dirichlet BCs for a
component of a variable in an elliptic equation
template <int dim>
void MatrixFreePDE<dim>::setRigidBodyModeConstraints( std::vector<int>
rigidBodyModeComponents, ConstraintMatrix * constraints, DoFHandler<dim>*
dof_handler){
if ( rigidBodyModeComponents.size() > 0 ){
// Choose the point where the constraint will be placed. Must be the
coordinates of a vertex.
dealii::Point<dim> target_point(0,0);
unsigned int vertices_per_cell=GeometryInfo<dim>::vertices_per_cell;
// Loop over each locally owned cell
typename DoFHandler<dim>::active_cell_iterator cell=
dof_handler->begin_active(), endc = dof_handler->end();
for (; cell!=endc; ++cell){
if (cell->is_locally_owned()){
for (unsigned int i=0; i<vertices_per_cell; ++i){
// Check if the vertex is the target vertex
if (target_point.distance (cell->vertex(i)) <
1e-2 * cell->diameter()){
// Loop through the list of components with
rigid body modes and add an inhomogeneous constraint for each
for (unsigned int component_num = 0;
component_num < rigidBodyModeComponents.size(); component_num++){
unsigned int
nodeID=cell->vertex_dof_index(i,component_num);
constraints->add_line(nodeID);
constraints->set_inhomogeneity(nodeID,0.0);
}
}
}
}
}
}
}
On Saturday, October 22, 2016 at 6:49:41 AM UTC-4, Stephen DeWitt wrote:
>
> 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.
