Re: [deal.II] Issue with convergence of iterative linear solver for system matrix in modified step-57 with no-normal flux constraints

[Bruno Blais]

I have found my problem. A boundary condition in my mesh was ill-defined 
and this lead to my error.
Sorry for the confusion.
Thanks


On Monday, 19 February 2018 13:41:35 UTC-5, Bruno Blais wrote:
>
> Hello,
> Sorry I feel I have not explained myself correctly. Here is a drawing of 
> the case:
>
>
> 
>
> With Ux may be a profile or a constant.
> Initially I had set-up that case with Ux as a Poiseuille flow profile and 
> the Top and Bottom as no slip boundary conditions. In that case I got the 
> solution I was expecting. Convergence was good and the solution was what I 
> had in mind.
> Then what I wish to do is set Ux as a constant and apply slip at the top 
> and bottom boundary conditions. When I do this is when I am unable to reach 
> convergence in my iterative solution, yet I do not understand why because 
> the problem should be well posed.
> I don't think that leaving my right boundary as a free traction boundary 
> condition is problematic since it is by all mean an outlet boundary. This 
> worked well for instance in a backward facing step case.
> Thank you for your time, sorry for the confusion.
> Best regards
>
>
>
> On Monday, 19 February 2018 10:50:48 UTC-5, Timo Heister wrote:
>>
>> > Does the order in which I apply the nonzero and zero constraints 
>> matter? 
>>
>> These are two independent objects, so no. 
>>
>> > Currently I apply the inlet and then the no-slip in the 
>> nonzero_constraints, 
>> > thus the bottom and top wall appear after the inlet. Afterward the 
>> cylinder 
>> > is put in the zero constraints. 
>>
>> I don't understand. You need boundary conditions in 
>> nonzero_constraints and zero_constraints for all boundaries. The 
>> reason for these two sets is that one is used for the initial solve, 
>> while the second one is used for the Newton updates. 
>>
>>
>> -- 
>> Timo Heister 
>> http://www.math.clemson.edu/~heister/ 
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: step-22 partial boundary conditions

[Jane Lee]

Hi Wolfgang, thanks so much for getting back to me

1. I was wondering re dimensions because I could find a component mask 
function or something similar to fe_values[velocities], eg, when using 
fe_face_values which you need to apply the neumann conditions. 

2. This is my fault. I meant working in 2d with a solution that only acts 
vertically, hence choosing a solution that is only dependent on y and that 
is 0 in the xdirection. 

3. indeed i was showing the full system right hand side component

4. for the top and bottom boundaries, you are right, i do mean the stress 
as mentioned in previous messages. Either way, I am now implementing things 
in the neumann BC way as in step-22. with the solution as before, i add the 
following to my top boundary:

for (unsigned int face_number=0; 
face_numberface(face_number)->at_boundary()

&&

(cell->face(face_number)->boundary_id() == 1))

