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 >> <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/307218a3-3cd0-4ff7-8b67-928a20281e1fo%40googlegroups.com.