implemention of periodic boundary condition

[chucu...@gmail.com]

Hi All,

I am having a problem of implementing periodic boundary condition, I have 4 variables to solve, and I make them into a blockvector:

each component of it is a scalar. All of them have periodic boundary condition: 

on domain:

and my code is:

  template <int dim>

  void Problem<dim>::setup_dofs ()

  {

    

    std::vector<unsigned int> block_component (4);

    block_component[0] = 0;

    block_component[1] = 1; 

    block_component[2] = 2;

    block_component[3] = 3;

    

    

    dof_handler.distribute_dofs (fe);

    //DoFRenumbering::component_wise (dof_handler);

    DoFRenumbering::component_wise (dof_handler, block_component);

    

     

      constraints.clear ();

      

      FullMatrix<double> rotation_matrix(dim);

      rotation_matrix[0][0]=1.;

      rotation_matrix[1][1]=1.;

      Tensor<1,dim> offset;

      std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >

      periodicity_vector;

    const unsigned int direction = 0;

           const ComponentMask component_mask = ComponentMask();

           

           std::vector<unsigned int> first_vector_components;

           first_vector_components.push_back(0);

           first_vector_components.push_back(1);

           first_vector_components.push_back(2);

           first_vector_components.push_back(3);

           

           GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 0,

                                  /*b_id2*/ 1,

                                  /*direction*/ 0,

                                  periodicity_vector);

           GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 2,

                                  /*b_id2*/ 3,

                                  /*direction*/ 0,

                                  periodicity_vector);

           DoFTools::make_periodicity_constraints<DoFHandler<dim> >

      (periodicity_vector, constraints);

    constraints.close ();    

    

[...]  

  template <int dim>

  void Problem<dim>::run ()

  {

   

        const unsigned int n_cycles = 1;

      for (unsigned int cycle=0; cycle<n_cycles; ++cycle)

      {

        if (cycle == 0)

          {

            GridGenerator::hyper_cube (triangulation, 0, 2.25);

            

            triangulation.refine_global (7);

            for (typename Triangulation<dim>::active_cell_iterator

                cell = triangulation.begin_active();

                cell != triangulation.end();

                ++cell)

                for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)

                   if (cell->face(f)->at_boundary())

                          if (std::fabs(cell->face(f)->center()(0) - (2.25)) < 1e-12)

                          cell->face(f)->set_boundary_id (1);

                          else if(std::fabs(cell->face(f)->center()(1) - (0)) < 1e-12)

                               cell->face(f)->set_boundary_id (2);

                                else if(std::fabs(cell->face(f)->center()(1) - (2.25)) < 1e-12)

                                     cell->face(f)->set_boundary_id (3);

                                     else if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)

                                          cell->face(f)->set_boundary_id (0);

           std::vector<

           GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator>>

           periodicity_vector;

            FullMatrix<double> rotation_matrix(dim);

            rotation_matrix[0][0]=1.;

            rotation_matrix[1][1]=1.;

           

           const ComponentMask component_mask = ComponentMask();

           

           std::vector<unsigned int> first_vector_components;

           first_vector_components.push_back(0);

           first_vector_components.push_back(1);

           first_vector_components.push_back(2);

           first_vector_components.push_back(3);

                      GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 0,

                                  /*b_id2*/ 1,

                                  /*direction*/ 0,

                                  periodicity_vector);

          GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 2,

                                  /*b_id2*/ 3,

                                  /*direction*/ 0,

                                  periodicity_vector); 

           DoFTools::make_periodicity_constraints<DoFHandler<dim> >

      (periodicity_vector, constraints);
 }

        }

        

        setup_dofs ();
[...]

but when make run it, got the following error:

An error occurred in line <3751> of file <.../boost/deal.ii-8.5.1/src/source/grid/grid_tools.cc> in function

    void dealii::GridTools::collect_periodic_faces(const MeshType&, dealii::types::boundary_id, dealii::types::boundary_id, int, std::vector<dealii::GridTools::PeriodicFacePair<typename MeshType::cell_iterator> >&, const dealii::Tensor<1, MeshType:: space_dimension>&, const dealii::FullMatrix<double>&) [with MeshType = dealii::DoFHandler<2>; dealii::types::boundary_id = unsigned char; typename MeshType::cell_iterator = dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2>, false> >]

