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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1309f302-7a9a-4033-90d7-a7b91c4a48d6n%40googlegroups.com.
