Step 13 & 14 -- Different BCs on different boundary indicators

[Elena Greco]

Hello,

In Steps 13 and 14, for primal problem, the boundary condition is calculate based on an analytical solution, therefore it is easy to set different boundary indicators and apply the BC on different boundary indicators.

In my problem I have a cylinder (boundary indicator of R_in is 0 and boundary indicator of R_out is 1).
How can I implement two different BC on these BCs?
for example on R_in the displacement U_x = x + y and on R_out displacement U_x = x^2 + y.

I don't want to destroy the extensibility of the code, therefore I would like to modify template <int dim> double CurvedRidges<dim>::BoundaryValues::value (const Point<dim>   &p, const unsigned int  /*component*/) const {...} part NOT the Solver<dim>::assemble_linear_system part.

I was thinking of passing boundary indicator to the boundary function or something similar; like:

VectorTools::interpolate_boundary_values (dof_handler,
                                                0,
                                                *boundary_values(0),
                                                boundary_value_map);

and then relate the BC in CurvedRidges<dim>::BoundaryValues::value(...) to boundary indicator, but I am not sure if it is possible and which parts should I modify.

Thanks
Elena

[Uwe Köcher]

Elena,

in your case you do not really need a second boundary indicator.

You have the position argument "Point<dim> x" in the value function and

from that you can compute the radius and decide if the first or second

boundary condition holds and return the respective value.

Then you have not to change any other parts of the example to include that.

Kind regards

  Uwe

[Elena Greco]

Dear Uwe

Thanks for you prompt reply.

But if I import a mesh from gmesh and for a not regular geometry I would face a problem.
Do you have any suggestion for this case?

Thanks.
Elena

On Tuesday, May 16, 2017 at 12:50:36 PM UTC+2, Elena Greco wrote:

[Elena Greco]

May be it would be better to edit my question and remove the "cylinder" part, but I do not have access to modify the question.
Please consider my question for a general case of a geometry with multiple boundary indicators and different BCs on them.

Thanks.
Elena

[Bruno Turcksin]

Elena,

On Tuesday, May 16, 2017 at 7:19:08 AM UTC-4, Elena Greco wrote:

May be it would be better to edit my question and remove the "cylinder" part, but I do not have access to modify the question.
Please consider my question for a general case of a geometry with multiple boundary indicators and different BCs on them.

You should first tag the boundaries with boundary indicators and then use http://dealii.org/8.5.0/doxygen/deal.II/group__constraints.html#gaa96810fcf5053e0bf08785aa53bfd6cc The map associates a boundary ID to a function and you can have as many functions as you want.

Best,

Bruno

[Elena Greco]

Dear Bruno
Thanks for your reply.

I know how to use interpolate_boundary_values but in steps 13 & 14, if I want to keep Solver<dim>::assemble_linear_system in the general form, I don't know how to implement new boundary conditions.

In steps 13 & 14, Solver<dim>::assemble_linear_system assembles the linear system and implements BC for primal and dual problems of several cases like CurvedRidges, Exercise_2_3, and three new cases which I have added. I do not want to change Solver<dim>::assemble_linear_system, and I want to keep it in general form and I only want to modify template <int dim> double My_case_2<dim>::BoundaryValues::value(...)

Thanks.
Elena

[Bruno Turcksin]

Elena,

I am not sure I understand exactly what you want to do. Can you write
a simple pseudo code to show the problem? The way I understand what
you want to do, you need to replace interpolate_boundary_values
function from step 13 & 14 (which allows you to use only one BC) and
instead used the one that I showed you which uses a map. Yous also
need to add a function that fills this map and where you do the switch
between the different cases. This way assemble_linear_system is still
general and everything is encapsulated in the filling function.

Best,

Bruno

> --
> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/Nt9BTAyX6Pk/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Elena Greco]

Dear Bruno,

I would like to somehow pass the boundary indicator to the boundary value function.
I don't know if it is possible to write *boundary_values(0) or not?
I want to do something like below code, but I don't know how to write it correctly:

template <int dim>
    void
    Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
    {
        ....

        VectorTools::interpolate_boundary_values (dof_handler,
                                                                      0,
                                                                     *boundary_values(0),
                                                                     boundary_value_map);

        VectorTools::interpolate_boundary_values (dof_handler,
                                                                      1,
                                                                      *boundary_values(1),
                                                                      boundary_value_map);
        ... 
    } 


template<int dim>
    void My_Case_2<dim>::BoundaryValues::vector_value(const Point<dim> &p,
                                                                                    Vector<double> &values) const
    {
        ...
        if (boundary_indicator == 0) {
           BC = p(0) + p(1);
        } else if (boundary_indicator == 1) {
           BC = p(0)*p(0) + p(1);
        }
        ...
    }
             
Thanks.
Elena

[Daniel Arndt]

Dear Elena,

I want to do something like below code, but I don't know how to write it correctly:

template <int dim>
    void
    Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
    {
        ....
        VectorTools::interpolate_boundary_values (dof_handler,
                                                                      0,
                                                                     *boundary_values(0),
                                                                     boundary_value_map);

        VectorTools::interpolate_boundary_values (dof_handler,
                                                                      1,
                                                                      *boundary_values(1),
                                                                      boundary_value_map);
        ... 
    } 


template<int dim>
    void My_Case_2<dim>::BoundaryValues::vector_value(const Point<dim> &p,
                                                                                    Vector<double> &values) const
    {
        ...
        if (boundary_indicator == 0) {
           BC = p(0) + p(1);
        } else if (boundary_indicator == 1) {
           BC = p(0)*p(0) + p(1);
        }
        ...
    }

You can do something like this if you define a member variable current_boundary_indicator and set it correctly before calling interpolate_boundary_values(), i.e. you would do

boundary_values.current_boundary_indicator=0;

VectorTools::interpolate_boundary_values (dof_handler,0,boundary_values, boundary_value_map);

boundary_values.current_boundary_indicator=1; 

VectorTools::interpolate_boundary_values (dof_handler,1,boundary_values, boundary_value_map);

 

If you have a look at the documentation of VectorTools::interpolate_boundary_value [1], you see that there are multiple functions with this name.

In particular, the one with the signature

void VectorTools::interpolate_boundary_values

  (const DoFHandlerType< dim, spacedim > &                                                       dof,

   const std::map< types::boundary_id, const Function< spacedim, number > * > & function_map,

   std::map< types::global_dof_index, number > &                                                 boundary_values,

   const ComponentMask &                                                                                 component_mask = ComponentMask() )

might be a good alternative for you as also Bruno (with a ConstraintMatrix instead of std::map< types::global_dof_index, number > &) already pointed out.

Then, you would define one Function for each boundary_indicator and call VectorTools::interpolate_boundary_values only once.

Best,

Daniel

[1] https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceVectorTools.html#a9f3e3ae1396811f998cc35f94cbaa926