Re: [deal.II] Re: Setting up Periodic Boundary Condition for Heat Transfer

[JT Hu]

Hi Daniel,

I am using 8.3 version but the step-45 in the example file is
outdated. I think we will still get to the older step 45 version if we
google "deal ii" and go to the tutorial instead of the updated one. Is
that intentional or a bug to fix?
It figured out this specific problem - I forgot to include the
GridTools class.. I do have some new questions.
I noticed that the compiler gave me "no matching function for call to
‘make_periodicity_constraints" when using dof_handler as the first
input to the function, but had no error if I used triangulation
instead. According to the step 45 tutorial, it seems that I need both.
Is that right? Do you have any idea what the problem is? I always
assumed the default face masks and offset.

Best,
Jingtian

On Thu, Mar 3, 2016 at 3:58 AM, Daniel Arndt
<d.a...@math.uni-goettingen.de> wrote:
> Jingtian,
>
> "collect_period_faces" is a member of GridTools since version 8.1.0.
> Therefore, you don't need the current developer version or any pre-release.
> Which version are you using so far? Does the current step-45 compile for
> you?
>
> Yes, in step-45 the standard boundary_ids are used, i.e. the coincide with
> the face numbers in GeometryInfo.
>
> Best,
> Daniel
>
>
>
> Am Donnerstag, 3. März 2016 00:36:23 UTC+1 schrieb JT Hu:
>>
>> Hi Daniel,
>>
>> I received the compile error ‘collect_periodic_faces’ is not a member
>> of ‘dealii::GridTools’. Is that because I did not upgrade to the 8.4
>> pre-release?
>>
>> Also, I am a little confused when we need to use the
>> parallel::distributed::Triangulation::add_periodicity function.
>> I notice from step 45 that when calling:
>> GridTools::collect_periodic_faces(dof_handler, 2, 3, direction,...)
>> The boundary_ids are assumed to be 2 and 3, but I did not find where
>> those are defined in the function. Are those ids assumed to be the
>> same as the face numbers described in the structGeometryInfo page?
>>
>> Best,
>> Jingtian
>>
>> On Tue, Mar 1, 2016 at 10:22 AM, Daniel Arndt
>> <d.a...@math.uni-goettingen.de> wrote:
>> > Hey Jingtian,
>> >
>> > You can call collect_periodic_faces multiple times with the same
>> > periodicity_vector.
>> > Periodicity constraints are constructed for each PeriodicFacePair
>> > consecutively. Therefore it doesn't matter,
>> > whether you first use make_periodicity_constraints multiple times for
>> > each
>> > periodic direction separately or once with
>> > PeriodicFacePairs for both periodic directions.
>> >
>> > Best,
>> > Daniel
>> >
>> > Am Dienstag, 1. März 2016 16:35:09 UTC+1 schrieb JT Hu:
>> >>
>> >> Hi Daniel,
>> >>
>> >> Thanks for the prompt response! Yeah it looks much easier now. I have
>> >> one question for extending to 3D. Do we need to create two
>> >> periodicity_vector variables, each for one periodic direction, add
>> >> them separately by triangulation.add_periodicity, and finally call
>> >> DoFTools::make_periodicity_constraints?
>> >>
>> >> Best,
>> >> Jingtian
>> >>
>> >> On Fri, Feb 26, 2016 at 4:07 AM, Daniel Arndt
>> >> <d.a...@math.uni-goettingen.de> wrote:
>> >> > Jingtian,
>> >> >
>> >> > the version of step-45 you are referring to is outdated. If I were
>> >> > you,
>> >> > I
>> >> > would just use GridTools::collect_periodic_faces
>> >> > and DoFTools::make_periodicity_constraints as explained in the
>> >> > current
>> >> > step-45
>> >> > https://www.dealii.org/developer/doxygen/deal.II/step_45.html
>> >> > This is more flexible and not restricted to Q1 elements.
>> >> >
>> >> > Best regards,
>> >> > Daniel
>> >> >
>> >> >
>> >> > Am Freitag, 26. Februar 2016 01:12:47 UTC+1 schrieb JT Hu:
>> >> >>
>> >> >> Hi All,
>> >> >>
>> >> >> I am trying to build a simple 3D thermal model combining step 26
>> >> >> with
>> >> >> step 45. I pretty much just added the make_periodicity_constraints
>> >> >> function from step 45 to step 26, and made the changes for 3D
>> >> >> rectangular grid. The code compiles, but will run in to errors when
>> >> >> calling the constraints.add_entry functions:
>> >> >>
>> >> >> line != column
>> >> >>
>> >> >> and the additional message was: cannot constrain a degree of freedom
>> >> >> to
>> >> >> itself.
>> >> >>
>> >> >> Do you have any suggestions to debug? Did I misuse the add_line
>> >> >> function?
>> >> >>
>> >> >> I attached the .cc file and the specific function.
>> >> >>
>> >> >> Thanks for any suggestion!
>> >> >>
>> >> >> Best,
>> >> >> Jingtian
>> >> >>
>> >> >>
>> >> >>
>> >> >>
>> >> >> template<int dim>
>> >> >> void HeatEquation<dim>::make_periodicity_constraints()
>> >> >> {
>> >> >> Assert (dim == 3, ExcNotImplemented()); //PBC implemented only
>> >> >> for
>> >> >> 3D
>> >> >> case
>> >> >>
>> >> >>
>> >> >> std::map<unsigned int, double> dof_locations_finite_dim1;
>> >> >> //First
>> >> >> number is the infinite direction index, Second number is the
>> >> >> coordinate value to store
>> >> >> std::map<unsigned int, double> dof_locations_finite_dim2;
>> >> >> std::map<unsigned int, double> dof_locations_inf_dim1; //First
>> >> >> number is the infinite direction index, Second number is the
>> >> >> coordinate value to store
>> >> >> std::map<unsigned int, double> dof_locations_inf_dim2;
>> >> >> for (typename DoFHandler<dim>::active_cell_iterator cell =
>> >> >> dof_handler.begin_active ();
>> >> >> cell != dof_handler.end (); ++cell)
>> >> >> if (cell->at_boundary ())
>> >> >> {
>> >> >> if(cell->face( 2*INFINITE_DIRECTION_1 + 1 )->at_boundary
>> >> >> ())
>> >> >> // 'Face' a n dimensional hypercube has 2^(n-1) vertices
>> >> >> for (int i_vertex = 0; i_vertex < std::pow(2,dim-1);
>> >> >> i_vertex++)
>> >> >> {
>> >> >> dof_locations_finite_dim1[cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 + 1 )->vertex_dof_index(i_vertex, 0)]
>> >> >> = cell->face( 2*INFINITE_DIRECTION_1 + 1
>> >> >> )->vertex(i_vertex)[FINITE_DIRECTION];
>> >> >> // The cross numbering here is not an error. For faces on
>> >> >> periodic cell boundary in INFINITE_DIRECTION_1, the
>> >> >> // INFINITE_DIRECTION_2 coordinate value should be the same
>> >> >> for the matching face.
>> >> >> dof_locations_inf_dim1[cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 + 1 )->vertex_dof_index(i_vertex, 0)]
>> >> >> = cell->face( 2*INFINITE_DIRECTION_1 + 1
>> >> >> )->vertex(i_vertex)[INFINITE_DIRECTION_2];
>> >> >> }
>> >> >> else if(cell->face( 2*INFINITE_DIRECTION_2 + 1
>> >> >> )->at_boundary
>> >> >> ())
>> >> >> for (int i_vertex = 0; i_vertex < std::pow(2,dim-1);
>> >> >> i_vertex++)
>> >> >> {
>> >> >> dof_locations_finite_dim2[cell->face(
>> >> >> 2*INFINITE_DIRECTION_2 + 1 )->vertex_dof_index(i_vertex, 0)]
>> >> >> = cell->face( 2*INFINITE_DIRECTION_2 + 1
>> >> >> )->vertex(i_vertex)[FINITE_DIRECTION];
>> >> >> // The cross numbering here is not an error. For faces on
>> >> >> periodic cell boundary in INFINITE_DIRECTION_2, the
>> >> >> // INFINITE_DIRECTION_1 coordinate value should be the same
>> >> >> for the matching face.
>> >> >> dof_locations_inf_dim2[cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 + 1 )->vertex_dof_index(i_vertex, 0)]
>> >> >> = cell->face( 2*INFINITE_DIRECTION_1 + 1
>> >> >> )->vertex(i_vertex)[INFINITE_DIRECTION_1];
>> >> >> }
>> >> >> }
>> >> >>
>> >> >> for (typename DoFHandler<dim>::active_cell_iterator cell =
>> >> >> dof_handler.begin_active ();
>> >> >> cell != dof_handler.end (); ++cell)
>> >> >> if (cell->at_boundary ())
>> >> >> {
>> >> >>
>> >> >> if (cell->face (2*INFINITE_DIRECTION_1)->at_boundary ())
>> >> >> {
>> >> >>
>> >> >> for (unsigned int face_vertex = 0;
>> >> >> face_vertex<std::pow(2,dim-1); ++face_vertex)
>> >> >> {
>> >> >> constraints.add_line (cell->face(
>> >> >> 2*INFINITE_DIRECTION_1
>> >> >> )->vertex_dof_index (face_vertex, 0));
>> >> >>
>> >> >> std::map<unsigned int, double>::const_iterator p =
>> >> >> dof_locations_finite_dim1.begin();
>> >> >> // Look through the map dof_locations_finite_dim1 to
>> >> >> find a vertex on the other side
>> >> >> // with matching coordinates ( finite direction and
>> >> >> the
>> >> >> other infinite direction.
>> >> >> double temp_cell_finite_coordinate = cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 )->vertex(face_vertex)[FINITE_DIRECTION];
>> >> >> double temp_cell_inf_coordinate = cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 )->vertex(face_vertex)[INFINITE_DIRECTION_2];
>> >> >> for (; p != dof_locations_finite_dim1.end(); ++p)
>> >> >> {
>> >> >> double mapped_vertex_index = p -> first;
>> >> >> double mapped_vertex_finite_coordinate = p ->
>> >> >> second;
>> >> >> double mapped_vertex_inf_coordinate =
>> >> >> dof_locations_inf_dim1[mapped_vertex_index];
>> >> >> if (std::fabs(mapped_vertex_finite_coordinate -
>> >> >> temp_cell_finite_coordinate) < INFINITESIMAL_ERROR
>> >> >> && std::fabs( mapped_vertex_inf_coordinate -
>> >> >> temp_cell_inf_coordinate ) < INFINITESIMAL_ERROR )
>> >> >> {
>> >> >> constraints.add_entry (cell->face(
>> >> >> 2*INFINITE_DIRECTION_1 )->vertex_dof_index (face_vertex, 0),
>> >> >> p->first, 1.0);
>> >> >> break;
>> >> >> }
>> >> >> Assert (p != dof_locations_finite_dim1.end(),
>> >> >> ExcMessage ("No corresponding degree of
>> >> >> freedom
>> >> >> was found!"));
>> >> >> }
>> >> >> }
>> >> >> }
>> >> >> else if (cell->face (2*INFINITE_DIRECTION_2)->at_boundary
>> >> >> ())
>> >> >> {
>> >> >>
>> >> >> for (unsigned int face_vertex = 0;
>> >> >> face_vertex<std::pow(2,dim-1); ++face_vertex)
>> >> >> {
>> >> >> constraints.add_line (cell->face(
>> >> >> 2*INFINITE_DIRECTION_2
>> >> >> )->vertex_dof_index (face_vertex, 0));
>> >> >>
>> >> >> std::map<unsigned int, double>::const_iterator p =
>> >> >> dof_locations_finite_dim2.begin();
>> >> >> // Look through the map dof_locations_finite_dim2 to
>> >> >> find a vertex on the other side
>> >> >> // with matching coordinates ( finite direction and
>> >> >> the
>> >> >> other infinite direction.
>> >> >> double temp_cell_finite_coordinate = cell->face(
>> >> >> 2*INFINITE_DIRECTION_2 )->vertex(face_vertex)[FINITE_DIRECTION];
>> >> >> double temp_cell_inf_coordinate = cell->face(
>> >> >> 2*INFINITE_DIRECTION_2 )->vertex(face_vertex)[INFINITE_DIRECTION_1];
>> >> >> for (; p != dof_locations_finite_dim2.end(); ++p)
>> >> >> {
>> >> >> double mapped_vertex_index = p -> first;
>> >> >> double mapped_vertex_finite_coordinate = p ->
>> >> >> second;
>> >> >> double mapped_vertex_inf_coordinate =
>> >> >> dof_locations_inf_dim1[mapped_vertex_index];
>> >> >> if (std::fabs(mapped_vertex_finite_coordinate -
>> >> >> temp_cell_finite_coordinate) < INFINITESIMAL_ERROR
>> >> >> && std::fabs( mapped_vertex_inf_coordinate -
>> >> >> temp_cell_inf_coordinate ) < INFINITESIMAL_ERROR )
>> >> >> {
>> >> >> constraints.add_entry (cell->face(
>> >> >> 2*INFINITE_DIRECTION_2 )->vertex_dof_index (face_vertex, 0),
>> >> >> p->first, 1.0);
>> >> >> break;
>> >> >> }
>> >> >> Assert (p != dof_locations_finite_dim2.end(),
>> >> >> ExcMessage ("No corresponding degree of
>> >> >> freedom
>> >> >> was found!"));
>> >> >> }
>> >> >> }
>> >> >> }
>> >> >> }
>> >> >>
>> >> >> }
>> >> >
>> >> > --
>> >> > 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+un...@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+un...@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+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Daniel Arndt]

