Dear Martin,
what Luca suggested seems to be a simply solution. Of course you can
copy data and do everything by yourself, too, but that is a bit messy,
see below...
> Essentially, if I could extract the matrix entries that correspond to
> the columns for the non-zero boundary constraints, I could compute the
> right hand side contributions in a more straightforward way.
You can store these columns with one double vector "column_values" and
two unsigned int vectors, "row_indices" and "column_start" (basically
you create a compressed column storage with the columns you need). As
you are solving a Stokes system, the matrix should be symmetric, so you
can extract rows instead of columns:
1. Assemble the matrix, but do not yet call apply_boundary_values
2. Use interpolate_boundary_values with std::map argument. This gives
you the numbers of columns that you need to extract
3. Set the length of row_start to the size of boundary_values + 1. It
will store at which position in the vector the individual columns start,
according to the compressed row storage formant, but now applied to the
columns
4. Go through the map. Record the following:
a) the length of the row (matrix_rowstarts[current_dof
+1]-matrix_rowstarts[current_dof])
b) copy all indices from matrix_colnum with indices between
matrix_rowstarts[current_dof] and matrix_rowstarts[current_dof+1] to
row_indices
c) copy the data into column_values by accessing it
sparse_matrix.global_entry (index), index runs starts from
matrix_rowstarts[current_dof]
The data is accessed as follows:
const std::size_t * matrix_rowstarts =
sparse_matrix.get_sparsity_pattern().get_rowstart_indices();
const unsigned int * matrix_colnums =
sparse_matrix.get_sparsity_pattern().get_colnum_numbers();
5. Call apply_boundary_values to eliminate the entries from the matrix.
Now this gives you the data of columns in the matrix that you need.
In each time step do:
- interpolate_boundary_values fills a std::map
- Go through the entries of the map, and apply the matrix columns (Gauss
elimination, but unmodified diagonal (needs to be nonzero!)):
for (unsigned int i=column_start[map_index]+1;
index<column_start[map_index+1]; ++i)
rhs_vector(row_indices[i]) -= map_iterator->second * column_values[i];
rhs_vector(map_iterator->first) = map_iterator->second;
map_index increments with map_iterator, and you also use that
column_values[column_start[map_index]] stores the diagonal entry, which
we skip here (since it is set later on).
Best,
Martin
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
