Hello all, I am solving what is essentially a Navier-Stokes system with deal.II. However, when I attempt to implement the optimizations suggested in the "Handling vector-valued problems" module <https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html>. However, I am not getting the same matrix when I assemble in the optimized fashion.
My understanding thus far: Given that my finite element space is just a stack of primitive finite elements (Q_k for velocity P_{k-1} discontinuous for pressure), I can recover the nonzero component of the current shape function with something like: component_i = fe.system_to_component_index(i).first; I also understand that mass and Laplace matrices are zero except when the components are equal. Thus I cannot figure out why this code does not work properly. All of the other usual infrastructure is there such as a loop over the cells and quadrature points, multiplying by JxW[q], and so on. for (unsigned int i = 0; i < dofs_per_cell; ++i){ // We don't need to even do anything if the component doesn't correspond to a velocity dof. const unsigned int component_i = fe.system_to_component_index(i).first; for (unsigned int j = 0; j < dofs_per_cell; ++j){ // We don't need to even do anything if the component doesn't correspond to a velocity dof. const unsigned int component_j = fe.system_to_component_index(j).first; if (component_i < 3 && component_j < 3){ // Just do some checking... if (component_i!=component_j){ if ((1.0 / dt) * phi_u[j] * phi_u[i] + (1.0 / dt) * delta2 * scalar_product(grad_phi_u[j], grad_phi_u[i]) + 0.5 * viscosity * scalar_product(grad_phi_u[j], grad_phi_u[i]) != 0.0){ std::cout << "components " << (int)component_i << " and " << (int)component_j << '\n'; std::cout << "(" << phi_u[i][0] << "," << phi_u[i][1] << "," << phi_u[i][2] << ")\n"; std::cout << "(" << phi_u[j][0] << "," << phi_u[j][1] << "," << phi_u[j][2] << ")\n"; std::cout << "dot product is " << phi_u[i] * phi_u[j] << std::endl; } } local_matrix(i, j) += ( (component_i!=component_j) ? 0.0 : ( (1.0/dt) * phi_u[j] * phi_u[i] + (1.0/dt) * delta2 * scalar_product(grad_phi_u[j], grad_phi_u[i]) + 0.5 * viscosity * scalar_product(grad_phi_u[j], grad_phi_u[i]) ) .... and so on with other terms which I have not changed in my attempts to optimize. My confusion is amplified by the block which I have labeled "Just do some checking..." There is no output coming from those print statements, which seems to support my understanding of what is happening. Any help you can offer would be appreciated. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/c6cca5ac-193b-42d5-931f-baa6dd542e1en%40googlegroups.com.