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/fe8fc430-ae2b-405e-a7ea-29a088ddd6c5o%40googlegroups.com.
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/grid/grid_out.h>
#include <fstream>
#include <iostream>
using namespace dealii; // Need this
// My Class
class ASTLEC
{
public:
ASTLEC();
//~ASTLEC();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<3> triangulation;
FE_Q<3> fe;
DoFHandler<3> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
ASTLEC::ASTLEC() : fe(3), dof_handler(triangulation) {}
void ASTLEC::make_grid(){
GridGenerator::cylinder(triangulation, 1, 1);
triangulation.refine_global(2);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
// As the name implies, here is where we set up all objects and data structures that set up the problem as needed.
void ASTLEC::setup_system(){
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
void ASTLEC::assemble_system(){
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);
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() == 1))
{
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
// top face
}
}
}
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);
}
void ASTLEC::solve(){
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void ASTLEC::output_results() const {
DataOut<3> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
void ASTLEC::run() {
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
int main() {
deallog.depth_console(2);
ASTLEC objectofmine;
objectofmine.run();
std::cout << "Worked again.\n";
return 0;
}
