Problems with Constraintmatrix and hanging nodes on corners in 3D

[Bernhard Endtmayer]

Dear deal.II Developer Team,

I' am trying to solve the stationary Navier-Stokes benchmark Problem with configuration 2 given in

https://ganymed.math.uni-heidelberg.de/Paper/braackrichter-200410.pdf

The initial mesh is attached.

Here i used to following code fragment to set the boundary conditions for the newton update:

****************************************************************************************************************

template<int dim>
void set_newton_bc(DoFHandler<dim>& dof_handler, ConstraintMatrix & constraints, Dirichlet_Boundary dirichlet_doundary_enforcement =
                         Dirichlet_Boundary::strong)
{

   std::vector<bool> component_mask(dof_handler.get_fe(0).n_components(), true);
   component_mask[dim] = false;

   VectorTools::interpolate_boundary_values(dof_handler, 0, ZeroFunction<dim>(dof_handler.get_fe(0).n_components()), constraints,
                                            component_mask);
   VectorTools::interpolate_boundary_values(dof_handler, 2, ZeroFunction<dim>(dof_handler.get_fe(0).n_components()), constraints,
                                            component_mask);



   VectorTools::interpolate_boundary_values(dof_handler, 80, ZeroFunction<dim>(dof_handler.get_fe(0).n_components()), constraints,
                                            component_mask);
   IndexSet indexset;
   std::set<unsigned int> boundaryids;
   boundaryids.clear();
   boundaryids.insert(80);
   DoFTools::extract_boundary_dofs(dof_handler,ComponentMask(component_mask),indexset,boundaryids);
   indexset.print(std::cout);
}

****************************************************************************************************************

The last part of the code fragment is just to get an output of the index set of the DoF which belong to boundary 80.

The set is given by

{31, [34,36], 90, [162,163], 202, 208, 252, [255,257], 277, 279, [281,282], 286, 288, 291, 305, 307, 309, [327,328], [331,332], [343,344], [356,357], 399, [406,409], 417, 419, 430, [433,434], 438, 441, 458, 463, [469,470], [473,474], [505,506], 508, [511,512], 524, [544,545], [551,552], [556,557], [564,565], 586, 590, 593, 605, 611, 625, 636, 640, 642, [655,656], 658, [660,661], 665, 701, 703, 705, 707, 709, 712, 714, 716, 718, 720, [738,739], 743, 745, 747, 777, 779, [781,782], 786, 816, 818, 820, [822,824], [880,881], [884,886], 891, [894,895], [901,903], [930,931], [933,934], 942, [974,975], 1001, 1006, [1015,1017], 1021, [1025,1026], 1049, 1051, 1064, 1085, 1093, [1099,1102], [1115,1116], 1130, [1143,1144], [1154,1155], 1159, 1161, 1163, [1191,1192], 1196, 1198, 1200, 1205, 1212, 1214, 1228, [1233,1234], 1245, 1250, 1259, 1261, [1287,1288], [1293,1294], 1848, [1851,1853], 1907, [1979,1980], 2019, 2025, 2069, [2072,2074], 2094, 2096, [2098,2099], 2103, 2105, 2108, 2122, 2124, 2126, [2144,2145], [2148,2149], [2160,2161], [2173,2174], 2216, [2223,2226], 2234, 2236, 2247, [2250,2251], 2255, 2258, 2275, 2280, [2286,2287], [2290,2291], [2322,2323], 2325, [2328,2329], 2341, [2361,2362], [2368,2369], [2373,2374], [2381,2382], 2403, 2407, 2410, 2422, 2428, 2442, 2453, 2457, 2459, [2472,2473], 2475, [2477,2478], 2482, 2518, 2520, 2522, 2524, 2526, 2529, 2531, 2533, 2535, 2537, [2555,2556], 2560, 2562, 2564, 2594, 2596, [2598,2599], 2603, 2633, 2635, 2637, [2639,2641], [2697,2698], [2701,2703], 2708, [2711,2712], [2718,2720], [2747,2748], [2750,2751], 2759, [2791,2792], 2818, 2823, [2832,2834], 2838, [2842,2843], 2866, 2868, 2881, 2902, 2910, [2916,2919], [2932,2933], 2947, [2960,2961], [2971,2972], 2976, 2978, 2980, [3008,3009], 3013, 3015, 3017, 3022, 3029, 3031, 3045, [3050,3051], 3062, 3067, 3076, 3078, [3104,3105], [3110,3111], 3665, [3668,3670], 3724, [3796,3797], 3836, 3842, 3886, [3889,3891], 3911, 3913, [3915,3916], 3920, 3922, 3925, 3939, 3941, 3943, [3961,3962], [3965,3966], [3977,3978], [3990,3991], 4033, [4040,4043], 4051, 4053, 4064, [4067,4068], 4072, 4075, 4092, 4097, [4103,4104], [4107,4108], [4139,4140], 4142, [4145,4146], 4158, [4178,4179], [4185,4186], [4190,4191], [4198,4199], 4220, 4224, 4227, 4239, 4245, 4259, 4270, 4274, 4276, [4289,4290], 4292, [4294,4295], 4299, 4335, 4337, 4339, 4341, 4343, 4346, 4348, 4350, 4352, 4354, [4372,4373], 4377, 4379, 4381, 4411, 4413, [4415,4416], 4420, 4450, 4452, 4454, [4456,4458], [4514,4515], [4518,4520], 4525, [4528,4529], [4535,4537], [4564,4565], [4567,4568], 4576, [4608,4609], 4635, 4640, [4649,4651], 4655, [4659,4660], 4683, 4685, 4698, 4719, 4727, [4733,4736], [4749,4750], 4764, [4777,4778], [4788,4789], 4793, 4795, 4797, [4825,4826], 4830, 4832, 4834, 4839, 4846, 4848, 4862, [4867,4868], 4879, 4884, 4893, 4895, [4921,4922], [4927,4928]}

