Hi, thanks all for the suggestions.
As I am doing an all-(time-steps)-at-once approach the Filtered matrix is probably not useful as I would have to store N_T times Filtered matrix with different constraints, which for that format would mean to store N_T saddle point systems in advance. It looks as I will follow Martin's and Luca's suggestion. Thank you all so much for your help. Best, Martin On Wed, Jan 12, 2011 at 1:23 PM, Martin Kronbichler < [email protected]> wrote: > 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 > > -- *Martin Stoll* *Postdoctoral Research Fellow* Computational Methods in Systems and Control Theory Max Planck Institute for Dynamics of Complex Technical Systems Sandtorstr. 1 D-39106 Magdeburg Germany Email: [email protected] URL : http://www.mpi-magdeburg.mpg.de/people/stollm Tel :+49 391 6110 384
_______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