Jingtiang,

I think we will still get to the older step 45 version if we
google "deal ii" and go to the tutorial instead of the updated one. Is
that intentional or a bug to fix?

Searching for "step-45 deal.II" I get to the developer version. But you are right that in 8.3.0 we still used the old step-45.

I noticed that the compiler gave me "no matching function for call to
‘make_periodicity_constraints" when using dof_handler as the first
input to the function, but had no error if I used triangulation
instead. According to the step 45 tutorial, it seems that I need both.
Is that right? Do you have any idea what the problem is? I always
assumed the default face masks and offset.

Which version of make_periodicity_constraints are you referring to?
For the one using periodic_faces the type of periodic_faces has to match with the one you use in collect_periodic_faces.
In fact, you need an object that knows about the degrees of freedoms for calling make_periodicity_constraints.
If you successfully use collect_periodic_faces, make_periodicity_constraints should also be successful assuming you are using continuous elements.

Best,
Daniel

[JT Hu]

Hi Daniel,

Yeah you are right - I messed up the types. Now it compiles with no
problem, but sometimes encounter a run time error: " dof_index !=
line->line" and "Cycle in constraints detected!". Do you have any
ideas why? It happened at the 35th time step, this time, but did not
happen last time while running with different parameters.

Best,
Jingtian


