Non Homogenous Boundary Conditions on a Cylinder

[Samuel Rodriguez]

Hello Everyone,

I am a masters student barely getting into working with the very heavy theoretical foundations of using deal.ii. I have been working on a research project that uses a code that creates a cylinder with Dirichlet boundary conditions on the hull/bottom faces and a Neumann boundary condition on the top face Recently I have been trying to change the Dirichlet boundary condition on the hull of the cylinder to a Neumann boundary condition. Our code can be found here inside of QuasistaticBrownianThermalNoise.cpp:

Numerical Coating Thermal Noise

I have applied the boundary condition but it does not seem to be working. So what I did instead was go to the examples in the tutorial section to learn ow to apply non-homogenous boundary conditions on a cylinder. I have solved the laplacian on the this cylinder using step 3. I have tried to use step 7 to implement non-homogeneous boundary conditions but step 7 requires using a known solution. I would like to be able to set the derivative of the function on the hull to zero and the derivative on the top face to some known function. A Dirichlet boundary condition would be set on the bottom face. The assemble function is down below. My question is how exactly do I apply the Neumann boundary condition to the faces of the cylinder. In the code below I have comments that explain my question in more detail. My question begins where I loop over the faces of the cylinder. However, I included the entire assembly function for completeness in case it was needed. I also attached the whole cpp file in case that was also needed. Any help would be greatly appreciated. If I can provide any more information I would be happy to do so.

QGauss<3> quadrature_formula(fe.degree + 1);

QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);

FEValues<3> fe_values(fe,

quadrature_formula,

update_values | update_gradients | update_JxW_values);

FEFaceValues<3> fe_face_values(fe,

face_quadrature_formula,

update_values | update_quadrature_points |

update_normal_vectors |

update_JxW_values);

const unsigned int dofs_per_cell = fe.dofs_per_cell;

const unsigned int n_q_points = quadrature_formula.size();

const unsigned int n_face_q_points = face_quadrature_formula.size();

FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);

Vector<double> cell_rhs(dofs_per_cell);

std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

for (const auto &cell : dof_handler.active_cell_iterators())

{

fe_values.reinit(cell);

cell_matrix = 0;

cell_rhs = 0;

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

{

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

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

cell_matrix(i, j) +=

(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)

fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)

fe_values.JxW(q_index)); // dx

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

cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)

1 * // f(x_q)

fe_values.JxW(q_index)); // dx

}

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)

system_matrix.add(local_dof_indices[i],

local_dof_indices[j],

cell_matrix(i, j));

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

system_rhs(local_dof_indices[i]) += cell_rhs(i);

// I know that the hull of the cylinder has boundary id 0, the top face has boundary id 1, and the bottom of the face has boundary id 2.

for(unsigned int face_number = 0; face_number < GeometryInfo<3>::faces_per_cell; ++face_number)

{

if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->boundary_id() == 0))

{

fe_face_values.reinit(cell, face_number);

// If we come in here we have found a face that belongs to the boundary condtion of the hull

// I know that I am supposed to do something like the code in green below, but I don't know the exact solution.

// What I would like to do is set the derivative of my function to zero. My thinking is that it would entail

// taking the gradient of fe_face_values to ZeroFunction<3>(). I think that if I could understand how to apply

// the boundary condition to the hull of the cylinder, I could understand how to apply the the boundary condition

// to the top of the face just as easily. Here is the code:

for (unsigned int q_point = 0; q_point < n_face_q_points;

++q_point) { const double neumann_value = (exact_solution.gradient( fe_face_values.quadrature_point(q_point)) * fe_face_values.normal_vector(q_point)); for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_rhs(i) += (neumann_value * // g(x_q) fe_face_values.shape_value(i, q_point) * // phi_i(x_q) fe_face_values.JxW(q_point)); // dx }

        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)

              system_matrix.add(local_dof_indices[i],

                                local_dof_indices[j],

                                cell_matrix(i, j));

            system_rhs(local_dof_indices[i]) += cell_rhs(i);

          }

}

}

}

std::map<types::global_dof_index, double> boundary_values;

VectorTools::interpolate_boundary_values(dof_handler,

0,

Functions::ZeroFunction<3>(),

boundary_values);

MatrixTools::apply_boundary_values(boundary_values,

system_matrix,

solution,

system_rhs);

[Simon Sticko]

Hi,
I'm not sure I understand your question correctly, but if you want to impose an inhomogeneous Neumann boundary condition you would typically write a small class representing the boundary condition you want:

class NeumannBoundaryCondition : public Function<3>
{
public:
  double
  value(const Point<3> &point, const unsigned int component = 0) const
  {
    // Implement the boundary condition you want here.
    return 1;
  }
};

Then just use this when you assemble:

const NeumannBoundaryCondition boundary_condition;

