Assemble the rhs of time-dependent linear elasticity

[Jie Cheng]

Hi everyone

I have a time-dependent linear elasticity solver which is different from step-17 in which I assume small-deformation all the time. Due to the time discretization, I have a -K*u^n term at the rhs of the system where u^n is the known displacement (at previous time step) and K is the stiffness matrix. I tried two different approaches to deal with this term:

1. The first approach is simple: I assemble the system_rhs as if there is no such term, and then explicitly compute system_stiffness.vmult(previous_displacement), and

subtract the result from the system_rhs. It worked out nicely and I got my time-dependent results.

2. In the second approach I deal with this term in the assembly, namely I add this term to the local_rhs at cell level

Loop over cells:

  local_matrix = 0; // this is stiffness + mass term

  local_stiffness = 0;

  local_rhs = 0;

  Loop over quadrature points:

      compute local_matrix

      compute local_stiffness

      compute local_rhs // here I compute the rhs due to external load

  end of quadrature loop

  // Now compute the rhs due to time-discretization

  cell->get_dof_indices(local_dof_indices);

  for (unsigned int i = 0; i < dofs_per_cell; ++i)

  {

    for (unsigned int j = 0; j < dofs_per_cell; ++j)

    {

       local_rhs[i] -= local_stiffness[i][j] * dt * previous_displacement[local_dof_indices[j]];

    }

  }  

  constraints.distribute_local_to_global(local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs);

end of cell loop

Because the global stiffness matrix is the sum of local matrices, ideally this approach should work as the first one. But the solution turned out to be garbage. Could anybody see why the second approach is wrong?

Thank you

Jie

      

      

[Wolfgang Bangerth]

On 11/28/2017 02:33 PM, Jie Cheng wrote:
>
> Because the global stiffness matrix is the sum of local matrices, ideally this
> approach should work as the first one. But the solution turned out to be
> garbage. Could anybody see why the second approach is wrong?

Can you be more specific what "garbage" means? Are the boundary values wrong?
The difference between your two methods is what happens to non-zero boundary
values in the rhs vector.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jie Cheng]

Hi Wolfgang

Thanks for your reply. I carefully tested these two methods and found that they produce two different rhs vectors given identical previous solution. There is no non-zero boundary values in my case, and the constrained values are indeed zero. It seems that the unconstrained components in the rhs vector are wrong. 

Here is what I suspect: Method 1 uses the stiffness_matrix that has been modified during the assembly, while method 2 is using the original stiffness. Although method 2 will eventually distribute the rhs vector as well as the matrix, it is distributed along with system_matrix, not stiffness_matrix. In the time-dependent simulation, those two are different. I don't know if the modifications to the rhs vector depends on the matrix or not.

Thanks

Jie

[Wolfgang Bangerth]

On 11/29/2017 09:39 PM, Jie Cheng wrote:
>
> Here is what I suspect: Method 1 uses the stiffness_matrix that has been
> modified during the assembly, while method 2 is using the original
> stiffness. Although method 2 will eventually distribute the rhs vector
> as well as the matrix, it is distributed along with system_matrix, not
> stiffness_matrix. In the time-dependent simulation, those two are
> different. I don't know if the modifications to the rhs vector depends
> on the matrix or not.

For hanging nodes they do. For non-zero boundary values the same applies
(see video lectures 21.6 and 21.65 to understand why).

I would have to think for a while to figure out what exactly is going on
in your case, and how to make sure you get the same results, but don't
have the necessary time for this. I think that if you look through the
video lectures references above, you'll see what one of the two methods
you use is doing, and that may help you understand what you need to do
to get things right.