Hanging node constraints

[Deepak Gupta]

Dear All,  

Below is a simple piece of code where one finite element has a different p-order compared to the rest. Thus, I expect certain hanging_node_constraints (which is 4) in the output. What I cannot figure out is why I get zero constraints when I check the hanging constraints a second time? More can be understood from the code below which forms part of the setup_system() function in my code. The associated output follows the code.

 boundary_values.clear();
    dof_handler.distribute_dofs(fe_collection);
    hanging_node_constraints.clear();
    DoFTools::make_hanging_node_constraints(dof_handler,
            hanging_node_constraints);
    boundary_info();  //created by  me for updating rhe boundary indicators
    //Applying the boundary conditions
    BoundaryValues<dim> boundary_v;
    VectorTools::interpolate_boundary_values(dof_handler,
            42,
            boundary_v,
            boundary_values);
    //hanging_node_constraints.close();
      std::cout<< "No.of degrees of freedom: " << dof_handler.n_dofs() << "\n";
    std::cout<<"No. of hanging node constraints : "<<hanging_node_constraints.n_constraints()<<std::endl;



    boundary_values.clear();
    hanging_node_constraints.clear();
    DoFTools::make_hanging_node_constraints(dof_handler,
            hanging_node_constraints);
    boundary_info(); //created by  me for updating rhe boundary indicators
    //Applying the boundary conditions
   VectorTools::interpolate_boundary_values(dof_handler,
            42,
            boundary_v,
            boundary_values);
    hanging_node_constraints.close();
      std::cout<< "No.of degrees of freedom: " << dof_handler.n_dofs() << "\n";
    std::cout<<"No. of hanging node constraints : "<<hanging_node_constraints.n_constraints()<<std::endl;

OUTPUT
No.of degrees of freedom: 6652
No. of hanging node constraints : 4
No.of degrees of freedom: 6652
No. of hanging node constraints : 0

--

Deepak K. Gupta

[Wolfgang Bangerth]

On 09/27/2016 06:44 AM, Deepak Gupta wrote:
> Dear All,
>
> Below is a simple piece of code where one finite element has a different
> p-order compared to the rest. Thus, I expect certain hanging_node_constraints

> (which is 4) in the output. What I cannot figure out is why I get *zero
> *constraints when I check the hanging constraints a second time? More can be

> *OUTPUT

> No.of degrees of freedom: 6652
> No. of hanging node constraints : 4
> No.of degrees of freedom: 6652
> No. of hanging node constraints : 0

Deepak, this looks wrong but since this is such a well tested part of the
library, I can only assume that your code is doing something wrong. But then,
I've been known to suspect the bug somewhere outside the library :-) Can you
create a small, self-contained testcase that shows the problem?

Best
W.


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

[Deepak Gupta]

Hi Wolfgang,

I created a simple example but it worked correctly and is attached (I modified a test code given earlier my Jean-Paul). The difference between this simple code and my code stated earlier is that between the two checks of hanging nodes, the interpolate_boundary_values function is not called in hp_check.cpp. However, in the piece of code sent earlier, there is this line:

   VectorTools::interpolate_boundary_values(dof_handler,
            42,
            boundary_v,
            boundary_values);

This is somehow causing the problem. If I comment this one, the hanging constraints are correctly. I am trying to put this as well in my simplified example. In the meantime, may you be you have an idea why this is messing up. The online documentation is down, so could not check details of the function.

Best regards

Deepak

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 09/27/2016 08:52 AM, Deepak Gupta wrote:
>
> I created a simple example but it worked correctly and is attached (I modified
> a test code given earlier my Jean-Paul). The difference between this simple
> code and my code stated earlier is that between the two checks of hanging
> nodes, the interpolate_boundary_values function is not called in hp_check.cpp.
> However, in the piece of code sent earlier, there is this line:
>
> VectorTools::interpolate_boundary_values(dof_handler,
> 42,
> boundary_v,
> boundary_values);
>
> This is somehow causing the problem. If I comment this one, the hanging
> constraints are correctly. I am trying to put this as well in my simplified
> example. In the meantime, may you be you have an idea why this is messing up.
> The online documentation is down, so could not check details of the function.

Well, VectorTools::interpolate_boundary_values() doesn't touch the
'constraints' variable, so it really shouldn't have any influence.

When trying to find issues such as this, I usually try to start with something
that doesn't work, and remove stuff as far as I can. So, start with your
non-working program and make it smaller and smaller until the issue goes away.
The part that you removed last was then the offending statement. Put it back
in, and that is your minimal testcase.

[Deepak Gupta]

Thanks Wolfgang,

I did as suggested and figured out the error. Although, I have also solved it, I would like to know why this was a problem. For the same, I have attached the new simple version. In that code, there is a piece of sub code added which is actually the content of boundary_info() function which existed earlier.

    for (typename Triangulation<dim>::active_cell_iterator
            cell = tria.begin();
            cell != tria.end();
            ++cell ){
        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f){
            Point<dim> point1 = cell->face(f)->center();
            std::vector<double> x(dim);
            for(unsigned int  i = 0; i < dim; ++i){
                x[i] = point1(i);
            }
            if (fabs(point1(0) - 0.0) < 1e-12)
                cell->face(f)->set_boundary_indicator(42);
            else
                cell->face(f)->set_boundary_indicator(0);

        }
    }

If I remove the else part of the if-condition, it works, else the constraints become zero. This implies, I should not try to set boundary indicator for anything which is not at the boundary, which is obvious. I know it was a mistake from my side, but I would like to know why it was a problem (what goes on in the background as per the implementation?)

Best

Deepak

[Wolfgang Bangerth]

On 09/27/2016 10:05 AM, Deepak Gupta wrote:
>
>
> If I remove the else part of the if-condition, it works, else the constraints
> become zero. This implies, I should not try to set boundary indicator for
> anything which is not at the boundary,

Correct.

> which is obvious. I know it was a
> mistake from my side, but I would like to know why it was a problem (what goes
> on in the background as per the implementation?)

In essence you are saying that some interior face is a boundary. If you come
from each of the two adjacent cells, the interface will then look like an
external boundary (at which there can be no neighbor) and so no constraints
will be computed.

[Deepak Gupta]

Thanks, after this explanation, it sounds extremely silly on my side :)

Best,

Deepak

[Wolfgang Bangerth]

On 09/27/2016 10:12 AM, Deepak Gupta wrote:
> Thanks, after this explanation, it sounds extremely silly on my side :)

You're not the first one to make this mistake.

[Timo Heister]

>> Thanks, after this explanation, it sounds extremely silly on my side :)
>
> You're not the first one to make this mistake.

That raises the question if there is ever a valid reason to call
set_boundary_indicator on an interior face or if we should put an
Assert() into place.

--
Timo Heister
http://www.math.clemson.edu/~heister/

[Wolfgang Bangerth]

On 09/27/2016 11:18 AM, Timo Heister wrote:
> That raises the question if there is ever a valid reason to call
> set_boundary_indicator on an interior face

No.

> or if we should put an
> Assert() into place.

That's difficult because a face itself doesn't know whether it's at the
boundary or not. Only cells know whether they're at the boundary. In other
words, once you write
cell->face(f)->set_boundary(42)
you've lost the information at the first ->. We would need to add a flag to
all faces that says whether they're at the boundary. We could then also write a
face->at_boundary()
function, and do
Assert (!at_boundary(), ...)
in the implementation of set_boundary. But we will have to have that flag first.

Feel free to copy this into a github issue if you think it would be a
worthwhile project.

Best