Hi Huan, The distribute_local_to_global modifies both the right hand side and the matrix itself in how it preforms the algorithm. Hence it is not enough just to give it a dummy matrix (in which case the modifications will be incorrect) or even to give it a copy of the actual matrix, but not use that matrix later.
If you look at step-22 in the deal-6.1.0 tutorial, you will notice that the boundary conditions are implemented in a completely different way, using VectorTools::interpolate_boundary_values and MatrixTools::apply_boundary_values. This method is still available I believe however the constraint matrix class seems to be preferred in general. Regards, Michael On Tue, May 18, 2010 at 8:12 AM, Huan Sun <[email protected]> wrote: > Hi Michael, > > Thank you for the suggestions. I will look at the filtered matrix > class. > I understand that distribute_local_to_global() usually requires matrix > and righthand side at the same time. That's why I need an auxiliary > matrix matrix_for_bc. > I was basically following tutorial-34 and my Poisson test is much > simpler than the hydrodynamic system solved there. So I guess I am using > ConstraintMatrix in a legitimate way ... somehow something goes wrong. > > Also I have one more question on the inhomogeneous boundary condition: > in this case is VectorTools::project() forbidden? > > Thanks, > > Huan > > > > > > > > > On Tue, 2010-05-18 at 04:12 +0200, Michael Rapson wrote: >> Hi Huan, >> >> I also sometimes separate the matrix and right hand side assembly. >> However, care must be taken because the method >> distribute_local_to_global uses both the right hand side matrix and >> the system matrix. Basically, the algorithm scales the Dirichlet >> boundary condition relative to the size of the terms in the matrix. In >> generally it will not work to apply the method to matrix and right >> hand side separately (hence the method which takes both) or even to >> apply the method again on a matrix that has had some of the terms >> assembled already. >> >> One option in the library is to use the FilteredMatrix class. (Which >> should work as long as you don't need to add terms to the matrix.) >> Another is to keep a copy of the matrix without the boundary >> conditions applied which you keep initializing from, completing with >> any extra terms and then use to apply the boundary conditions to. >> >> Regards, >> Michael >> >> On Mon, May 17, 2010 at 8:18 PM, Huan Sun <[email protected]> wrote: >> > Hi Markus, >> > >> > Yes, it works to use the original local stiffness matrix. But I want to >> > later separate the assemblage of the matrix and the right-hand side, >> > similar to what tutorial 32 does. So I was testing this easier case. >> > Somehow it doesn't work with this auxiliary matrix. I am just wondering >> > if there is anything wrong or missing in the code. >> > >> > Thanks, >> > >> > Huan >> > >> > >> > On Mon, 2010-05-17 at 08:29 +0200, Markus Bürg wrote: >> >> Hello Huan Sun, >> >> >> >> why do you create another matrix for the boundary conditions? You have >> >> already created the problem's matrix for this cell right above. Thus >> >> you can hand over just this matrix: >> >> constraints.distribute_local_to_global (local_rhs, local_dof_indices, >> >> rhs, local_matrix); >> >> >> >> Best Regards, >> >> Markus >> >> >> >> >> >> >> >> Am 17.05.10 08:12, schrieb Huan Sun: >> >> > Hi all, >> >> > >> >> > I was following tutorial 22 and 34 to implement inhomogeneous Dirichlet >> >> > boundary conditions using the ConsraintMatrix class but couldn't get it >> >> > right. I think the problem mainly lies in the assemblage part (i was >> >> > solving the basic Poisson's eq.) >> >> > >> >> > if(constraints.is_inhomogeneously_constrained >> >> > (local_dof_indices[i])){ >> >> > for(unsigned int j=0; j <dofs_per_cell; ++j){ >> >> > matrix_for_bc(j,i) += >> >> > grad_basis_phi[i] * grad_basis_phi[j] * >> >> > fe_values.JxW(q); >> >> > }//j-loop >> >> > }//if >> >> > ... >> >> > constraints.distribute_local_to_global(local_rhs, >> >> > local_dof_indices, >> >> > rhs, >> >> > matrix_for_bc); >> >> > >> >> > However, when I comment out the if statement, the program seems to work >> >> > correctly. The code is attached in this email. >> >> > Your help will be greatly appreciated. >> >> > >> >> > Thanks! >> >> > >> >> > Huan >> >> > >> >> > >> >> > _______________________________________________ >> >> > dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii >> >> > >> >> _______________________________________________ >> >> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii >> > >> > >> > _______________________________________________ >> > dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii >> > > > > _______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