An error occurred in line <289> of file
</home/jhh254/dealii-source/source/lac/constraint_matrix.cc> in
function
void dealii::ConstraintMatrix::close()
The violated condition was:
dof_index != line->line
The name and call sequence of the exception was:
ExcMessage ("Cycle in constraints detected!")
Additional Information:
Cycle in constraints detected!

Stacktrace:
-----------
#0 /home/jhh254/FEM_Softwares/Deall_ii_8.3.0/lib/libdeal_II.g.so.8.3.0:
dealii::ConstraintMatrix::close()
#1 ./step-26: Step26::HeatEquation<3>::setup_system()
#2 ./step-26: Step26::HeatEquation<3>::refine_mesh(unsigned int, unsigned int)
#3 ./step-26: Step26::HeatEquation<3>::run()
#4 ./step-26: main

[Daniel Arndt]

Now it compiles with no problem, but sometimes encounter a run 

time error: " dof_index != line->line" and "Cycle in constraints 

detected!". Do you have any ideas why? It happened at the 35th 

time step, this time, but did not happen last time while running 

with different parameters.

I can't reproduce your run time error. Neither version 8.3.0 nor a current developer version

throw an error when running the program in the attached file. Are you using different parameters?

Best,

Daniel

[jingt...@gmail.com]