If manually look in the set of active cells  for the DoFs  which belong to one of the "problematic" points i find

231 : comp =0   at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 0
636 : comp =0   at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 1
2048 : comp =1  at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 0
2453 : comp =1  at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 1
3865 : comp =2  at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 0
4270 : comp =2   at the node located in 4.5000000000e-01 2.5000000000e-01 2.0500000000e-01 on level 1

However the bold ones are not contained in the indexset, an do have troubles at this points, when i would like to apply constraintmatrix.distribute(newton_update);

Has somebody a clue what is wrong, or how i can resolve the problem.

With best Regard,

Bernhard Endtmayer

[Bernhard Endtmayer]

The problem appears in the marked node.

Thanks in advance for your help!

[Wolfgang Bangerth]

Bernhard,

So to summarize, you are interpolating boundary conditions on boundaries with
boundary_id 0, 2, 80, and then you are outputting all DoFs on boundary 80.

> The set is given by

> [...]

>
> If manually look in the set of active cells  for the DoFs  which belong to one
> of the "problematic" points i find
>
>

> *231* : comp =0   at the node located in 4.5000000000e-01 2.5000000000e-01

> 2.0500000000e-01 on level 0
> 636 : comp =0   at the node located in 4.5000000000e-01 2.5000000000e-01
> 2.0500000000e-01 on level 1

> *2048* : comp =1  at the node located in 4.5000000000e-01 2.5000000000e-01

> 2.0500000000e-01 on level 0
> 2453 : comp =1  at the node located in 4.5000000000e-01 2.5000000000e-01
> 2.0500000000e-01 on level 1

> *3865* : comp =2  at the node located in 4.5000000000e-01 2.5000000000e-01

> 2.0500000000e-01 on level 0
> 4270 : comp =2   at the node located in 4.5000000000e-01 2.5000000000e-01
> 2.0500000000e-01 on level 1
>
> However the bold ones are not contained in the indexset, an do have troubles
> at this points, when i would like to apply
> constraintmatrix.distribute(newton_update);

Even though I might have an idea where you are going, your message does not
actually contain a *question* :-) Can you describe what you expect, what you
see, and how the two differ? In other words, what the actual *problem* is?

Best
W.


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

[Bernhard Endtmayer]

Dear Wolfgang Bangerth,

Thank you for your answer.

The problem is that at this nodes the boundaryconditions do not apply.

So the question would be:

How i obtain the correct boundary values for nodes on the boundary?

or

How can i find all nodes which belong to the actual boundary and set the correct boundary conditions.

The problem is that interpolate_boundary_values does not find all nodes on the boundary after adaptive refinement.

If i consider the newton update, the correct values are not imposed on the update, and the solution gets wrong boundary conditions on the bold nodes.

However this just happens on a few different meshes during the adaptive process.

>So to summarize, you are interpolating boundary conditions on boundaries with
>boundary_id 0, 2, 80, and then you are outputting all DoFs on boundary 80.

Yes, this problem just happens on boundary 80.

After i gave them as output, i recognized that not all DoFs which schould be on the boundary are actually at the boundary.

Best Wishes,

Bernhard

[Wolfgang Bangerth]

Bernhardt,

> How i obtain the correct boundary values for nodes on the boundary?
> or
> How can i find all nodes which belong to the actual boundary and set the
> correct boundary conditions.
>
> The problem is that interpolate_boundary_values does not find all nodes on the
> boundary after adaptive refinement.

It's not impossible that you're right, but it would also surprise me that 20
years after this function was written, it would still have bugs :-) Do you
happen to know whether the missing DoFs happen to be those that are also
hanging nodes?

> >So to summarize, you are interpolating boundary conditions on boundaries with
> >boundary_id 0, 2, 80, and then you are outputting all DoFs on boundary 80.
>
> Yes, this problem just happens on boundary 80.
> After i gave them as output, i recognized that not all DoFs which schould be
> on the boundary are actually at the boundary.

I think you'll have to try to narrow down what exactly happens. If you call
VectorTools::interpolate (..., 80, ...);
then my expectation is that the output argument contains constraints for every
DoF located on boundary 80. You can find out which DoFs are located on
boundary 80 using DoFTools::extract_boundary_dofs(), and then you could loop
over all of these extracted DoFs and ask the AffineConstraints
(ConstraintMatrix) class whether these DoFs are constrained using
constraints.is_constrained(i);
If you get a `false` answer for any of these DoFs, then there is indeed
something wrong. You could then double-check whether these DoFs are maybe
hanging nodes using DoFTools::extract_hanging_nodes().