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.
