Hello everyone!
I'm trying to rewrite Step-35 (Navier-Stokes tutorial) without splitting
the implementation for component (In Step 35 one thread resolve the x
component and the other one the y component).
Proceeding in this way I can use "AffineConstraints" class to manipulate
constraints in the system matrices.
All went well until I tried to implement a mesh refinement.
In particular, when I assemble the gradient operator after the refinement
of the mesh (increasing the dof of the problem), this error occurs when
"distribute_local_to_global" function is called:
*The violated condition was: (this_cols[counter] == col_indices[i]) ||
(values[i] == number2())Additional information: You are trying to
access the matrix entry with index <35518,4519>, but this entry does not
exist in the sparsity pattern of this matrix.The most common cause for this
problem is that you used a method to build the sparsity pattern that did
not (completely) take into account all of the entries you will later try to
write into. An example would be building a sparsity pattern that does not
include the entries you will write into due to constraints on degrees of
freedom such as hanging nodes or periodic boundary conditions. In such
cases, building the sparsity pattern will succeed, but you will get errors
such as the current one at one point or other when trying to write into the
entries of the matrix.*
I found in the documentation the "add_entries" command, in such a way to
add the matrix entry in the sparsity pattern, but I don't know if it is the
right way to proceed. Any suggestions???
Thanks a lot
The piece of the code, where the gradient operator is assembled, is
reported below:
/* Gradient operator*/
DynamicSparsityPattern dsp3(dof_handler_velocity.n_dofs(),
dof_handler_pressure.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler_velocity,
dof_handler_pressure,
dsp3);
sparsity_pattern_pres_vel.copy_from(dsp3);
pres_Diff.reinit(sparsity_pattern_pres_vel);
unsigned int vel_dpc(fe_velocity.dofs_per_cell);
unsigned int pres_dpc(fe_pressure.dofs_per_cell);
FullMatrix<double> local_grad(fe_velocity.dofs_per_cell,
fe_pressure.dofs_per_cell);
std::vector<types::global_dof_index>
vel_local_dof_indices(fe_velocity.dofs_per_cell);
std::vector<types::global_dof_index>
pres_local_dof_indices(fe_pressure.dofs_per_cell);
unsigned int nqp(quadrature_velocity.size());
FEValues<dim> fe_val_vel(fe_velocity, quadrature_velocity,
update_values | update_gradients |
update_JxW_values);
FEValues<dim> fe_val_pres(fe_pressure, quadrature_velocity,
update_values);
auto vel_cell = dof_handler_velocity.begin_active();
const auto endc = dof_handler_velocity.end();
auto pres_cell = dof_handler_pressure.begin_active();
for (; vel_cell != endc; ++vel_cell, ++pres_cell) {
fe_val_vel.reinit(vel_cell);
fe_val_pres.reinit(pres_cell);
local_grad = 0;
for (unsigned int q_point = 0; q_point < nqp; ++q_point) {
for (unsigned int i = 0; i < vel_dpc; ++i) {
const unsigned int component_i =
fe_velocity.system_to_component_index(i).first;
for (unsigned int j = 0; j < pres_dpc; ++j) {
local_grad(i, j) +=
fe_val_pres.shape_value(j, q_point)*
fe_val_vel.shape_grad(i, q_point)[component_i] *
(-fe_val_vel.JxW(q_point));
}
}
}
vel_cell->get_dof_indices(vel_local_dof_indices);
pres_cell->get_dof_indices(pres_local_dof_indices);
velocity_hanging_node_constraints.distribute_local_to_global(local_grad,
vel_local_dof_indices,
pressure_hanging_node_constraints,
pres_local_dof_indices, pres_Diff);
}
--
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/01ca2aa2-120d-4912-bbf8-7c33b7f32176%40googlegroups.com.
