Assemble system with variational right hand side

[Jaekwang Kim]

Hi. I have a question on how we can impose external force term, or right hand side, in Stokes flow. 

My governing equation is same with that of step-22, stokes flow 

-nabla^2 u + grad p = f 

-grad u = 0 

but I want to impose 'f', 

To do this, I made function that gives me f value on each point. 

For example, if I want function of f thats looks like f=(x^2,y^2); 

I would make function like this...

  template <int dim>

  class RightHandSide : public Function<dim>

  {

    public:

      RightHandSide () : Function<dim>(dim+1) {}

      virtual void vector_value (const Point<dim> &p,

Vector<double>   &value) const;

  };

  template <int dim>

  void

  RightHandSide<dim>::vector_value (const Point<dim> &p,

    Vector<double>   &values) const

  {

    values(0) =  pow(p[0],2);

    values(1) =  pow(p[1],2);

    values(2) =  0.0;

  }

After that I may consider this 'f' term when I assemble system. 

but I don't see a way to how I should try modify following read lines of the code from step-22 

  for (; cell!=endc; ++cell)

      {

fe_values.reinit (cell);

local_matrix = 0;

local_rhs = 0;

right_hand_side.vector_value_list(fe_values.get_quadrature_points(),

  rhs_values);

for (unsigned int q=0; q<n_q_points; ++q)

  {

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

      {

phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);

div_phi_u[k]   = fe_values[velocities].divergence (k, q);

phi_p[k]       = fe_values[pressure].value (k, q);

      }

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

      {

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

  {

    local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]

  - div_phi_u[i] * phi_p[j]

  - phi_p[i] * div_phi_u[j]

  + phi_p[i] * phi_p[j])

* fe_values.JxW(q);

  }

const unsigned int component_i =

  fe.system_to_component_index(i).first;

local_rhs(i) += fe_values.shape_value(i,q) *

rhs_values[q](component_i) *

fe_values.JxW(q);

      }

  }

My 'f' is not '0' any more but it has multiple components. 

in fact, the weak form of the governing equation '-nabla^2 u + grad p = f ' also should have two components,..(x and y) ....

How can I achieve this....?

Thank you 

[Daniel Arndt]

Jaekwang,

 

[...]

My governing equation is same with that of step-22, stokes flow

-nabla^2 u + grad p = f

-grad u = 0

but I want to impose 'f',

To do this, I made function that gives me f value on each point. 

For example, if I want function of f thats looks like f=(x^2,y^2); 

I would make function like this...

[...]

After that I may consider this 'f' term when I assemble system. 

but I don't see a way to how I should try modify following read lines of the code from step-22 

[...]

const unsigned int component_i = fe.system_to_component_index(i).first;
local_rhs(i) += fe_values.shape_value(i,q) * rhs_values[q](component_i) * fe_values.JxW(q);  

[...]

My 'f' is not '0' any more but it has multiple components. 

in fact, the weak form of the governing equation '-nabla^2 u + grad p = f ' also should have two components,..(x and y) ....

What you are doing already looks correct. In fact, this is how the right-hand side is treated in step-22 and is also suitable for a non-zero 'f'.
The only thing you need to change is the implementation of 'RightHandSide' and that looks okay as well (provided you are only considering dim=2).
What is the problem you are observing?

Best,
Daniel

[Jaekwang Kim]

Thanks for the reply, 

I just wanted to make sure that whether I implemented it right, 

but I still cannot sure whether I am doing correctly..... Let me be with more specific questions. 

First, Yes I am considering 2-dimensional flow.

if I writes a code that looks as previous, 

template <int dim>

  void

  RightHandSide<dim>::vector_value (const Point<dim> &p,

    Vector<double>   &values) const

  {

    values(0) =  pow(p[0],2);  

    values(1) =  pow(p[1],2);

    values(2) =  0.0;  // I assumed this always zero. because f is 2 dimensional force, and we will not have third component of f. 

  }

1) whether my interpretation of 'values(2) is correct' - marked with red line 

2) If we considered a vector value test function phi=(v,q)= (v_x, v_y , q), the right hand side term in weak form should look like this 

RHS = (v, f ) at domain o mega. 

since both v and f is vector that has two direction , (x and y) , the calculation should be conducted in both direction 

but in the following code, it seems that I am multiplying only the first component v and f... 

Do I understand this appropriately?

In this sense, I don't understand my code well... cuz it looks like...

const unsigned int component_i =

  fe.system_to_component_index(i).first;   <- Here I only said the first component...

local_rhs(i) += fe_values.shape_value(i,q) *

rhs_values[q](component_i) *

fe_values.JxW(q);

Shouldn't this be some thing like this?  

const unsigned int component_i =

  fe.system_to_component_index(i).first; 

const unsigned int component_j = 

                  fe.system_to_component_index(j).second; 

local_rhs(i) += (fe_values.shape_value(i,q) *rhs_values[q](component_i)) 

                                (fe_values.shape_value(i,q) *rhs_values[q](component_j)) 

fe_values.JxW(q); 

[Bruno Turcksin]

Hi,

On Monday, January 23, 2017 at 2:06:17 PM UTC-5, Jaekwang Kim wrote:

I just wanted to make sure that whether I implemented it right, 

but I still cannot sure whether I am doing correctly..... Let me be with more specific questions.

This is not the right way to do it. You cannot write your code and then, ask people if what you did is correct. Maybe it will look correct but something that you don't show is wrong or the person looking at your code misses something obvious. On top of that if everybody on the mailing list starts asking if their code is correct, we won't have time to answer other questions. What you need to do is to write a test, i.e. a small code that use the function you wrote, and you need to make sure that the code works like it should on this simple example. This way not only you will be able to convince yourself that the code works like it should but you can make sure six months from now that everything is still working. If your test doesn't work like expected, then you may want to ask questions on the mailing list.

Best,

Bruno

[Jaekwang Kim]

you are right. 

Sorry about that. 

I may test and learn more by myself. 

2017년 1월 23일 월요일 오후 1시 21분 56초 UTC-6, Bruno Turcksin 님의 말: