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 
> <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 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);
>
>
>
>

-- 
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 the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/2d770452-ec8c-4c8f-bffc-9f385f21009bo%40googlegroups.com.

Reply via email to