multiply constrained dofs (hanging nodes+periodic) fails a simple test case

[Sambit Das]

Hello dealii-team,

I have created a failing test case in serial when hanging nodes+periodic constraints are mixed. 

The algorithm of the test is follows: (attached minimalExampleHangingNodesPeriodicConstraintsBug.cc)

1) Create a hypercube (-20,20) with origin at the center

2) Set periodic boundary conditions on the hypercube

3) Perform two levels of mesh refinement:- a) one step uniform refinement hypercube to get 8 cells ,and b) then pick one of the 8 cells and refine only that cell which creates hanging nodes on the faces. Finally I get 15 cells (see attached image).

4) Create two constraintMatrices: constraints and onlyHangingNodeConstraints.

                constraints: both hanging node and periodic constraints

                 onlyHangingNodeConstraints: just hanging node constraints

5) Create two nodal vectors vec1 and vec2 with nodal value = nodal_coordinate.norm(). Set and distribute vec1 using  constraints. Set and distribute vec2 using onlyHangingNodeConstraints.

6) Print l2 norm of vec1 and vec2.

The issue:

Since the nodal value = nodal_coordinate.norm() is intrinsically periodic (the domain being a hypercube with origin at the center), I would expect vec1 and vec2 to have the same l2 norm. However I get different l2 norm values:

L2 Norm vec1 (distributed using combined constraints): 178.494

L2 Norm vec2 (distributed using only hanging node constraints): 176.271

If I modify step 3) to: a)  two step uniform refinement hypercube to get 64 cells ,and b) then pick one of the interior cells and refine only that cell which creates hanging nodes only in the interior. Total (71 cells). The l2 norm of vec1 and vec2 matches in this case:

L2 Norm vec1 (distributed using combined constraints): 278.505

L2 Norm vec2 (distributed using only hanging node constraints): 278.505 

Thanks,

Sambit

[Denis Davydov]

Hi Sambit,

On Saturday, April 7, 2018 at 6:51:37 PM UTC+2, Sambit Das wrote:

Hello dealii-team,

I have created a failing test case in serial when hanging nodes+periodic constraints are mixed. 

The algorithm of the test is follows: (attached minimalExampleHangingNodesPeriodicConstraintsBug.cc)

1) Create a hypercube (-20,20) with origin at the center

2) Set periodic boundary conditions on the hypercube

3) Perform two levels of mesh refinement:- a) one step uniform refinement hypercube to get 8 cells ,and b) then pick one of the 8 cells and refine only that cell which creates hanging nodes on the faces. Finally I get 15 cells (see attached image).

4) Create two constraintMatrices: constraints and onlyHangingNodeConstraints.

                constraints: both hanging node and periodic constraints

                 onlyHangingNodeConstraints: just hanging node constraints

5) Create two nodal vectors vec1 and vec2 with nodal value = nodal_coordinate.norm(). Set and distribute vec1 using  constraints. Set and distribute vec2 using onlyHangingNodeConstraints.

6) Print l2 norm of vec1 and vec2.

The issue:

Since the nodal value = nodal_coordinate.norm() is intrinsically periodic (the domain being a hypercube with origin at the center), I would expect vec1 and vec2 to have the same l2 norm. However I get different l2 norm values:

I don't think that's the case. The domain is indeed periodic, but this is completely detached from location of support/nodal points. 

Same applies to geometry, you will have different coordinates of vertices across the PBC so

Point<3> nodalCoor = cellDof->vertex(i);

is expected to give different values on the opposite boundaries.

Regards,

Denis.

[Sambit Das]

Hi Denis,

I don't think that's the case. The domain is indeed periodic, but this is completely detached from location of support/nodal points. 

Same applies to geometry, you will have different coordinates of vertices across the PBC so

I agree, the location of nodal points is detached from the periodicity of the domain, but in this case the origin is at the center of the hypercube. 

This artificially enforces that the nodal_coordinate.norm() is periodic. 

Best,

Sambit

[Denis Davydov]

true, but I don't see why you would have the same norms if you distribute with constraints from hanging nodes only or constraints from hanging nodes+ PBC.

I think we can agree that the two ConstraintMatrix objects should be different as in the case of PBC you additionally need to make sure that FE space on the refined boundary matches that on the opposite, non-refined side. 

If you suspect that there is a bug in constraints, you could check this by simply choosing some more-or-less random vector, distribute and plot-over-line in Paraview / Visit. 

More cumbersome comparison would be to evaluate random field at the opposite points.

You can use FEField function and then choose   L/2-\delta  and -L/2+\delta   with \delta = 1e-8 or so for X coordinate and then 

whatever you want to Y/Z. This should give you the same value anywhere on the two periodically matching points for a random input vector after constraints are distributed.

 

Denis.

[Sambit Das]

true, but I don't see why you would have the same norms if you distribute with constraints from hanging nodes only or constraints from hanging nodes+ PBC.

I think we can agree that the two ConstraintMatrix objects should be different as in the case of PBC you additionally need to make sure that FE space on the refined boundary matches that on the opposite, non-refined side. 

Yes I agree the ConstraintMatrix objects are different but the coefficients (a_{ij}) of the hanging node constraint equations 

x_{i} = a_{ij} x_{j}

 would be the same in both cases, only x_{j}'s would be different in both cases. Now x_{j}'s are nodes without any constraints which are set to the correct values explicitly in both cases:

    if(!constraints.is_constrained(globalDofIndex))

       vec1[globalDofIndex]=nodalCoor.norm();

    if(!onlyHangingNodeConstraints.is_constrained(globalDofIndex))

       vec2[globalDofIndex]=nodalCoor.norm();

So the hanging nodes in both cases should have the same value after calling distribute. 

 

 

If you suspect that there is a bug in constraints, you could check this by simply choosing some more-or-less random vector, distribute and plot-over-line in Paraview / Visit. 

More cumbersome comparison would be to evaluate random field at the opposite points.

You can use FEField function and then choose   L/2-\delta  and -L/2+\delta   with \delta = 1e-8 or so for X coordinate and then 

whatever you want to Y/Z. This should give you the same value anywhere on the two periodically matching points for a random input vector after constraints are distributed.

 I will try doing this. 

Best,

Sambit

[Sambit Das]

I understand my mistake now: there are extra hanging nodes constraints in the constraint matrix with (hanging nod constraints + PBC) compared to the constraint matrix with only hanging node constraints.

That is why the minimal example fails.

Sambit