The easiest way to do what you want is by using matrix.zero_rows(). The way it
works is that you fill up a vector of dof indices that correspond to rows that
you want to zero and put a 1 on the diagonal. Then you just call
matrix.zero_rows(dofs_to_zero, 1). The second argument to the function is the
number you want to place on the diagonal (zero by default).
So.... fill up your matrix like normal, keeping track of rows you want to wipe
out. Close the matrix. Call zero_rows(). Close the matrix again.
A quick note on this... if you are using Petsc this will modify your nonzero
pattern (ie remove all off-diagonal nonzero entries for these rows). That's
not normally what you want because the next time you come around to fill the
matrix and you do some insertions on that row it will be incredibly slow
(because Petsc will have to reallocate those entries). You can get beyond this
by using something like:
MatSetOption(static_cast<PetscMatrix<Number> &>(my_matrix).mat(),
MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
That is the correct line for Petsc 3.1. If you are using a different version
there will most likely be a slightly different way to tell it to keep the
nonzero pattern. I recommend setting this each time before you start filling
the matrix up. Petsc seems to like to forget about this option.....
Derek
On May 10, 2011, at 5:10 PM, Omar Al-Abbasi wrote:
> Thank you all for all the suggestions you made.
> In fact I tried doing two seperate sweeps as Boyce suggested, but it seems
> that I am doing a mistake which leads to no change in mass matrix.
>
> This how I am doing that after the first sweep for assembling the matrices:
> .
> .
> .
> for ( ; el != end_el; ++el)
> {
> const Elem* elem = *el;
> dof_map.dof_indices (elem, dof_indices);
> fe->reinit (elem);
> Ke.resize (dof_indices.size(), dof_indices.size());
> Me.resize (dof_indices.size(), dof_indices.size());
>
> for (unsigned int side=0; side<elem->n_sides(); side++){
> if (elem->neighbor(side) == NULL)
> {
> const std::vector<std::vector<Real> >& phi_face =
> fe_face->get_phi();
> fe_face->reinit(elem, side);
> for (unsigned int i=0; i<phi_face.size(); i++){
> matrix_B.set(dof_indices[i],dof_indices[i],1.);
> } // end if (elem->neighbor(side) == NULL)
> } // end boundary condition section
> }
> }
>
> Also, what is the best way to loop over the nodes which live on the outer
> surface of the domain?
> Is it just the way I am doing it for (unsigned int side=0;
> side<elem->n_sides(); side++){
> if (elem->neighbor(side) == NULL)
> and then for (unsigned int i=0; i<phi_face.size(); i++){ ....
>
> Thanks again,
> Omar
> ------------------------------------------------------------------------------
> Achieve unprecedented app performance and reliability
> What every C/C++ and Fortran developer should know.
> Learn how Intel has extended the reach of its next-generation tools
> to help boost performance applications - inlcuding clusters.
> http://p.sf.net/sfu/intel-dev2devmay
> _______________________________________________
> Libmesh-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
------------------------------------------------------------------------------
Achieve unprecedented app performance and reliability
What every C/C++ and Fortran developer should know.
Learn how Intel has extended the reach of its next-generation tools
to help boost performance applications - inlcuding clusters.
http://p.sf.net/sfu/intel-dev2devmay
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