Hmm interesting.

For me it happened several times depending on parameters. Once it happened around step 175. Let me try again

Sincerely,

Jingtian Hu

--

[JT Hu]

HI Daniel,


I actually used exactly the same file and got the error at step 35.
Let me try a different computer...


Best,
Jingtian

[JT Hu]

I just tried another computer and had the same problem. There are all
running Ubuntu 14 on virtualbox though. Do you have any suggestions on
troubleshooting?

Thank you so much for helping all the way!

Best,
Jingtian

[Daniel Arndt]

I just rechecked with a clean 8.3.0 installation and can reproduce your run time error.

I would suggest to use a recent developer version if you can. I rechecked that the error vanishes there. 

Best,

Daniel

[JT Hu]

Okay I am trying that! Is there a specific developer version you would
recommend? I am just worrying about stability issues.
Also, did I do anything wrong in the code that caused the error? In
principle, I was not trying to anything that is very fancy, so I did
not expect it to be an internal problem of Deal.II.

Best,
Jingtian

[JT Hu]

Hmm, unfortunately, the output of the 8.4.0 rc3 is probably not
correct. The 8.3 version outputs the expected results before
encountering that error using the same code...
Do you have further suggestions? The model is about heating from
periodic particles.

Best,
Jingtian

[JT Hu]

Hi Daniel,

I figured out some bugs that might cause the problem. The code had too
large a time step at first, and did not constrain the 0 boundary
correctly at first. Now it runs okay with 8.3.0 but not on 8.4.0 rc3.
Do you want the file to figure out the bugs in the pre-built?

Best,
Jingtian

[Daniel Arndt]

Hi Jingtian,

this is indeed strange. Can you localize the error somehow?

Best,

Daniel

[JT Hu]

Hmm, I am okay with debugging, but kind of bad at FEM theory but let
me try changing the mesh size to see if it helps. I will keep you
updated!

Thanks so much for keeping track of this!

Best,
Jingtian