const Point<3> &point = fe_face_values.quadrature_point(q_point);

const double neumann_value = boundary_condition.value(point);

cell_rhs(i) += (neumann_value *
                fe_face_values.shape_value(i, q_point) *
                fe_face_values.JxW(q_point));


Does that answer your question?

Best,
Simon

[Samuel Rodriguez]

Hi Simon,

It seems I might not have understood my own question. You answered the question I was trying to ask. Would you be able to help me out with one more thing? I followed your advice but the boundary conditions do not seem bedo not seem to be getting the correct boundary conditions, so I don't think they are being applied correctly. Here is the simple piece of code that I wrote:

class NeumannBoundaryCondition : public Function<3>

{

public:

  double

  value(const Point<3> &point,  const unsigned int component = 0) const

  {

    return component;

  }

};

This neumann boundary condition is only meant to apply a constant value of zero on the sides. Here is how I implemented it:

QGauss<3> quadrature_formula(fe.degree + 1);
QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);

FEFaceValues<3> fe_face_values(fe,
       face_quadrature_formula,
       update_values | update_quadrature_points |
       update_normal_vectors |

       update_JxW_values);const NeumannBoundaryCondition boundary_condition;

const unsigned int n_face_q_points = face_quadrature_formula.size();

for(unsigned int face_number = 0; face_number < GeometryInfo<3>::faces_per_cell; ++face_number)
{
  if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->boundary_id() == 0))
  {
   fe_face_values.reinit(cell, face_number);

   for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {

     const double neumannvalue = boundary_condition.value(fe_face_values.quadrature_point(q_point));
     for (unsigned int i = 0; i < dofs_per_cell; ++i) {
       const unsigned int component_i = fe.system_to_component_index(i).first;
       cell_rhs(i) += (neumannvalue *
       fe_face_values.shape_value(i, q_point) *
       fe_face_values.JxW(q_point));
       }
     }
 
  }
}

This seems correct to me, but the code is showing a cyldiner that looks like this:

which doesn't appear to be the right solution. Even when I change the value of component to something ridiculous such as 10,000,000 the solution on the cylinder does not change. Would you be able to tell me what I am doing wrong in the implementation?

- Samuel

[Simon Sticko]

Hi,
From the figure you can't really tell whether the boundary condition is fulfilled or not. If you want to check if the normal derivative is reasonable (almost equal to 0) at the boundary, the easiest is probably to plot the solution on a line in the radial direction. This can be done by using
Filters -> Data Analysis -> Plot Over Line in Paraview.

I don't quite understand what you mean when you say that you have tried changing the value of component. Component is intended to be used for vector valued problems. If you are solving Poisson, it doesn't make sense to set it to something different than 0.

Note that, right now you are calling the function NeumannBoundaryCondition::value with a single argument (of type Point<3>). This means that you get the default argument for component: component=0.

Best,
Simon

[Samuel Rodriguez]

Hi Simon,

I am a slow learner but you were correct. Here is the first image of the function plotted along the radius. As you can see there is only noise but value of the solution is constant. I can also change the value of the solution by changing the value of the boundary condition as you can also see below:

I was able to change my boundary condition to something ridiculous like:

class NeumannBoundaryCondition : public Function<3>

{

public:

  double

  value(const Point<3> &point,  const unsigned int component = 0) const

  {

    double test = 1;

    for(int i = 0; i < 3; i++)

      test += 20 + std::pow(point(i), 20) + std::pow(point(i), 10) + exp(point(i));

    std::cout << test << std::endl;

    return test;

  }

};

The reason this was not initially working was because when I implemented the boundary condition below:

      for(unsigned int face_number = 0; face_number < GeometryInfo<3>::faces_per_cell; ++face_number)

      {

       if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->boundary_id() == 0))

        {

         fe_face_values.reinit(cell, face_number);

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

          {

           double neumannvalue = boundary_condition.value(fe_face_values.quadrature_point(q_point));

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

            {

             const unsigned int component_i = fe.system_to_component_index(i).first;

             cell_rhs(i) += (neumannvalue *

                             fe_face_values.shape_value(i, q_point) *

                             fe_face_values.JxW(q_point));

            }

          }

        

       }

     }

I never added it the right hand side (I think that I am using that terminology correctly). I added the previous lines of code right above these

      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)

          system_matrix.add(local_dof_indices[i],

                            local_dof_indices[j],

                            cell_matrix(i, j));

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

        system_rhs(local_dof_indices[i]) += cell_rhs(i);

and the code seems to be working fine! Is there anything that I am missing in your opinion?

Best,

Samuel Rodriguez

[Simon Sticko]

Hi,

> and the code seems to be working fine! Is there anything that I am

> missing in your opinion?

No. If you're happy, I'm happy.

Best,

Simon