# [deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

```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

On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>
> 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 - 1> face_quadrature_formula(fe.degree + 1);
> FEFaceValues<3> fe_face_values(fe,
>        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.
>      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:
>
>
> [image: Screen Shot 2020-06-29 at 5.50.19 PM.png]
> 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
>
>
>
> On Friday, June 26, 2020 at 2:21:19 AM UTC-7, Simon Sticko wrote:
>>
>> 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));
>>
>>
>>
>> Best,
>> Simon
>>
>> On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> 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
>>> <https://git.ligo.org/geoffrey-lovelace/NumericalCoatingThermalNoise>
>>>
>>> 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
>>>
>>> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
>>>
>>> FEValues<3> fe_values(fe,
>>> update_values | update_gradients | update_JxW_values);
>>> FEFaceValues<3> fe_face_values(fe,
>>> 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.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)
>>> 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.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)
>>>                                 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);
>>>
>>>
>>>
>>>

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see