The violated condition was: 

    pairs1.size() == pairs2.size()

Additional information: 

    Unmatched faces on periodic boundaries

Stacktrace:

-----------

#0  .../boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1: void dealii::GridTools::collect_periodic_faces<dealii::DoFHandler<2, 2> >(dealii::DoFHandler<2, 2> const&, unsigned char, unsigned char, int, std::vector<dealii::GridTools::PeriodicFacePair<dealii::DoFHandler<2, 2>::cell_iterator>, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::DoFHandler<2, 2>::cell_iterator> > >&, dealii::Tensor<1, dealii::DoFHandler<2, 2>::space_dimension, double> const&, dealii::FullMatrix<double> const&)

#1  ./Problem13-2: Step22::StokesProblem<2>::run()

#2  ./Problem13-2: main

--------------------------------------------------------

make[3]: *** [CMakeFiles/run] Aborted (core dumped)

make[2]: *** [CMakeFiles/run.dir/all] Error 2

make[1]: *** [CMakeFiles/run.dir/rule] Error 2

make: *** [run] Error 2

I wander what is the reason for this issue, and I don't need to implement it in parallel. 

Thanks in advance.

Chucui

[chucu...@gmail.com]

Hi All,

I am working on this project and find some problems new:

1. Why the source code of function "GridTools::collect_periodic_faces" use "cell_iterator" rather than "active_cell_iterator"? As I need to refine my domain by 

triangulation.refine_global (7);

in https://www.mail-archive.com/dea...@googlegroups.com/msg00220.html

I find the same problem, what "Periodic boundaries can only have a difference of 1 refinement level between pairs of faces." means?  And how can I implement my project for refine global (not in parallel)?

2. If I use "Triangulation<dim>" rather than "parallel::distributed::Triangulation<dim>", shall I use the function "GridTools::collect_periodic_faces"? Or just use only "DoFTools::make_periodicity_constraints"?

3. If use "DoFTools::make_periodicity_constraints" only, my code is

  template <int dim>

  void StokesProblem<dim>::setup_dofs ()

  {

    

    std::vector<unsigned int> block_component (4);

    block_component[0] = 0;

    block_component[1] = 1; 

    block_component[2] = 2;

    block_component[3] = 3;

    

    

    dof_handler.distribute_dofs (fe);

    //DoFRenumbering::component_wise (dof_handler);

    DoFRenumbering::component_wise (dof_handler, block_component);

    

     

    constraints.clear ();

    const unsigned int n_dofs_per_face = fe.dofs_per_face;

    FullMatrix<double> rotation_matrix(n_dofs_per_face);

    for (unsigned int d=0; d<n_dofs_per_face; ++d)

         {

           rotation_matrix[d][d]=1.;

         }

    Tensor<1,dim> offset;

    const ComponentMask component_mask = ComponentMask();

           

           std::vector<unsigned int> first_vector_components;

           first_vector_components.push_back(0);

           first_vector_components.push_back(1);

           first_vector_components.push_back(2);

           first_vector_components.push_back(3);

           

             const typename DoFHandler<dim>::active_face_iterator face_0; 

             const typename DoFHandler<dim>::active_face_iterator face_1;

             const typename DoFHandler<dim>::active_face_iterator face_2;

             const typename DoFHandler<dim>::active_face_iterator face_3;

     for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();

          cell != dof_handler.end();

          ++cell)

       {     

         for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)

           {

             const typename DoFHandler<dim>::active_face_iterator face = cell->face(i);

             

             if (face->at_boundary() && face->boundary_id() == 0)

               {

                const typename DoFHandler<dim>::active_face_iterator face_0 = face; 

               }

 

             else if (face->at_boundary() && face->boundary_id() == 1)

               {

                 const typename DoFHandler<dim>::active_face_iterator face_1 = face; 

               }

             else if (face->at_boundary() && face->boundary_id() == 2)

               {

                const typename DoFHandler<dim>::active_face_iterator face_2 = face; 

               }

 

            else if (face->at_boundary() && face->boundary_id() == 3)

               {

                 const typename DoFHandler<dim>::active_face_iterator face_3 = face; 

               }   

               

           }

         

             DoFTools::make_periodicity_constraints (face_0, face_1, constraints, component_mask, true, false, false, rotation_matrix, first_vector_components);

             DoFTools::make_periodicity_constraints (face_2, face_3, constraints, component_mask, true, false, false, rotation_matrix, first_vector_components);

       }
