[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder
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
On Wednesday, July 1, 2020 at 7:51:11 AM UTC+2, Samuel Rodriguez wrote:
>
> 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:
>
> [image: Screen Shot 2020-06-30 at 2.20.11 PM.png]
>
> [image: Screen Shot 2020-06-30 at 10.22.07 PM.png]
>
> 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
>
>
>
>
>
>
>
>
>
>
>
> On Tuesday, June 30, 2020 at 3:33:25 AM UTC-7, Simon Sticko wrote:
>>
>> 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> 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);c
[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> 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:
>
>
> [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));
>>
>>
>> Does that answer your question?
>>
>> 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
>>>
[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder
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:
[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));
>
>
> Does that answer your question?
>
> 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
>>
>>
>> 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 ne
[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder
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
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
>
>
> 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 cell_matrix(dofs_per_cell, dofs_per_cell);
> Vector cell_rhs(dofs_per_cell);
>
> std::vector 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