{

fe_face_values.reinit (cell, face_number);



topstress.vector_value_list(fe_face_values.get_quadrature_points(),

topstress_values);


for (unsigned int q_point=0; q_point
>
> Jane, 
>
> > Firstly, I decided not to use the normal vector for now. Since the 
> > normal vector is 2D, i wasn't sure how to implement the rest so that it 
> > is a 'double' since my g (neumann condition vector) has 3 components 
> > when the normal vector will only have 2? 
>
> I'm not sure I understand this -- in 2 space dimensions, the normal 
> vector has 2 components; in 3 space dimensions, it has 3 components. 
>
>
> > I went back to the beginning - step-22 tutorial and started from 
> > absolute basics. Testing with 1D, i have started with exact solution 
> > u=(0,-y^2)^T and p = y^3, ie purely vertical. 
>
> I also don't understand this, I'm afraid: in d space dimensions, the 
> velocity vector has d components. So if you're in 1d, how come your 
> velocity vector has 2 components? 
>
> Similarly... 
>
> > Doing a manufactured 
> > solution type, i get the right hand side vector f = (0, 4+3z^2, 2z)^T. 
>
> the right hand side has d components. I suspect you are here showing the 
> right hand side of the whole system, which should then have d+1 
> components (which in 1d should be 2, not 3). 
>
> Or do you want to suggest in all of this that you really have a 
> two-dimensional domain but that the flow is one-dimensional? 
>
>
> > with my solution, on a rectangular domain, i have no tangential stresses 
> > anywhere, no flux on sides (unit normal being (/pm 1,0)), and plus and 
> > minus (0, z^3 + 4z) at the top. 
>
> Your domain is [0,1]^2, so I understand how there is no flux through the 
> bottom because there u=[0,0]. But how do you arrive at the flux you show 
> here? The velocity field is quadratic in z (y?), and so should the flux, 
> but you show something that is cubic with the distance. 
>
> If instead of 'flux' you really mean 'stress', then it should actually 
> only be linear in z since it involves the *derivative* of the velocity. 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received 

Re: [deal.II] project_boundary_values for Stokes equations

[Timo Heister]

Stephen,

I looked at the source code and it looks to me like the documentation
is correct: component_map is just to remap components and can not be
used to do what you want. It is just passed to
create_boundary_mass_matrix().

On Wed, Feb 21, 2018 at 7:03 PM, Stephen Shannon
 wrote:
> I am solving the 2D Stokes equations with Dirichlet BCs, so I have 3
> unknowns (u1, u2, p) with Dirichlet BCs on u1 and u2.  Rather than
> interpolate, I'd like to project the Dirichlet boundary values.
>
> Looking at the documentation for VectorTools::project_boundary_values(), the
> description of the "std::vector< unsigned int > component_map" input
> describes exactly this scenario, seemingly indicating that using a
> component_map will allow me to project values to 2 of the 3 components of my
> finite element.  However, reading the remainder of the description, and
> attempting to use the component_map, this seems to be not the case.  What
> I'd need is a component_map of {0,1}that is, the first boundary_function
> should apply to u1, and the second boundary_function should apply to u2.
> However, the map goes in the opposite directionthat is, it is required
> to have 3 components {0,1,?}: u1 gets the first boundary_function applied to
> it, u2 gets the second boundary_function applied to it, and p gets the ?
> boundary function applied to it.
>
> Is this part of the documentation incorrect when it mentions the Stokes
> problem as a use for the component_map?  In any case, what is a recommended
> method to projecting the boundary values for the Stokes problem?
>
> Thanks,
> Stephen
>
> --
> The deal.II project is located at 
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_=DwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=8teVs2Dl-uJLObDYEkzG-7e8YqH4XvNB-oZpCqBfVmg=fLQrhRDOCkvYUHpuoybOnpBqGVnhHifsf10HLYawDXo=
>  
> For mailing list/forum options, see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_forum_dealii-3Fhl-3Den=DwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=8teVs2Dl-uJLObDYEkzG-7e8YqH4XvNB-oZpCqBfVmg=3QTXEsJeXmRBshYVvxI8WlN7LZCI_N7PdahRDVXrSQc=
>  
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> For more options, visit 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_optout=DwIBaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=8teVs2Dl-uJLObDYEkzG-7e8YqH4XvNB-oZpCqBfVmg=EKokJBGcOix21OVFY2QpQDK3YBNsy-x11HCOtLvKLis=
>  .



-- 
Timo Heister
http://www.math.clemson.edu/~heister/

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: increment Dirichlet boundary condition

[Denis Davydov]

Hi IIiya,

On Thursday, February 22, 2018 at 11:55:33 AM UTC+1, Bryukhanov Ilya wrote:
>
> Denis, thanks a lot for the answer!
>
> I think that on the first step I can use usual 
> "interpolate_boundary_values + apply_boundary_values" functions,
>

the outline I gave above is agnostic to "first step"/"other steps", you 
would call VectorTools::interpolate_boundary_values() regardless 
stage to get constraints into the ConstraintMatrix object. If you do 
Newton-Raphson and solve for increments, at the first non-linear step you 
should have all constraints (hanging nodes + non-zero Dirichlet for a given 
step in quasi-static(?) loading), and for consequent steps only hanging 
nodes plus zero Dirichlet. 

Please read this module:
https://www.dealii.org/developer/doxygen/deal.II/group__constraints.html

This code-gallery may help in figuring out how to do 
this: 
https://www.dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html
 

 

> that transform the global matrix: columns of the constrained nodes become 
> zero columns with diagonal non-zero entries.
>
> In a time loop I have to account boundary condition without modifing 
> global matrix. While looping over dirichlet constrained cells
> I call the function  
> "hanging_node_constraints.distribute_local_to_global(cell_rhs, 
> local_dof_indices, system_rhs, cell_matrix)"
> which modifies the system_rhs vector by adding columns from constrained 
> nodes, which are actually zero in global matrix.
>
 
Nope, if your "hanging_node_constraints" indeed contains only constraints 
from hanging nodes, you won't see any contribution to the RHS
from it as constraints are homogeneous (ie. you wont' have  u_{123} = 456 
from hanging nodes).
 

>  
> Also I need to add components to rhs vector from the force boundary 
> conditions.
>

that's a different story, just do integration over faces and distribute as 
usual.
 

>
> However, I have to somehow account Dirichlet boundary conditions as in 
> previous step I didn't use any boundary values. The function
>

use ConstraintMatrix which has non-zero Dirichlet BC constraints.
 

> "interpolate_boundary_values(dof_handler, boundary_id, DirichletHandler, 
> boundary_values)"
> gives me "boundary_values" variables that maps node numbers to dirichle 
> boundary. 
>
> But I don't know how should I use it. What should I do with that variable 
> to account boundary conditions together with "distribute_local_to_global"?
>
> Thanks a lot in advance. Sorry for my misunderstanding.
>