[...]

but get the error:

Linking CXX executable Problem13-2

/.../13-4-pbc-for-run/together-pbc.cc:647: error: undefined reference to 'void dealii::DoFTools::make_periodicity_constraints<dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > >(dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > const&, dealii::identity<dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > >::type const&, dealii::ConstraintMatrix&, dealii::ComponentMask const&, bool, bool, bool, dealii::FullMatrix<double> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&)'

/.../13-4-pbc-for-run/together-pbc.cc:648: error: undefined reference to 'void dealii::DoFTools::make_periodicity_constraints<dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > >(dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > const&, dealii::identity<dealii::TriaActiveIterator<dealii::DoFAccessor<1, dealii::DoFHandler<2, 2>, false> > >::type const&, dealii::ConstraintMatrix&, dealii::ComponentMask const&, bool, bool, bool, dealii::FullMatrix<double> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&)'

collect2: error: ld returned 1 exit status

make[3]: *** [Problem13-2] Error 1

make[2]: *** [CMakeFiles/Problem13-2.dir/all] Error 2

make[1]: *** [CMakeFiles/run.dir/rule] Error 2

make: *** [run] Error 2

I wander what is the reason for this issue.

Thanks very much!

Best,

Chucui

在 2018年8月5日星期日 UTC+8下午3:28:01，chucu...@gmail.com写道：

[Jean-Paul Pelteret]

Dear Chucui,

Since you haven’t provided a minimal working example, I could only speculate what might be the issue. From the incomplete example that you’ve shown here, I would guess that this might be the problem:

On 05 Aug 2018, at 09:28, chucu...@gmail.com wrote:

                      GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 0,

                                  /*b_id2*/ 1,

                                  /*direction*/ 0,

                                  periodicity_vector);

          GridTools::collect_periodic_faces(dof_handler,

                                  /*b_id1*/ 2,

                                  /*b_id2*/ 3,

                                  /*direction*/ 0,   <——— Should be a "1"

                                  periodicity_vector);

For the unit cube coarse face pairs {0,1} should have their normal in the x-direction, and coarse face pairs {2,3} would then have their normal vectors orientated in the y-direction.

If not then it *might* be that your coarsest grid needs to have more than one cell. However, looking at the documentation to step-45 (specifically the introduction), I don’t believe this to be the case. You can try to take the coarsest grid, refine it globally once, and then use the GridTools::flatten_triangulation() function to make this the new coarsest grid.

I hope that this helps you solve the problem.

Best,

Jean-Paul

[Jean-Paul Pelteret]

Hi Chucui,

I’ll try to answer some of your questions quickly:

1. Why the source code of function "GridTools::collect_periodic_faces" use "cell_iterator" rather than "active_cell_iterator”?

Because your mesh should be periodic from the coarsest representation all the way down to the finest representation. This is important for very large problems, as deal.II can then exploit the parent-child relationship to find matching DoF pairs rather than traversing the entire fine-scale problem to do the same.

As I need to refine my domain by 

triangulation.refine_global (7);

in https://www.mail-archive.com/dea...@googlegroups.com/msg00220.html

I find the same problem, what "Periodic boundaries can only have a difference of 1 refinement level between pairs of faces." means?  

Consider a mesh on which boundary with id 0 is considered the periodic match to a boundary with id 1. This means that there are faces that have boundary id 0 that are considered the periodic neighbours to faces with boundary id 1. In the simplest view of things, one might wish to ensure that refinement of cells along boundary 0 is exactly the same as that along boundary 1, and the number of DoFs along boundary 0 exactly matches boundary 1. However, it is possible to implement constraints for “hanging nodes” between boundaries 0 and 1, i.e. the refinement of cells along boundary 0 is not the same as that along boundary 1, and the number of DoFs does not necessarily match either. To accomplish this, the same arguments for normal internal hanging nodes applies, and therefore periodic neighbours can only have 1 level of refinement difference between them.

And how can I implement my project for refine global (not in parallel)?

There should be no problem with this. I would suggest looking at step-45 if you haven’t done so already.

