Hi Bastien,
There are a few potential problems here.
1. You're only refining cells along the top edge of your domain.
Periodic boundaries can only have a difference of 1 refinement level
between pairs of faces.
2. You may not be colouring any of your boundaries. This should be a
tolerant check, something like "if (std::abs(face_center[1] -1.0) <=
1e-6)" to take into account computational errors.
3. The offset vector should probably be all zero since the vertices
align perfectly on either side of a hypercube. With an offset of 0.1 and a
regular global refinement, if I'm not mistaken its impossible to get two
vertices that are 0.1 units offset from one another. From the
documentation of the collect_periodic_faces function:
The offset is a vector tangential to the faces that is added to the
> location of vertices of the 'first' boundary when attempting to match them
> to the corresponding vertices of the 'second' boundary
I hope that this helps,
J-P
On Tuesday, June 7, 2016 at 6:47:51 PM UTC+2, Bastien Lauras wrote:
>
> Hi Daniel.
>
> Thank you for your answer, it worked well to enforce my periodic boundary
> conditions on the 1D beam.
> Thus, I'm back to the 2D neo-hookean case. Here is the code of my
> *make_grid()* function :
>
> template <int dim>
> void Solid<dim>::make_grid()
> {
> GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
>
> // We will refine the grid in five steps towards the inner circle of
> // the domain. See step 1 for details.
> for (unsigned int step=0; step<parameters.local_refinement_cycles;
> ++step)
> {
> typename Triangulation<dim>::active_cell_iterator
> cell_ref = triangulation.begin_active(),
> endc = triangulation.end();
> for (; cell_ref!=endc; ++cell_ref)
> {
> for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
> if (cell_ref->face(f)->at_boundary())
> {
> const Point<dim> face_center = cell_ref->face(f)->center();
> if (face_center[1] == 1) //Top
> cell_ref->set_refine_flag ();
> }
> }
>
> // Now that we have marked all the cells that we want refined, we let
> // the triangulation actually do this refinement.
> triangulation.execute_coarsening_and_refinement ();
> triangulation.clear_user_flags();
> }
>
> // We refine our mesh globally, at least once, to not too have a poor
> // refinement at the bottom of our square (like two cells)
> triangulation.refine_global(std::max (1U,
> parameters.global_refinement));
>
> // We mark the surfaces in order to apply the boundary conditions after
> typename Triangulation<dim>::active_cell_iterator cell =
> triangulation.begin_active(), endc = triangulation.end();
> for (; cell != endc; ++cell)
> for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
> if (cell->face(f)->at_boundary())
> {
> const Point<dim> face_center = cell->face(f)->center();
> if (face_center[1] == 0) //Bottom
> cell->face(f)->set_boundary_id (0);
> else if (face_center[1] == 1) //Top
> cell->face(f)->set_boundary_id (2);
> else if (face_center[0] == 0) //Left
> cell->face(f)->set_boundary_id (1);
> else if (face_center[0] == 1) //Right
> cell->face(f)->set_boundary_id (3);
> else //No way
> cell->face(f)->set_boundary_id (4);
> }
>
> std::vector<GridTools::PeriodicFacePair<typename
> Triangulation<dim>::cell_iterator> >
> periodicity_vector;
>
> const unsigned int direction = 0;
>
> FullMatrix<double> rotation_matrix(dim);
> rotation_matrix[0][0]=1.;
> rotation_matrix[1][1]=1.;
>
> Tensor<1, dim> offset;
> offset[0]=0.1;
>
> GridTools::collect_periodic_faces(triangulation, 1, 3, direction,
> periodicity_vector, offset,
> rotation_matrix);
>
> std::cout << "\nSize of periodicity_vector : " <<
> periodicity_vector.size();
> }
>
> The last output gives me "Size of periodicity_vector : 0". I don't
> understant why. Am I misusing the collect_periodic_faces function?
> Thank you very much.
>
> Best,
>
> Bastien
>
>
> On Friday, June 3, 2016 at 7:29:50 AM UTC-5, Daniel Arndt wrote:
>>
>> You are right, Bastien. The function is not implemented in 1d.
>> However, you can just ask for the boundary dofs using
>> DoFTools::extract_boundary_dofs [1].
>> This gives you two DoFs and you can just set the constraint using
>> constraint_matrix.add_line(dof_1)
>> constraint_matrix.add_entry(dof_1,dof_2,-1.)
>>
>> Best,
>> Daniel
>>
>> [1]
>> https://www.dealii.org/8.4.1/doxygen/deal.II/namespaceDoFTools.html#ad3067f335a97de429176178689222f3a
>>
>> Am Freitag, 3. Juni 2016 01:42:58 UTC+2 schrieb Bastien Lauras:
>>>
>>> Hi,
>>>
>>> Ok, after reviewing it, I am going to impose a constraint of the form :
>>> u(0,y) = u(1,y) + lambda and v(0,y) = v(1,y) (lambda is the total
>>> compression of my square). If I'm not mistaken, that is what Jean-Paul was
>>> advicing.
>>> I have two example, the first one is the 2D incompressible neo-hookean
>>> square, where I would like to implement the constraint just above. I don't
>>> think it will solve the problem I had but I'll try and will get back to you.
>>>
>>> The second example is a 1D-beam on a linear elastic foundation. The only
>>> degree of freedom is the vertical displacement. So everything is in one
>>> dimension. I want to impose the constraint y(0) = y(1) (only vertical
>>> displacement y allowed).
>>> Here is the declaration of my FESystem :
>>> fe(FE_Q<dim>(parameters.poly_degree), dim)
>>>
>>> This needs to be changed to Hermitian cubic elements, but this is
>>> another topic. My triangulation is :
>>> Triangulation<dim> triangulation;
>>> with dim = 1. I then construct it, refine it, and mark the "faces". But
>>> what are the faces in 1D?
>>> GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
>>> triangulation.refine_global(std::max (1U,
>>> parameters.global_refinement));
>>>
>>> const double length = GridTools::volume(triangulation);
>>> std::cout << "Grid:\n\t Length: " << length << std::endl;
>>>
>>> // We mark the sides of the beam in order to apply the boundary
>>> conditions after
>>> typename Triangulation<dim>::active_cell_iterator cell =
>>> triangulation.begin_active(), endc = triangulation.end();
>>> for (; cell != endc; ++cell)
>>> for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
>>> ++f)
>>> if (cell->face(f)->at_boundary())
>>> {
>>> const Point<dim> face_center = cell->face(f)->center();
>>> if (face_center[0] == 0) //Left
>>> cell->face(f)->set_boundary_id (1);
>>> else if (face_center[0] == 1) //Right
>>> cell->face(f)->set_boundary_id (2);
>>> else //No way
>>> cell->face(f)->set_boundary_id (3);
>>> // Is this working in 1D ?
>>> }
>>>
>>> Then I make my sparsity pattern and all the like, and begin my
>>> timesteps, and so on, my Newton-Raphson iterations. The function
>>> make_constraints is as follows :
>>> template <int dim>
>>> void Solid<dim>::make_constraints(const int &it_nr)
>>> {
>>> // Since the constraints are different at different Newton
>>> iterations, we
>>> // need to clear the constraints matrix and completely rebuild
>>> // it. However, after the first iteration, the constraints remain
>>> the same
>>> // and we can simply skip the rebuilding step if we do not clear it.
>>> if (it_nr > 1)
>>> return;
>>>
>>> constraints.clear();
>>> std::vector<GridTools::PeriodicFacePair<typename
>>> DoFHandler<dim>::cell_iterator> >
>>> periodicity_vector;
>>>
>>> const int direction = 0;
>>> GridTools::collect_periodic_faces(dof_handler_ref, 1, 2, direction,
>>> periodicity_vector);
>>>
>>> const FEValuesExtractors::Scalar y_displacement(0);
>>>
>>> //DoFTools::make_periodicity_constraints<DoFHandler<dim> >
>>> // (periodicity_vector, constraints/*,
>>> fe.component_mask(y_displacement)*/);
>>> DoFTools::make_periodicity_constraints (dof_handler_ref,1,2,
>>> direction,
>>> constraints,
>>> fe.component_mask(y_displacement));
>>> // I tried both implementations of this function, but noone is
>>> working.
>>> constraints.close();
>>> }
>>>
>>> When I try to compile my file, I get the following error :
>>> error: undefined reference to 'void
>>> dealii::DoFTools::make_periodicity_constraints<dealii::DoFHandler<1, 1>
>>> >(dealii::DoFHandler<1, 1> const&, unsigned char, unsigned char, int,
>>> dealii::ConstraintMatrix&, dealii::ComponentMask const&)'
>>>
>>> I think the problem comes from the fact I am working in 1D, doesn't it?
>>>
>>> Many thanks,
>>>
>>> Bastien
>>>
>>> On Wednesday, June 1, 2016 at 11:55:31 AM UTC-5, Daniel Arndt wrote:
>>>>
>>>> It is quite difficult to say what is going wrong without seeing the
>>>> code.
>>>> You probably want to check first if you really want to impose that kind
>>>> of boundary conditions and follow Jean-Paul's advice.
>>>>
>>>> If the problem persists, it would be really helpful if you can give us
>>>> a minimal compiling code that reproduces the problem.
>>>>
>>>> Best,
>>>> Daniel
>>>>
>>>> Am Mittwoch, 1. Juni 2016 18:05:32 UTC+2 schrieb Bastien Lauras:
>>>>>
>>>>> Hi Daniel,
>>>>>
>>>>> So, I'm modelling a unit square, with an incompressible neo-hookean
>>>>> material. The formulation is (for the density of energy) : c_1 * (I_1 -
>>>>> 3)
>>>>> - p (I_3 - 1) where c_1 = mu / 2, p is a lagrangian parameter, and I_1
>>>>> and
>>>>> I_3 are the first and third invariant of C = F^T * F.
>>>>> For me to see that my periodic boundary conditions are not applied, I
>>>>> put a non-uniform mu aver the material (graded material). In the example
>>>>> attached, I blocked the horizontal displacement (x_displacement) over
>>>>> left
>>>>> and right faces, and blocked vertical displacement (y_displacement) over
>>>>> the bottom face. And I (try to) apply a periodic boundary condition for
>>>>> the
>>>>> vertical displacement over left and right faces.
>>>>> Attached is the output given to me after the first timestep (all
>>>>> timesteps are the same since I'm not applying any displacement nor force
>>>>> constraint).
>>>>>
>>>>> We can obviously see that vertical displacement is not the same on the
>>>>> left and right faces. Any idea of where it could come from?
>>>>>
>>>>> Many thanks!
>>>>> Best,
>>>>>
>>>>> Bastien
>>>>>
>>>>>
>>>>> On Wednesday, June 1, 2016 at 5:43:48 AM UTC-5, Daniel Arndt wrote:
>>>>>>
>>>>>> Bastian,
>>>>>>
>>>>>> can you provide us with a minimal example that shows your problem?
>>>>>> What you are trying to do should be possible even if it might not be the
>>>>>> right
>>>>>> thing to do in your situation.
>>>>>>
>>>>>> Best,
>>>>>> Daniel
>>>>>>
>>>>>> Am Mittwoch, 1. Juni 2016 02:10:36 UTC+2 schrieb Bastien Lauras:
>>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> I am working on a 2D square, a incompressible neo-hookean material.
>>>>>>> I've used step 44 and modified it to work on a two field formulation (
>>>>>>> *u* and p).
>>>>>>> I am trying to implement periodic boundary conditions for the
>>>>>>> vertical displacement on lateral faces (right and left).
>>>>>>> My left face has boundary indicator 1, and my right face has
>>>>>>> boundary indicator 3.
>>>>>>>
>>>>>>> In the make_constraints function, I've added, after
>>>>>>> constraints.clear() :
>>>>>>>
>>>>>>> DoFTools::make_hanging_node_constraints (dof_handler_ref,
>>>>>>> constraints); // I have local refinement in my code
>>>>>>>
>>>>>>> std::vector<GridTools::PeriodicFacePair<typename
>>>>>>> DoFHandler<dim>::cell_iterator> >
>>>>>>> periodicity_vector;
>>>>>>>
>>>>>>> const unsigned int direction = 0; // I've tried 0 and 1, I
>>>>>>> didn't understan what it really was
>>>>>>>
>>>>>>> GridTools::collect_periodic_faces(dof_handler_ref, 1, 3,
>>>>>>> direction,
>>>>>>> periodicity_vector);
>>>>>>>
>>>>>>> const FEValuesExtractors::Scalar x_displacement(0);
>>>>>>> const FEValuesExtractors::Scalar y_displacement(1);
>>>>>>>
>>>>>>> DoFTools::make_periodicity_constraints<DoFHandler<dim> >
>>>>>>> (periodicity_vector, constraints,
>>>>>>> fe.component_mask(y_displacement));
>>>>>>>
>>>>>>>
>>>>>>> And then I apply Dirichlet boundary conditions on face 1 and 3 on
>>>>>>> the horizontal displacement (x_displacement).
>>>>>>> I've also added the command :
>>>>>>>
>>>>>>> DoFTools::make_hanging_node_constraints (dof_handler_ref,
>>>>>>> constraints);
>>>>>>>
>>>>>>> before making my sparsity pattern.
>>>>>>> My triangulation is not a parallel::distributed one. I'm working on
>>>>>>> multithreads with the workstream::run of step 44.
>>>>>>>
>>>>>>>
>>>>>>> The problem is that my periodic boundary conditions are not
>>>>>>> enforced. The make_periodicity_constraints command does not change the
>>>>>>> number of constrained degrees of freedom (that I compute using a loop
>>>>>>> over
>>>>>>> the DOFs and calling constraints.is_constrained(i)).
>>>>>>> And moreover, I can see on my outputs that there is definitely no
>>>>>>> periodicity on my vertical displacement.
>>>>>>>
>>>>>>> Where can the problem come from?
>>>>>>>
>>>>>>> Thanks a lot for your help beforehand!
>>>>>>>
>>>>>>> Bastien Lauras
>>>>>>>
>>>>>>>
>>>>>>
--
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.
