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.

Reply via email to