2. If I use "Triangulation<dim>" rather than "parallel::distributed::Triangulation<dim>", shall I use the function "GridTools::collect_periodic_faces"? Or just use only "DoFTools::make_periodicity_constraints"?

I guess that depends on what you’re trying to achieve… Based on your first query with your problem statement, DoFTools::make_periodicity_constraints() should do what you are wanting.

I don’t feel too inclined to look into this in any detail, but it could be that DoFTools::make_periodicity_constraints() is only instantiated for a dealii::TriaIterator rather than a dealii::TriaActiveIterator. I guess that you would need to pass it a typename DoFHandler<dim>::face_iterator.

I hope that I’ve been able to give you a little useful insight into what you’ve asked here.

Best,

Jean-Paul

[chucu...@gmail.com]

Dear Jean-Paul,

Thank you very much! I had read step-45, but still have some questions:

1. How to use function "DoFTools::make_periodicity_constraints() "  only? if my code need not to compute parallel case.  If I choose

void DoFTools::make_periodicity_constraints	(	const FaceIterator &	face_1,
		const typename identity< FaceIterator >::type &	face_2,
		AffineConstraints< double > &	constraints,
		const ComponentMask &	component_mask = ComponentMask(),
		const bool	face_orientation = true,
		const bool	face_flip = false,
		const bool	face_rotation = false,
		const FullMatrix< double > &	matrix = FullMatrix<double>(),
		const std::vector< unsigned int > &	first_vector_components = std::vector<unsigned int>()
	)

how to choose face1 and face2 we need there correctly? If I use 

template<typename DoFHandlerType >

void DoFTools::make_periodicity_constraints	(	const std::vector< GridTools::PeriodicFacePair< typename DoFHandlerType::cell_iterator >> &	periodic_faces,
		AffineConstraints< double > &	constraints,
		const ComponentMask &	component_mask = ComponentMask(),
		const std::vector< unsigned int > &	first_vector_components = std::vector<unsigned int>()
	)

how to make a 

const std::vector< GridTools::PeriodicFacePair< typename DoFHandlerType::cell_iterator >> &

periodic_faces

maybe I still need to use "GridTools::collect_periodic_faces" to input information of faces into PeriodicFacePair. If I do this,  I am not use "DoFTools::make_periodicity_constraints() "  only....

2. In step-45, the complete implemention is parallel case. For every step, it use functions used in parallel case. If I want to implement a case not in parallel, shall I rewrite all of them and remove the parallel parts (still not succeed), or use it directly, but in shell run my program via

mpirun -np 1 ./myCase

Thank you very much!

Best,

Chucui

[Daniel Arndt]

Chucui,

1. How to use function "DoFTools::make_periodicity_constraints() "  only? if my code need not to compute parallel case.  If I choose [...] how to choose face1 and face2 we need there correctly? If I use [...] maybe I still need to use "GridTools::collect_periodic_faces" to input information of faces into PeriodicFacePair. If I do this,  I am not use "DoFTools::make_periodicity_constraints() "  only....

In your case, it should be sufficient to use the overload that just takes boundary ids. Alternatively, you can use GridTools::collect_periodic_faces and pass the return value to the respective overload of GridTools::make_periodicity_constraints.

 

2. In step-45, the complete implemention is parallel case. For every step, it use functions used in parallel case. If I want to implement a case not in parallel, shall I rewrite all of them and remove the parallel parts (still not succeed), or use it directly, but in shell run my program via

mpirun -np 1 ./myCase

 

As explained in the introduction of step-45 there are only very few extra steps you need to do when you are dealing with a parallel::distributed::Triangulation object.

In the end, it is sufficient to call GridTools::make_periodicity_constraints (in case you are using continuous ansatz spaces).

Note that you have to take care yourself that there is at most a one level difference between periodic faces.

Best,

Daniel

[chucu...@gmail.com]

Dear Daniel,