Regards,
Denis.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: increment Dirichlet boundary condition

[Bryukhanov Ilya]

Denis, thanks a lot for the answer!

I think that on the first step I can use usual "interpolate_boundary_values 
+ apply_boundary_values" functions,
that transform the global matrix: columns of the constrained nodes become 
zero columns with diagonal non-zero entries.

In a time loop I have to account boundary condition without modifing global 
matrix. While looping over dirichlet constrained cells
I call the function  
"hanging_node_constraints.distribute_local_to_global(cell_rhs, 
local_dof_indices, system_rhs, cell_matrix)"
which modifies the system_rhs vector by adding columns from constrained 
nodes, which are actually zero in global matrix. 
Also I need to add components to rhs vector from the force boundary 
conditions.

However, I have to somehow account Dirichlet boundary conditions as in 
previous step I didn't use any boundary values. The function
"interpolate_boundary_values(dof_handler, boundary_id, DirichletHandler, 
boundary_values)"
gives me "boundary_values" variables that maps node numbers to dirichle 
boundary. 

But I don't know how should I use it. What should I do with that variable 
to account boundary conditions together with "distribute_local_to_global"?

Thanks a lot in advance. Sorry for my misunderstanding.

On Thursday, February 22, 2018 at 1:14:41 AM UTC+3, Denis Davydov wrote:
>
> Hi Iliya,
>
> there are numerous things that are wrong
>
> On Wednesday, February 21, 2018 at 5:32:55 PM UTC+1, Bryukhanov Ilya wrote:
>>
>> Denis, thanks a lot for your answer and links that you provided. 
>>
>> However, my solution becomes in a way that the displacement values on 
>> that boundary abnormally grow
>> even without assigning new boundary condition in a time loop. I attach 
>> figures on the zero step (correct) and on the first step (wrong).
>>
>> I do the following these things on each iteration: 
>> 1) recompute rhs vector - system_rhs (which is zero for my problem),
>> 2) compute matrix terms that are on the dirichlet boundary:
>> for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); 
>> ++cell) {
>> cell_matrix = 0; cell_rhs = 0; fe_values.reinit(cell);
>> for (unsigned int face_number=0; 
>> face_number> if (cell->face(face_number)->at_boundary()) {
>> // boundary index
>> int id = cell->face(face_number)->boundary_id();
>> if (boundary_displacement.find(id) != 
>> boundary_displacement.end()) {
>> //compute cell_matrix
>>
>> Then call function
>> hanging_node_constraints.distribute_local_to_global(cell_matrix, 
>> cell_rhs, local_dof_indices, system_matrix, system_rhs);
>>
>
> ^^ that should be hanging_node_constraints + Dirichlet 
>
>>
>> 3) Then  apply functions to account Dirichlet boundary condition:  
>>  VectorTools::interpolate_boundary_values(dof_handler, boundary_id, 
>> DirichletHandler, boundary_values);
>>
> ^^ that you would use during your setup
>  
>
>> MatrixTools::apply_boundary_values (boundary_values,
>> system_matrix,
>> solution,
>> system_rhs);
>>
>
> ^^ this is already what constraints.distribute_local_to_global() does, so 
> no need to do that.
>
> Regards
> Denis.
>  
>
>>
>>
>> Can you give some advice how to make it correct?
>>
>>
>> 
>>  
>> 
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Tuesday, February 20, 2018 at 10:51:17 PM UTC+3, Denis Davydov wrote:
>>>
>>> Hi Ilya,
>>>
>>> I am not away of principally simple solution, the RHS simply contains 
>>> your bilinear form times interpolated values at constrained DoFs. 
>>> What you can do is to loop over a subset of cells (those which are at 
>>> the boundary), still assemble the local matrix and do 
>>> constraints.distribute_local_to_global() which takes global matrix and 
>>> RHS.
>>>
>>> Or you can use constrained_linear_operator 
>>> ,
>>>  
>>> see this 
>>>  forum 
>>> discussion.
>>>
>>> Please also study this 
>>>  
>>> thread, which covers similar question.
>>>
>>> Regards,
>>> Denis.
>>>
>>> On Tuesday, February 20, 2018 at 4:46:19 PM UTC+1, Bryukhanov Ilya wrote:

 Hi,

 I consider elastic problem and forces and displacements on boundaries 
 are changing with time. I don't want to compute FEM matrix
 on each step. How to correctly simulate the task? I know that there is 
 a 

[deal.II] Re: increment Dirichlet boundary condition

[Bryukhanov Ilya]

Denis, thanks a lot for the answer!

I think that on the first step I can use usual "interpolate_boundary_values 
+ apply_boundary_values" functions,
that transform the global matrix: columns of the constrained nodes become 
zero columns with diagonal non-zero entries.

In a time loop I have to account boundary condition without modifing global 
matrix. While looping over dirichlet constrained cells
I call the function  
"hanging_node_constraints.distribute_local_to_global(cell_rhs, 
local_dof_indices, system_rhs, cell_matrix)"
which modifies the system_rhs vector by adding columns from constrained 
nodes, which are actually zero in global matrix. 

However, I have to somehow account boundary conditions. The function
"interpolate_boundary_values(dof_handler, boundary_id, DirichletHandler, 
boundary_values)"
gives me "boundary_values" variables that maps node numbers to dirichle 
boundary. 
But I don't know how should I use it. What should I do with that variable 
to account boundary conditions together with "distribute_local_to_global"?

Thanks a lot in advance. Sorry for my misundersting.

On Thursday, February 22, 2018 at 1:14:41 AM UTC+3, Denis Davydov wrote:
>
> Hi Iliya,
>
> there are numerous things that are wrong
>
> On Wednesday, February 21, 2018 at 5:32:55 PM UTC+1, Bryukhanov Ilya wrote:
>>
>> Denis, thanks a lot for your answer and links that you provided. 
>>
>> However, my solution becomes in a way that the displacement values on 
>> that boundary abnormally grow
>> even without assigning new boundary condition in a time loop. I attach 
>> figures on the zero step (correct) and on the first step (wrong).
>>
>> I do the following these things on each iteration: 
>> 1) recompute rhs vector - system_rhs (which is zero for my problem),
>> 2) compute matrix terms that are on the dirichlet boundary:
>> for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); 
>> ++cell) {
>> cell_matrix = 0; cell_rhs = 0; fe_values.reinit(cell);
>> for (unsigned int face_number=0; 
>> face_number> if (cell->face(face_number)->at_boundary()) {
>> // boundary index
>> int id = cell->face(face_number)->boundary_id();
>> if (boundary_displacement.find(id) != 
>> boundary_displacement.end()) {
>> //compute cell_matrix
>>
>> Then call function
>> hanging_node_constraints.distribute_local_to_global(cell_matrix, 
>> cell_rhs, local_dof_indices, system_matrix, system_rhs);
>>
>
> ^^ that should be hanging_node_constraints + Dirichlet 
>
>>
>> 3) Then  apply functions to account Dirichlet boundary condition:  
>>  VectorTools::interpolate_boundary_values(dof_handler, boundary_id, 
>> DirichletHandler, boundary_values);
>>
> ^^ that you would use during your setup
>  
>
>> MatrixTools::apply_boundary_values (boundary_values,
>> system_matrix,
>> solution,
>> system_rhs);
>>
>
> ^^ this is already what constraints.distribute_local_to_global() does, so 
> no need to do that.
>
> Regards
> Denis.
>  
>
>>
>>
>> Can you give some advice how to make it correct?
>>
>>
>> 
>>  
>> 
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Tuesday, February 20, 2018 at 10:51:17 PM UTC+3, Denis Davydov wrote:
>>>
>>> Hi Ilya,
>>>
>>> I am not away of principally simple solution, the RHS simply contains 
>>> your bilinear form times interpolated values at constrained DoFs. 
>>> What you can do is to loop over a subset of cells (those which are at 
>>> the boundary), still assemble the local matrix and do 
>>> constraints.distribute_local_to_global() which takes global matrix and 
>>> RHS.
>>>
>>> Or you can use constrained_linear_operator 
>>> ,
>>>  
>>> see this 
>>>  forum 
>>> discussion.
>>>
>>> Please also study this 
>>>  
>>> thread, which covers similar question.
>>>
>>> Regards,
>>> Denis.
>>>
>>> On Tuesday, February 20, 2018 at 4:46:19 PM UTC+1, Bryukhanov Ilya wrote:

 Hi,

 I consider elastic problem and forces and displacements on boundaries 
 are changing with time. I don't want to compute FEM matrix
 on each step. How to correctly simulate the task? I know that there is 
 a simple solution, but my solution is incorrect :)
 I have two variables for the system_matrix and rhs_vector. I try to 
 recompute rhs_vector on 

[deal.II] Re: increment Dirichlet boundary condition

[Bryukhanov Ilya]

Denis, thanks a lot for the answer! I think that on the first step I can 
use usual "interpolate_boundary_values + apply_boundary_values" functions,
and this procedures transform the global matrix so that the columns of 
constrained nodes become zero columns with diagonal non-zero entries.
In a time loop I have to account boundary condition so the calling of 
"hanging_node_constraints.distribute_local_to_global(cell_rhs, 
local_dof_indices, system_rhs, cell_matrix)"
while looping over dirichlet constrained cells modify the system_rhs vector 
by adding columns from constrained nodes,
which are actually zero in global matrix. However, I have to somehow 
account boundary conditions. The function
"interpolate_boundary_values(dof_handler, boundary_id, DirichletHandler, 
boundary_values)"
gives me "boundary_values" variables that maps node numbers to dirichle 
boundary. But I don't know how should I use it. What should I 
do with the variable "boundary_values" to account boundary conditions in 
"distribute_local_to_global"?

