The problem with using VectorTools::interpolate_boundary_values is that it takes a function f(x,y,..), so if my input boundary condition is discrete, as one value for each original node, i would have to test the nodes positions in order to apply the correct boundary values inside this function. Also, if i refine the mesh, it wouldn't be possible to apply boundary values to the new nodes because they weren't in the initial list.
One way to solve this would be to keep two separate meshes, one to hold the initial boundary values and another to keep the refined mesh. So each step i would make an interpolation on the first one to get boundary values for the new nodes created at the refined mesh. 2016-10-24 19:51 GMT-02:00 Daniel Arndt <d.ar...@math.uni-goettingen.de>: > Tulio, > > If you want to prescribe Dirichlet boundary conditions on a part of the > boundary of your domain, you can just set some boundary_id to the > respective cell boundary faces[1] and call > VectorTools::interpolate_boundary_values[2]. > The boundary_ids are inherited when you refine the mesh. > If however you need to constrain an internal DoF, you have to do this > manually. > > Best, > Daniel > > [1] https://dealii.org/8.4.1/doxygen/deal.II/group__boundary.html# > ga09db1c877a02021474d1e7eaa71a2958 > [2] https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html# > ga09828fb0c70a95462718f994f8b8ab5a > > > Am Montag, 24. Oktober 2016 12:01:59 UTC+2 schrieb Tulio Ligneul: >> >> Hi, taking advantage of this topic, i'd like to ask a related question. >> >> In a 2D mesh, suppose i want to enforce a dirichlet boundary for which i >> have its nodal values. Edge's values could off course be obtained trough >> its nodes interpolation. >> >> In this context, I guess i could use an approach like Stephen's, but if i >> have too many boundary nodes it will probably be slow. Also, if i use an >> adaptive mesh and a new node is created on the boundary, it will be >> difficult to find this new node and apply a boundary constraint to it. Is >> there any easy way to do this? >> >> Em domingo, 23 de outubro de 2016 13:55:15 UTC-2, Hamed Babaei escreveu: >>> >>> Dear Stephen, >>> Thank you very much for your incredible help. >>> >>> On Saturday, October 22, 2016 at 6:43:29 AM UTC-5, Stephen DeWitt wrote: >>>> >>>> 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 a topic in the > Google Groups "deal.II User Group" group. > To unsubscribe from this topic, visit https://groups.google.com/d/ > topic/dealii/iVrxwM1b-Cc/unsubscribe. > To unsubscribe from this group and all its topics, send an email to > dealii+unsubscr...@googlegroups.com. > For more options, visit 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.