Thank you very much for your rapid reply! I use overload as you say, and it can run before assemble_system, the relevant part of my code is

  template <int dim>

  void Problem<dim>::run ()

  {

   

        const unsigned int n_cycles = 1;

        double entropy = 0;

        

    for (unsigned int cycle=0; cycle<n_cycles; ++cycle)

      {

        if (cycle == 0)

          {

            GridGenerator::hyper_cube (triangulation_active, 0, 2.25);

            triangulation_active.refine_global (7);

            GridGenerator::flatten_triangulation (triangulation_active, triangulation);

            for (typename Triangulation<dim>::active_cell_iterator

                cell = triangulation.begin_active();

                cell != triangulation.end();

                ++cell)

            {    

                for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)

                {

                   if (cell->face(f)->at_boundary())

                   {                      

                          if (std::fabs(cell->face(f)->center()(0) - (2.25)) < 1e-12)

                          

                          cell->face(f)->set_boundary_id (1);

                          else if(std::fabs(cell->face(f)->center()(1) - (0)) < 1e-12)

                               

                               cell->face(f)->set_boundary_id (2);

                               else if(std::fabs(cell->face(f)->center()(1) - (2.25)) < 1e-12)

                                     

                                     cell->face(f)->set_boundary_id (3);

                                     else if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)

                                          {

                                          cell->face(f)->set_boundary_id (0);

                                          }  

                   }                     

               }

            }

             

         }

        else

        {

        //triangulation.refine_global (1);

        }

           setup_dofs ();

           std::cout << "start" << cycle << std::endl;  

                               

        VectorTools::project (dof_handler_phi, constraints, QGauss<dim>(degree+2),

                        InitialValues_phi<dim>(),

                        phi_0);                       

        VectorTools::project (dof_handler_e, constraints, QGauss<dim>(degree+2),

                        InitialValues_e<dim>(),

                        e_0);  

        std::cout << "phi_0" << cycle << std::endl;

        //VectorTools::project (dof_handler_U, constraints_U, QGauss<dim>(degree+2),

        //                InitialValues_U<dim>(),

        //                U_n);                        

        //std::cout << "U_0" << cycle << std::endl;  

        VectorTools::project (dof_handler_q, constraints_q, QGauss<dim>(degree+2),

                        InitialValues_q<dim>(),

                        q_0);                  

        std::cout << "q_0" << cycle << std::endl;                 

        assemble_system_1 (phi_0, e_0);
[...]

and my setup_dofs ():

  template <int dim>

  void Problem<dim>::setup_dofs ()

  {

    std::vector<unsigned int> block_component (4);

    block_component[0] = 0;

    block_component[1] = 1; 

    block_component[2] = 2;

    block_component[3] = 3;

    

    dof_handler.distribute_dofs (fe);

    //DoFRenumbering::component_wise (dof_handler);

    DoFRenumbering::component_wise (dof_handler, block_component);

    constraints.clear ();

    const ComponentMask component_mask = ComponentMask();

    DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, component_mask);

    DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, component_mask);

    constraints.close (); 
    

    std::vector<types::global_dof_index> dofs_per_block (4); 

    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);

    

    const unsigned int n_phi = dofs_per_block[0],

                       n_e = dofs_per_block[1],

                       n_U = dofs_per_block[2],

                       n_q = dofs_per_block[3];

    system_matrix_1.clear ();

   

    BlockDynamicSparsityPattern dsp (4,4);

    dsp.block(0,0).reinit (n_phi, n_phi);

    dsp.block(1,0).reinit (n_e, n_phi);  

    dsp.block(2,0).reinit (n_U, n_phi);

    dsp.block(3,0).reinit (n_q, n_phi); 

     

    dsp.block(0,1).reinit (n_phi, n_e);

    dsp.block(1,1).reinit (n_e, n_e);

    dsp.block(2,1).reinit (n_U, n_e);

    dsp.block(3,1).reinit (n_q, n_e); 

    

    dsp.block(0,2).reinit (n_phi, n_U);

    dsp.block(1,2).reinit (n_e, n_U);  

    dsp.block(2,2).reinit (n_U, n_U);

    dsp.block(3,2).reinit (n_q, n_U); 

     

    dsp.block(0,3).reinit (n_phi, n_q);

    dsp.block(1,3).reinit (n_e, n_q);

    dsp.block(2,3).reinit (n_U, n_q);

    dsp.block(3,3).reinit (n_q, n_q); 

    

    dsp.collect_sizes();

    DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);

    sparsity_pattern.copy_from (dsp);

    

    system_matrix_1.reinit (sparsity_pattern);

[...]

the output and error is

start0

phi_00

q_00

--------------------------------------------------------

