Hi Markus and Michael,
Thank you for the replies.
I tested using constraints.is_constrained but the results doesn't
improve. It only works if I comment out the whole IF statement, i.e.
assemble the dummy matrix on every DOF.
The main reason I want to separate the matrix and rhs assembly is I want
to do sth. similar to tutorial 34, computing some 3 dimenional
hydrodynamic system in parallel. The the stiffness matrix does not need
to be assembled at every time step. I guess in this case
MatrixTools::apply_boundary_values() does not work in this case. (please
correct me if I'm wrong.)
Also when I used VectorTools:project() for inhomogeneous boundary
condition, I got some exception message as follows, not sure what's
going on.
--------------------------------------------------------
An error occurred in line <162> of file
</gpfs/home/hus126/work/download/deal.II/lac/include/lac/constraint_matrix.templates.h>
in function
void dealii::ConstraintMatrix::condense(VectorType&) const [with
VectorType = dealii::Vector<double>]
The violated condition was:
constraint_line->inhomogeneity == 0.
The name and call sequence of the exception was:
ExcMessage ("Inhomogeneous constraint cannot be condensed " "without
any matrix specified.")
Additional Information:
Inhomogeneous constraint cannot be condensed without any matrix
specified.
Stacktrace:
-----------
#0 /gpfs/home/hus126/work/download/deal.II/lib/liblac.g.so.6.2.1: void
dealii::ConstraintMatrix::condense<dealii::Vector<double>
>(dealii::Vector<double>&) const
#1 /gpfs/home/hus126/work/download/deal.II/lib/libdeal_II_2d.g.so.6.2.1: void
dealii::VectorTools::project<2, dealii::TrilinosWrappers::Vector,
2>(dealii::Mapping<2, 2> const&, dealii::DoFHandler<2, 2> const&,
dealii::ConstraintMatrix const&, dealii::Quadrature<2> const&,
dealii::Function<2> const&, dealii::TrilinosWrappers::Vector&, bool,
dealii::Quadrature<(2)-(1)> const&, bool)
#2 /gpfs/home/hus126/work/download/deal.II/lib/libdeal_II_2d.g.so.6.2.1: void
dealii::VectorTools::project<2, dealii::TrilinosWrappers::Vector,
2>(dealii::DoFHandler<2, 2> const&, dealii::ConstraintMatrix const&,
dealii::Quadrature<2> const&, dealii::Function<2> const&,
dealii::TrilinosWrappers::Vector&, bool, dealii::Quadrature<(2)-(1)> const&,
bool)
#3 ./angle_newton: MultiPhaseProblem<2>::run()
#4 ./angle_newton: main
this can be circumvented by solving a mass matrix inversion. So this
problem does not bother me so far. But I would appreciate any suggestion
for the above error message.
Thanks,
Huan
On Tue, 2010-05-18 at 09:10 +0200, Markus Bürg wrote:
> Hello Huan,
>
> I see. I did similar tests some time ago. The only difference I did was
> using constraints.is_constrained instead of
> constraints.is_inhomogeneously_constraint. Can you try, if this works?
>
> Anyway for me this method did not pay out. It took much longer than
> assemling all in one function and saving the matrix, because most of the
> time is - of course - spent on assembling the matrix. I think this
> implementation is only favourable, if one has homogeneous Dirichlet
> boundary conditions.
>
> VectorTools::project () should also work for inhomogeneous boundary
> conditions.
>
> Best Regards,
> Markus
>
>
>
> Am 18.05.10 08:51, schrieb Michael Rapson:
> > 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
> >
> >
> _______________________________________________
> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
