Reinit a 2D FEBase on a side of a 3D TET10

[rebpoli]

Hi,

I am trying to map the aperture of a fracture on the sides of a TET10.

That means that only the DOFs on that side really matter, all other have zero-weight on the mapping and some of the diagonal elements of the "M" matrix (see below) are zeroed out, breaking the lu_solve.

I added a dusty assignment to one whenever that happens, and that works.
However, I would be much nicer if I could iterate and map only through the side DOFs.

Is there a way to do that?

Thanks,
Renato

for (uint qp=0; qpn_points(); qp++) {
for(uint i=0; i<n_dofs; i++)
{
F(i) += jxw[qp] * vals_qp[qp] * phi[i][qp];
for(uint j=0; j<n_dofs; j++) M(i,j) += jxw[qp] * phi[i][qp] * phi[j][qp];
}
for(uint i=0; i<n_dofs; i++) if ( std::abs(M(i,i)) < 1e-10 ) M(i,i) = 1;
}
DenseVector projected_data;
M.lu_solve(F, projected_data);

[roystgnr]

Take a look at FEInterface::dofs_on_side()? That's what we use to do these sorts of problems in generic_projector.h. In that context we typically do still have to loop over every dof (we have forcing functions coming from vertex and edge dof values by the time we're projecting side dofs) but we at least only have to solve a linear system with side dofs alone.