An error occurred in line <1764> of file </.../boost/deal.ii-8.5.1/deal.II/include/deal.II/lac/sparse_matrix.h> in function

    void dealii::SparseMatrix<number>::add(dealii::SparseMatrix<number>::size_type, dealii::SparseMatrix<number>::size_type, number) [with number = double; dealii::SparseMatrix<number>::size_type = unsigned int]

The violated condition was: 

    (index != SparsityPattern::invalid_entry) || (value == number())

Additional information: 

    You are trying to access the matrix entry with index <0,1>, but this entry does not exist in the sparsity pattern of this matrix.

The most common cause for this problem is that you used a method to build the sparsity pattern that did not (completely) take into account all of the entries you will later try to write into. An example would be building a sparsity pattern that does not include the entries you will write into due to constraints on degrees of freedom such as hanging nodes or periodic boundary conditions. In such cases, building the sparsity pattern will succeed, but you will get errors such as the current one at one point or other when trying to write into the entries of the matrix.

I wander how to build the sparsity pattern completely (when I use the same part as step-45, i.e. the same setup-dof(), but receive the same error as above), maybe I need to add some information to the " BlockDynamicSparsityPattern dsp", but how to do it? I can't find a correct example to describe this.

[Daniel Arndt]

Chucui,

In first sight this looks reasonable. In particular, it is important to consider the constraints when building the sparsity pattern.

Can you build a minimal working example that shows the problem?

Best,

Daniel

[Jean-Paul Pelteret]

Dear Chucui,

I think that what Daniel was asking for was the absolute minimum size test case that reproduces the error in building the sparsity pattern. I interpreted (and he’s welcome to correct me) that wants to help you find out why you hit this error:

An error occurred in line <1764> of file </.../boost/deal.ii-8.5.1/deal.II/include/deal.II/lac/sparse_matrix.h> in function
    void dealii::SparseMatrix<number>::add(dealii::SparseMatrix<number>::size_type, dealii::SparseMatrix<number>::size_type, number) [with number = double; dealii::SparseMatrix<number>::size_type = unsigned int]
The violated condition was: 
    (index != SparsityPattern::invalid_entry) || (value == number())
Additional information: 
    You are trying to access the matrix entry with index <0,1>, but this entry does not exist in the sparsity pattern of this matrix.

So you can remove a lot of the code that would only be used after the sparsity pattern was successfully built. That would be a great help and would accelerate the process of getting an answer to your question.

Best,

Jean-Paul

On 08 Aug 2018, at 15:05, chucu...@gmail.com wrote:

Chucui

[Jean-Paul Pelteret]

Dear Chucui,

Thanks for giving this a go. This is better than before, but I think that you still misinterpret what a minimal working example is. This is a 1000 line piece of code, but I estimate only 100 lines of code are required to reproduce the problem. If the problem arises in the first call to setup_dofs(), which to the best of my understanding it does, then you really don’t need anything other than

- the set of calls used to create the Triangulation, and 

- the call to setup_dofs() itself. 

So could you please remove everything below the call to setup_dofs() in run(), and everything related to the assembly, solving the linear system, exact solution, etc. Forget about solving the original problem -- really try to make this compact while reproducing the error message hit when creating the sparse matrix. Then we will have a focussed piece of code that we can use to debug the issue, and if there is something to be changed in the library we will add this code to the test suite and would naturally give you attribution for it.

Best,

Jean-Paul

P.S. A huge plus would be to strip all of the remaining functional code out of the Problem class and have it all in the run() function. This would be nice, but is not entirely necessary.

On 08 Aug 2018, at 15:40, chucu...@gmail.com wrote:

Dear Daniel,

I'm sorry to make this mistake, (and thank Jean-Paul very much). The newest code focus on build the sparsity pattern in "setup_dof()", while other parts (assemble_system(), assemble_rhs(), run() ) need to be left to make the code work and maintain the structure of 3 components of the blockvector solution.

Thank you much!

Best,

Chucui

--
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.
<LittleCase0-3.cc>

[chucu...@gmail.com]

Dear Jean-Paul,

Thank you very much for your patience, when I print some information like this:

template <int dim>

void Problem<dim>::run ()

