Dirichlet BC on block matrices

[Мария Бронзова]

Deal all,

I have an issue implementing Dirichlet boundary conditions in a mixed pressure-displacement problem. 

In step-20 there is an example for implementing pressure boundary conditions like that:

template <int dim>

class PressureBoundaryValues : public Function<dim>

{

public:

PressureBoundaryValues()

: Function<dim>(1)

{}

virtual double value(const Point<dim> & p,

const unsigned int component = 0) const override;

};

template <int dim>

double RightHandSide<dim>::value(const Point<dim> & /*p*/,

const unsigned int /*component*/) const

{

return 0;

}

And later:

for (const auto &face : cell->face_iterators())

if (face->at_boundary())

{

fe_face_values.reinit(cell, face);

pressure_boundary_values.value_list(

fe_face_values.get_quadrature_points(), boundary_values);

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

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

local_rhs(i) += -(fe_face_values[velocities].value(i, q) *

fe_face_values.normal_vector(q) *

boundary_values[q] *

fe_face_values.JxW(q));

}

I was trying to follow the same logic, implementing displacement BC (want to set the displacement of certain nodes to zero). In this case I create the Boundary condition template as follows:

template <int dim>

      class DisplBoundaryValues : public Function<dim, std::complex<double>>

      {

      public:

      DisplBoundaryValues()

        : Function<dim, std::complex<double>>(dim)

          {}

      virtual std::complex<double> value(const Point<dim> &p,

               const unsigned int component = 0) const override;

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

      Vector<std::complex<double>> &  value) const override;

      };

    template <int dim>

  std::complex<double> DisplBoundaryValues<dim>::value(const Point<dim> & p,

     const unsigned int component) const

{

   const std::complex<double> i = {0,1};

   return 0.+0.*i;

}

    template <int dim>

    void DisplBoundaryValues<dim>::vector_value(const Point<dim> &p,

                                           Vector<std::complex<double>> & values) const

    {

      for (unsigned int c = 0; c < this->n_components; ++c)

      {

      const std::complex<double> i = {0,1};

      values(c) = 0.+0.*i;

      }

    }

After that I am trying to apply the boundary values to the right hand side. This time i should use the vector_value_list instead of value_list to get the boundary values assigned. But I am confused about the further proceeding while implementing the locar_rhs part. How can I apply these boundary values in vector form so that the algorithm understands these values need to be assigned exactly to the displacement components? 

for (const auto &face : cell->face_iterators())

  if (face->boundary_id()==2)

  {

  displ_boundary_values.vector_value_list(

  fe_face_values.get_quadrature_points(), d_boundary_values);

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

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

  {

  std::cout << d_boundary_values[q]<< "  ";

  local_rhs(i) += -porosity*(1.+Q[k]/R[k])

  *d_boundary_values[q]

  *fe_face_values[pressure].value(i,q)

  *fe_face_values.JxW(q);

  }

  }

Than you a lot for your time!

Kind regards,

Mariia Bronzova

[Wolfgang Bangerth]

On 9/28/21 8:28 AM, Мария Бронзова wrote:
>
> After that I am trying to apply the boundary values to the right hand

> side. This time i should use the *vector_value_list *instead of

> value_list to get the boundary values assigned. But I am confused about
> the further proceeding while implementing the locar_rhs part. How can I
> apply these boundary values in vector form so that the algorithm
> understands these values need to be assigned exactly to the displacement
> components?

The first step is always to write out what the weak formulation should
be. In your case, there will probably be a boundary integral term that
includes \vec g, the boundary conditions for your vector-valued
quantity. But how does the rest of the term look like?

Best
W.

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

[Мария Бронзова]

Dear Wolfgang,

Thank you a lot for the reply! I am trying to implement the Biot equations for the propagation of acoustic waves in porous media. The equations are attached in the "Biot .PNG" file. There are two boundary integrals there, containing information on work, done by external forces on the fluid and skeleton part of the media. I want to implement fixed displacement boundary conditions ("BC.PNG")  on a part of the geometry with ID = 2, so I need to set the displacement components to zero. I was wondering, whether there is a way to do it similarly to the way it"s done in step-20 example for the pressure boundary condition (so, in the local_rhs vector).

Please, let me know, if it works so.

Thank you for your support!

With kind regards,

Mariia Bronzova

среда, 29 сентября 2021 г. в 19:23:16 UTC+2, Wolfgang Bangerth:

[Wolfgang Bangerth]

On 9/30/21 2:48 AM, Мария Бронзова wrote:
>
> Thank you a lot for the reply! I am trying to implement the Biot equations for
> the propagation of acoustic waves in porous media. The equations are attached
> in the "Biot .PNG" file. There are two boundary integrals there, containing
> information on work, done by external forces on the fluid and skeleton part of
> the media. I want to implement fixed displacement boundary conditions
> ("BC.PNG")  on a part of the geometry with ID = 2, so I need to set the
> displacement components to zero. I was wondering, whether there is a way to do
> it similarly to the way it"s done in step-20 example for the pressure boundary
> condition (so, in the local_rhs vector).

Maria: so which of the terms in that formula are you interested in implementing?

[Мария Бронзова]

So, there are two boundary integrals in the formulation and I am trying to implement them for the case of fixed displacement boundary condition.  The first boundary integral falls to zero for such a BC, as no displacement variation is possible in this case. So, I am implementing the second integral from the second equation. The integral can be represented as in the BC.PNG file. There are those factors in brackets, assigned displacement values and variation of the pressure multiplied together:

  local_rhs(i) += -porosity*(1.+Q[k]/R[k])

  *d_boundary_values[q]

  *fe_face_values[pressure].value(i,q)

  *fe_face_values.JxW(q);

But the way it is written now it cannot work, as the d_boundary_values term is a vector of vectors (as we have three displacement components). So, the question is, whether there is a way to treat the displacement components seperately in this d_boundary_values term? Or maybe even a smarter way to make it work?

Please, let me know. Thank you a lot!

четверг, 30 сентября 2021 г. в 14:28:26 UTC+2, Wolfgang Bangerth:

[Wolfgang Bangerth]

On 9/30/21 6:49 AM, Мария Бронзова wrote:
> So, there are two boundary integrals in the formulation and I am trying to
> implement them for the case of fixed displacement boundary condition.  The
> first boundary integral falls to zero for such a BC, as no displacement
> variation is possible in this case. So, I am implementing the second integral
> from the second equation. The integral can be represented as in the BC.PNG
> file. There are those factors in brackets, assigned displacement values and
> variation of the pressure multiplied together:
>

> *local_rhs*(i) += -porosity*(1.+Q[k]/R[k])

>   *d_boundary_values[q]
>   *fe_face_values[pressure].value(i,q)
>   *fe_face_values.JxW(q);
>

> But the way it is written now it cannot work, as the *d_boundary_values* term

> is a vector of vectors (as we have three displacement components). So, the
> question is, whether there is a way to treat the displacement components

> seperately in this *d_boundary_values* term? Or maybe even a smarter way to
> make it work?

I think that your question is actually of mathematical nature, not one of
implementation. If I read the integral I_2 correctly in your previous email,
then what you prescribe there is
u^i_n
which I believe is not actually the displacement on the boundary (a vector)
but only the *normal component* of the velocity (a scalar). So you have two
options:
- You write a function that only returns the normal velocity (which is all
you can prescribe anyway)
- You write a function that returns the velocity at the boundary as a vector
and then in the bilinear form, you take the dot product with the normal
vector (which you can get from the fe_face_values object).

Both are reasonable, though if all you can prescribe is the normal component,
you might as well write your function in such a way that that is what it returns.

[Мария Бронзова]

Dear Wolfgang,

I had indeed a gap in understanding of some points, it is solved now! I thank you a lot for your time and effort in trying to understand my problem and to help me out! 

With kind regards,

Mariia

среда, 6 октября 2021 г. в 03:02:52 UTC+2, Wolfgang Bangerth: