Using interpolate_boundary_conditions

[Amy Kaczmarowski]

Hi,

I'm attempting to apply a displacement boundary condition on a surface in my problem.  However, I would like the direction of the displacement of each point on the surface to depend on its direction from some point (for example the origin).  Imagine a balloon inflating where the displacement direction is from the center of the balloon.  My thought was to use a function similar to the example boundary conditions in step 23 but extended for 3 dimensions.  I defined the following class in my code:

template <int dim>

class BoundaryValuesU : public Function<dim>{

  public:

    virtual void vector_value(const Point<dim> & p, Vector<double> &values) const override

      {

        double xpos, ypos, zpos;

        xpos = p[0];

        ypos = p[1];

        if (dim==3){zpos=p[2];}

        else {zpos=0.0;}

        values[0] = xpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        values[1] = ypos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        values[2] = zpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        return; 

    }

  };

Then I apply the boundary conditions in the following way to the surface with boundary id 2:

std::vector<bool> uBC (3, true);

BoundaryValuesU<dim> boundary_values_u_function;

VectorTools::interpolate_boundary_values (dof_handler, 2, boundary_values_u_function, constraints, uBC);

However, when I attempt to run with this case I get the following error:

An error occurred in line <2263> of file <src/deal.ii-candi/tmp/unpack/deal.II/include/deal.II/numerics/vector_tools.templates.h> in function

    void dealii::VectorTools::<unnamed>::do_interpolate_boundary_values(const M_or_MC<dim, spacedim> &, const DoFHandlerType<dim, spacedim> &, const std::map<unsigned char, const dealii::Function<spacedim, number> *, std::less<unsigned char>, std::allocator<std::pair<const unsigned char, const dealii::Function<spacedim, number> *>>> &, std::map<unsigned int, number, std::less<unsigned int>, std::allocator<std::pair<const unsigned int, number>>> &, const dealii::ComponentMask &) [with dim = 3, spacedim = 3, number = double, DoFHandlerType = dealii::DoFHandler, M_or_MC = dealii::Mapping]

The violated condition was:

    n_components == i->second->n_components

Additional information:

    Dimension 3 not equal to 1.

Is there a way I am supposed to be initializing n_components in the class to make this work?  Or is this just the wrong way to be doing this?

Thank you.

[Paras Kumar]

Hi Amy,

 I post below one possible solution which worked for me.

On Tuesday, October 1, 2019 at 4:44:25 AM UTC+2, Amy Kaczmarowski wrote:

Hi,

I'm attempting to apply a displacement boundary condition on a surface in my problem.  However, I would like the direction of the displacement of each point on the surface to depend on its direction from some point (for example the origin).  Imagine a balloon inflating where the displacement direction is from the center of the balloon.  My thought was to use a function similar to the example boundary conditions in step 23 but extended for 3 dimensions.  I defined the following class in my code:

template <int dim>

class BoundaryValuesU : public Function<dim>{

  public:

    virtual void vector_value(const Point<dim> & p, Vector<double> &values) const override

      {

        double xpos, ypos, zpos;

        xpos = p[0];

        ypos = p[1];

        if (dim==3){zpos=p[2];}

        else {zpos=0.0;}

        values[0] = xpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        values[1] = ypos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        values[2] = zpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;

        return; 

    }

  };

Then I apply the boundary conditions in the following way to the surface with boundary id 2:

std::vector<bool> uBC (3, true);

BoundaryValuesU<dim> boundary_values_u_function;

VectorTools::interpolate_boundary_values (dof_handler, 2, boundary_values_u_function, constraints, uBC);

 

const dealii::FEValuesExtractors::Vector & dispExtractor(0); // assuming we only have a vector field (having dim components) as unknown at each support point
 VectorTools::interpolate_boundary_values (dof_handler, 2, boundary_values_u_function, constraints, feSystem.component_mask(dispExtractor));

 

Here, the function called on the feSystem object is documented at https://www.dealii.org/8.5.0/doxygen/deal.II/classFiniteElement.html#a530f6c3f6326e516aeb006e7533a4532.

However, when I attempt to run with this case I get the following error:

An error occurred in line <2263> of file <src/deal.ii-candi/tmp/unpack/deal.II/include/deal.II/numerics/vector_tools.templates.h> in function

    void dealii::VectorTools::<unnamed>::do_interpolate_boundary_values(const M_or_MC<dim, spacedim> &, const DoFHandlerType<dim, spacedim> &, const std::map<unsigned char, const dealii::Function<spacedim, number> *, std::less<unsigned char>, std::allocator<std::pair<const unsigned char, const dealii::Function<spacedim, number> *>>> &, std::map<unsigned int, number, std::less<unsigned int>, std::allocator<std::pair<const unsigned int, number>>> &, const dealii::ComponentMask &) [with dim = 3, spacedim = 3, number = double, DoFHandlerType = dealii::DoFHandler, M_or_MC = dealii::Mapping]

The violated condition was:

    n_components == i->second->n_components

Additional information:

    Dimension 3 not equal to 1.

Is there a way I am supposed to be initializing n_components in the class to make this work?  Or is this just the wrong way to be doing this?

Thank you.

Best regards,

Paras

[Wolfgang Bangerth]

Amy,

> template <int dim>
> class BoundaryValuesU : public Function<dim>{
>   public:

You are correct: You need to tell the base class how many vector components
your function has. You should be able to do this by adding a constructor:

BoundaryValuesU ()
: Function<dim> (dim) // your function has 'dim' components
{}

>     virtual void vector_value(const Point<dim> & p, Vector<double> &values)
> const override
>       {
>         double xpos, ypos, zpos;
>         xpos = p[0];
>         ypos = p[1];
>         if (dim==3){zpos=p[2];}
>         else {zpos=0.0;}
>         values[0]
> = xpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;
>         values[1]
> = ypos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;
>         values[2]
> = zpos/std::pow(std::pow(xpos,2.0)+std::pow(ypos,2.0)+std::pow(zpos,2.0),0.5)*0.1;
>         return;
>     }

You can write this much simpler using the fact that Point objects are really
just vectors:

virtual void vector_value (...) const override
{
Tensor<1,dim> displacement = 0.1 * p / p.norm();
for (unsigned int d=0; d<dim; ++d)
values[d] = displacement[d];
}

Best
W.

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

[Amy Kaczmarowski]

Perfect!  Thank you this was exactly what I was missing.

Amy