{

[...]
  setup_dofs ();

std::cout << "test-1 "<< std::endl;                              

        VectorTools::project (dof_handler_0, constraints_phi, QGauss<dim>(degree+2),

                        SolutionsPhi<dim>(),

                        phi_0);      

std::cout << "test-2 " << std::endl;                                          

        VectorTools::project (dof_handler_2, constraints_q, QGauss<dim>(degree+2),

                        SolutionsQ<dim>(),

                        q_0);  

 std::cout << "test-3 " << std::endl;                           

        assemble_system_1 (phi_0);    

 std::cout << "test-4" << std::endl;                   

        assemble_rhs_1 (phi_0, q_0);

 std::cout << "test-5" << std::endl;         

and the output is

Linking CXX executable Case0-3

[ 50%] Built target Case0-3

[100%] Run with Debug configuration

   Number of active cells: 16384

   Number of degrees of freedom: 98818 (16641+16641+65536)

test-1 

test-2 

test-3 

--------------------------------------------------------

An error occurred in line <1764> of file </home/yucan/boost/deal.ii-8.5.1/deal.II/include/deal.II/lac/sparse_matrix.h> in function

    void dealii::SparseMatrix<number>::add(dealii::SparseMatrix<number>::size_type, dealii::SparseMatrix<number>::size_type, number) [with number = double; dealii::SparseMatrix<number>::size_type = unsigned int]

The violated condition was: 

    (index != SparsityPattern::invalid_entry) || (value == number())

Additional information: 

    You are trying to access the matrix entry with index <0,1>, but this entry does not exist in the sparsity pattern of this matrix.

So it can't run when arrive at "assemble_system_1", so "setup-dof()",  "assemble_system_1()", "run()" are left—— this is LittleCase0-3.cc

if only "setup_dof()", "run()" are left——this is LittleCase0-4.cc

Hope this is a minimal example. Thank you very much!

Best,

Chucui

[Daniel Arndt]

Chucui,

it is clearly suspicious that you don't use any of your ConstraintMatrix objects in the assembly loops.

Instead of calling

system_matrix.add(local_dof_indices[i],
                             local_dof_indices[j],
                             local_matrix(i, j));

you should use

constraints.distribute_local_to_global(local_matrix, local_dof_indices, system_matrix)

In case that you also have inhomogeneous boundary conditions, this is in general not sufficient

and you should asemble the respective right hand sides in the same loop you construct the matrix.

Then, call

constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);

The error you are getting is only slightly related to using periodic boundary conditions. By saying

DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);

you are telling the sparsity pattern (or the matrix) that you will never write to the constraint entries.

This assumption is true, when you use constraints.distribute_local_to_global, but not when you copy the matrix entries manually

(without taking the constraints into account at all).

Best,

Daniel

[Jean-Paul Pelteret]

Dear Chucui,

Thanks for the explanation - clearly I was wrong with regards to exactly how much complexity could be stripped from this example. Anyway, I think that Daniel has made the critical observation as to where the fault may lie. Hopefully his suggestion will sort this out for you.

Best,

Jean-Paul

--
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.

<LittleCase0-3.cc><LittleCase0-4.cc>

[chucu...@gmail.com]

Dear Daniel,

Thank you very much! I use code like this:
constraints.distribute_local_to_global(local_matrix, local_dof_indices, system_matrix) as you say, 

and use codes like this

constraints.distribute(solution);

after every "solve()" step.

and now, my LittleCase0-3.cc can run. 

Maybe a question more : how to make sure your period boundary coditions are added truly?

[Daniel Arndt]

Chucui,

Maybe a question more : how to make sure your period boundary coditions are added truly?

How would you do it for Dirichlet boundary conditions?

There are mutiple possibilities.

- You can just print the constraints (ConstraintMatrx.print(std::cout)),

print the locations of degrees of freedom via DoFTools::write_gnuplot_dof_support_point_info (https://www.dealii.org/9.0.0/doxygen/deal.II/namespaceDoFTools.html#a69d19d6d574269cc6e69fa5c5b2d89e2)

and compare.

- You can visually confirm that your solution looks periodic.

- You can calculate the solution for all the support points at the periodic faces and compare the values.

Best,

Daniel

[chucu...@gmail.com]

Dear Daniel,

Thank you very much! I can test my period boundary condition added truly. Your advision helps me a lot !

As my question is about periodic boundary condition like 

just like Dirichlet boundary condition.

But if "periodic boundary condition" means that all information of solution on periodic faces pair are equal, not only the value of solution (like Dirichlet boundary condition above), but also gradient value of solution (like Neuman boundary condition, as normal vector is opposite, so if grad_u_on_left = grad_u_on_right, for example, so n_left * grad_u_on_left + n_right * grad_u_on_right, top and bottom as well) . In my case, periodic boundary condition should include

Does deal.ii implement this kind of periodic boundary condition in fact? Is it different from what we have done above for periodic boundary condition like Dirichlet?

[Daniel Arndt]

Chucui,

But if "periodic boundary condition" means that all information of solution on periodic faces pair are equal, not only the value of solution (like Dirichlet boundary condition above), but also gradient value of solution (like Neuman boundary condition, as normal vector is opposite, so if grad_u_on_left = grad_u_on_right, for example, so n_left * grad_u_on_left + n_right * grad_u_on_right, top and bottom as well) . [...]

Does deal.ii implement this kind of periodic boundary condition in fact? Is it different from what we have done above for periodic boundary condition like Dirichlet?

Short answer: No.

Longer answer:

We identify the degrees of freedom on the respective faces with eachother. This implies treating them the same as internal faces.

In case, you have a finite element that is C1-conforming this means that the ansatz space will also be C1-conforming across periodic boundary faces.

However, for the standard C0-conforming FE_Q elements this only implies continuity with respect to values.

In case you need continuity of the gradients across periodic faces, I would expect that you require continuity of the gradients across (internal) faces as well.

Is that the case?

Best,

Daniel

and need to choose a

[Wolfgang Bangerth]

> But if "periodic boundary condition" means that all information of solution on
> periodic faces pair are equal, not only the value of solution (like Dirichlet
> boundary condition above), but also gradient value of solution (like Neuman

> boundary condition, as normal vector is opposite, so ifgrad_u_on_left =
> grad_u_on_right, for example, so n_left *grad_u_on_left + n_right *

> grad_u_on_right, top and bottom as well) . In my case, periodic boundary
> condition should include
>

> 46.JPG <about:invalid#zClosurez>

Do you mean that you want *both* Dirichlet and Neumann conditions to be
satisfied? Or that you don't want Dirichlet periodic boundary conditions and
only want the gradient to be periodic?

Daniel already answered the question from the *discrete* point of view. But
from the perspective of the PDE, if you have an elliptic PDE with smooth
solutions in the interior of the domain, then the solution of the PDE with
periodic boundary conditions is also smooth across the boundary. That means
that by requiring that the *value* of the solution is periodic across the
boundary, you also get that the *gradient* of the solution is periodic across
the boundary. This does not have to explicitly prescribed: It just happens, in
the same way as it is the case in the interior of the domain.

Of course, since we only approximate these solutions with piecewise
polynomials, continuity of value does not imply continuity of the gradient
across periodic boundaries. But this is the same as between faces of the mesh:
the solution is continuous across a face, but the gradient of the finite
element approximation is not.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[chucu...@gmail.com]

Dear Daniel and Wolfgang,

Thank you very much! In fact, I only want the gradient of periodic, not the Dirichlet periodic boundary. I thought they are same in periodic boundary conditions before, i.e. all informations of the solution on the face pairs are same. Now, I get that I have made mistake.

Why I need the gradient of periodic only? Because in my case, I have some boundary integrals in my cost functional, I need to make all of them equal to zero, so I can choose the gradient of periodic only.

But in the analysis of weak solution of original PDEs in my case, the space of trial and test funcution are (H^1) * (H^1) * (L^2) * (H^1) (sobolev space), as there are the gradient of 1st, 2nd, 4th component in the weak form. And I choose C0 * C0 * C0 * C0 FEM space

So I cannot make sure the continuity of the gradient in boundary or internal faces. But I don' know how to implement what I need? Maybe some penalty part like in Discontinuous Gradient Method to make up for continuity of the gradient?

[Wolfgang Bangerth]

On 08/12/2018 12:30 AM, chucu...@gmail.com wrote:
>
> So I cannot make sure the continuity of the gradient in boundary or internal
> faces. But I don' know how to implement what I need? Maybe some penalty part
> like in Discontinuous Gradient Method to make up for continuity of the gradient?

Yes, that's what I would do -- add a penalty term for the discontinuity of the
gradient across the periodic boundary.