[Libmesh-users] Non-symmetric constrained element matrices in case of recursive constraints

[Jonas B. <jon...@ep...>]

Dear all,

I have a question regarding constrained dofs that are themselves constrained by other dofs within the same element. Apparently, this leads to a non-symmetric constrained element matrix which I illustrate here by a small example:

I consider an element with dof indices = (0,1,2,3). I impose the following contraint rows

{
  constrained_dof_index = 0;
  constraint_row[1] = 1.0;
  constraint_row[2] = -1.0;
}
{
  constrained_dof_index = 1;
  constraint_row[0] = 1.0;
  constraint_row[3] = -1.0;
}

using

system.get_dof_map().add_constraint_row( constrained_dof_index, constraint_row, 0., true);

Since dof 1 is also constrained, these constraints get preprocessed so that

dof_map.print_dof_constraints();

yields

Constraints for DoF 0:      (0,1)     (2,-1)     (3,-1)    rhs: 0
Constraints for DoF 1:      (2,-1)     (3,-2)    rhs: 0

This means that the following constraint matrix is generated:

C =
     1     0    -1    -1
     0     0    -1    -2
     0     0     1     0
     0     0     0     1

Assume now that the (symmetric) element matrix K is given by

K =
     1     2     3     4
     2     3     4     5
     3     4     5     6
     4     5     6     7

Then

C^T*K*C =
     1     0     0    -1
     0     0     0     0
     0     0    -1    -1
    -1     0    -1     0

This means that

dof_map.constrain_element_matrix(K, dof_indices, /*asymmetric_constraint_rows*/false);

generates the (non-symmetric) constrained matrix

K_constrained =
       1        0        0        0
       0        1        0        0
       0        0       -1       -1
      -1        0       -1        0

because the first two rows are replaced by unit rows.

My question would be if this cannot be avoided so that I also have to use a non-symmetric solver in this case.

Best,
Jonas