Thanks a lot in advance. Sorry for my misundersting.

On Thursday, February 22, 2018 at 1:14:41 AM UTC+3, Denis Davydov wrote:
>
> Hi Iliya,
>
> there are numerous things that are wrong
>
> On Wednesday, February 21, 2018 at 5:32:55 PM UTC+1, Bryukhanov Ilya wrote:
>>
>> Denis, thanks a lot for your answer and links that you provided. 
>>
>> However, my solution becomes in a way that the displacement values on 
>> that boundary abnormally grow
>> even without assigning new boundary condition in a time loop. I attach 
>> figures on the zero step (correct) and on the first step (wrong).
>>
>> I do the following these things on each iteration: 
>> 1) recompute rhs vector - system_rhs (which is zero for my problem),
>> 2) compute matrix terms that are on the dirichlet boundary:
>> for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); 
>> ++cell) {
>> cell_matrix = 0; cell_rhs = 0; fe_values.reinit(cell);
>> for (unsigned int face_number=0; 
>> face_number> if (cell->face(face_number)->at_boundary()) {
>> // boundary index
>> int id = cell->face(face_number)->boundary_id();
>> if (boundary_displacement.find(id) != 
>> boundary_displacement.end()) {
>> //compute cell_matrix
>>
>> Then call function
>> hanging_node_constraints.distribute_local_to_global(cell_matrix, 
>> cell_rhs, local_dof_indices, system_matrix, system_rhs);
>>
>
> ^^ that should be hanging_node_constraints + Dirichlet 
>
>>
>> 3) Then  apply functions to account Dirichlet boundary condition:  
>>  VectorTools::interpolate_boundary_values(dof_handler, boundary_id, 
>> DirichletHandler, boundary_values);
>>
> ^^ that you would use during your setup
>  
>
>> MatrixTools::apply_boundary_values (boundary_values,
>> system_matrix,
>> solution,
>> system_rhs);
>>
>
> ^^ this is already what constraints.distribute_local_to_global() does, so 
> no need to do that.
>
> Regards
> Denis.
>  
>
>>
>>
>> Can you give some advice how to make it correct?
>>
>>
>> 
>>  
>> 
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Tuesday, February 20, 2018 at 10:51:17 PM UTC+3, Denis Davydov wrote:
>>>
>>> Hi Ilya,
>>>
>>> I am not away of principally simple solution, the RHS simply contains 
>>> your bilinear form times interpolated values at constrained DoFs. 
>>> What you can do is to loop over a subset of cells (those which are at 
>>> the boundary), still assemble the local matrix and do 
>>> constraints.distribute_local_to_global() which takes global matrix and 
>>> RHS.
>>>
>>> Or you can use constrained_linear_operator 
>>> ,
>>>  
>>> see this 
>>>  forum 
>>> discussion.
>>>
>>> Please also study this 
>>>  
>>> thread, which covers similar question.
>>>
>>> Regards,
>>> Denis.
>>>
>>> On Tuesday, February 20, 2018 at 4:46:19 PM UTC+1, Bryukhanov Ilya wrote:

 Hi,

 I consider elastic problem and forces and displacements on boundaries 
 are changing with time. I don't want to compute FEM matrix
 on each step. How to correctly simulate the task? I know that there is 
 a simple solution, but my solution is incorrect :)
 I have two variables for the system_matrix and rhs_vector. I try to 
 recompute rhs_vector on each step  and