The apply_boundary_values does several things, which you have to keep track of,
if you want to reuse the factorization:
1. eliminate the rows of the matrix corresponding to Dirichlet Data, putting a
value in the diagonal entry which is consistent with what has been taken away
(to preserve the condition number of your matrix), (we'll call these entries
a_ii)
2. put the correct values of the boundary conditions on your rhs ( set rhs[i] =
boundary_values[i] for each i which exists in the boundary_values map )
3. scale the newly inserted value with a_ii, so that the final solution would
be the same
4. (optionally) perform a Gauss elimination step, with the changed matrix and
the changed rhs to eliminate the column and possibly preserve the matrix
symmetry
Some of these operations are destructive. The most destructive operation is n.
4, since it modifies the rhs in a non trivial way, and you won't have a way to
do the same thing with another rhs, unless you start from the original matrix
(which you no longer have). Hence the reason the set the flag *not* to
eliminate the column.
To treat the others, as observed by Marcus, you have different choices. Mine is
the following:
1. I construct a boundary_values map interpolating the ConstantFunction<dim>
one(1) on my dirichlet portion of the boundary
2. I create an empty "diagonal_entries" vector
3. I call apply_boundary_values with "diagonal_entries" as rhs
in the code you sent me, your rhs is initialized with "solution". I'd do
diagonal_entries.reinit(dh.n_dofs()).
This makes an empty vector. (your rhs is my diagonal_entries).
After the call, in diagonal_entries you'll have the a_ii entries of the matrix
for each dirichlet dof (because you used a "1" boundary value). Now you can
factorize the matrix.
Suppose you have the b_values you want to apply, in b_values, and a map
iterator it:
for(it = b_values.begin(); it != b_values.end(); ++it) {
system_rhs(it->first) = it->second* diagonal_entries(it->first);
}
and you should be good to go, i.e.,
Ainv.vmult(solution, system_rhs)
should give you what you want.
Luca
--
Luca Heltai <[email protected]>
http://people.sissa.it/~heltai/
Scuola Internazionale Superiore di Studi Avanzati
Phone: +39 040 3787 449, Office: 732
--
There are no answers, only cross references
On Nov 25, 2011, at 6:06 AM, Marc Secanell Gallart wrote:
> Hello Markus and Luca,
>
> Thank you very much for your suggestions and explanation.
>
> It seems to me that both approaches are similar, however, Luca's approach
> will allow me to reuse the factorization with different RHSs which, please
> correct me if I am mistaken, I might not be able to do using Markus' approach.
>
> Luca, you mention "Apply boundary values to the matrix, with a fake
> boundary_values map, in which all entries which you want to later apply are
> set to one, a zero rhs" So, do you mean something like this
>
> BlockVector<double> rhs;
> rhs.reinit(solution);
> std::map<unsigned int, double> boundary_values;
> VectorTools::interpolate_boundary_values (dof,
>
> grid.get_boundary_id("c_CL/Membrane"),
>
> ConstantFunction<dim> PEM(1),
>
> boundary_values,
>
> comp_mask);
> MatrixTools::apply_boundary_values(boundary_values,
> matrix,
> solution,
> rhs,
> false);
>
> Then you mention "Factorize the matrix, then when you want to apply it to
> different bondary conditions, scale the values of the rhs corresponding to
> the dirichlet conditions with the previously stored rhs values" What do you
> mean by "scale the values of the rhs (...)"?
>
> Thank you,
>
> Marc
>
> _______________________________________________
> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
