Following up on this email.
John Peterson, Roy Stogner, any ideas?
On 02/04/2016 10:35 AM, Rob Blake wrote:
tl;dr: I think the 4 argument version of
DofMap::constrain_element_matrix is bugged, but I can't fix it without
breaking the moose tests. Help?
----
Is the 4 argument version of constrain_element_matrix bugged?
I tried to introduce incompressibility using a lagrange multiplier
with hydrostatic pressure. My code was very much like
http://libmesh.github.io/examples/systems_of_equations_ex3.html.
When I change
dof_map.constrain_element_matrix(Ke, dof_indices);
navier_stokes_system.matrix->add_matrix (Ke, dof_indices);
into something like this (paraphrasing):
auto _dof_indices_d = dof_indices_d;
dof_map.constrain_element_matrix(Kdd, _dof_indices_d);
navier_stokes_system.matrix->add_matrix (Kdd, _dof_indices_d);
auto _dof_indices_p = dof_indices_p;
dof_map.constrain_element_matrix(Kpp, _dof_indices_p);
navier_stokes_system.matrix->add_matrix (Kpp, _dof_indices_p);
_dof_indices_p = dof_indices_p;
_dof_indices_d = dof_indices_d;
dof_map.constrain_element_matrix(Kdp, _dof_indices_d, _dof_indices_p);
navier_stokes_system.matrix->add_matrix (Kpp, _dof_indices_d, _dof_indices_p);
_dof_indices_p = dof_indices_p;
_dof_indices_d = dof_indices_d;
dof_map.constrain_element_matrix(Kpd, _dof_indices_p, _dof_indices_d);
navier_stokes_system.matrix->add_matrix (Kpd, _dof_indices_p, _dof_indices_d);
the code no longer works. Dumping the element matrices, I see this
for the pressure row in Ke for elements with a dirichlet boundary
condition:
0.000000e+00
0.000000e+00
0.000000e+00
6.250000e+04
-6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
-6.250000e+04
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
6.250000e+04
6.250000e+04
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
versus this form for the asymmetric case:
-6.250000e+04
-6.250000e+04
-6.250000e+04
6.250000e+04
-6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
-6.250000e+04
-6.250000e+04
6.250000e+04
-6.250000e+04
-6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
6.250000e+04
6.250000e+04
-6.250000e+04
6.250000e+04
6.250000e+04
This suggests that the code is not applying the dirichlet BC correctly
for the displacement/pressure matrix when I use the Kdp version of
constrain element matrix. Reading the code, my suspicion is confirmed:
// It is possible that the matrix is not constrained at all.
if ((R.m() == matrix.m()) &&
(R.n() == row_dofs.size()) &&
(C.m() == matrix.n()) &&
(C.n() == col_dofs.size())) // If the matrix is constrained
{
// K_constrained = R^T K C
matrix.left_multiply_transpose (R);
matrix.right_multiply (C);
Indeed, in this case R != NULL but C == NULL, so the if statement is
false. We incorrectly do not multiply by R.
I tried to fix this in various ways, but the tests always fail:
https://github.com/libMesh/libmesh/pull/821
Any insights?
------------------------------------------------------------------------------
Site24x7 APM Insight: Get Deep Visibility into Application Performance
APM + Mobile APM + RUM: Monitor 3 App instances at just $35/Month
Monitor end-to-end web transactions and take corrective actions now
Troubleshoot faster and improve end-user experience. Signup Now!
http://pubads.g.doubleclick.net/gampad/clk?id=272487151&iu=/4140
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
------------------------------------------------------------------------------
Site24x7 APM Insight: Get Deep Visibility into Application Performance
APM + Mobile APM + RUM: Monitor 3 App instances at just $35/Month
Monitor end-to-end web transactions and take corrective actions now
Troubleshoot faster and improve end-user experience. Signup Now!
http://pubads.g.doubleclick.net/gampad/clk?id=272487151&iu=/4140
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
