Hi everyone,
I'm trying to solve a nonlinear ODE with Dirichlet conditions using
Newton's method. For this I give an initial guess for the solution ( y =
(1/2) x ) and then update each Newton step with an update that ought to be
zero at the boundaries (hence preserving the boundary conditions of the
initial guess). However, for some reason my update is taking on nonzero
values at the boundaries.
To enforce these boundary conditions, I initialize constraints with:
```
constraints.clear();
dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dealii::VectorTools::
interpolate_boundary_values(dof_handler,
0,
dealii::Functions::ZeroFunction<dim>(1),
constraints);
constraints.close();
```
I use the `distribute_local_to_global` in the assembly function, and then
once I've solved the linear system (directly, with UMFPACK), I use
`constraints.distribute`. Not sure what else I need to do to enforce those
boundary conditions.
I'm attaching the code (<300 lines including white space) in case that is
helpful. Any help is appreciated!
- Lucas
--
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/72a84737-055f-4fbc-b7b7-ad4400f687ebn%40googlegroups.com.
#ifndef DZYALOSHINSKII_SYSTEM_HPP
#define DZYALOSHINSKII_SYSTEM_HPP
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <string>
class DzyaloshinskiiSystem
{
public:
DzyaloshinskiiSystem(double eps_, unsigned int degree);
void make_grid(unsigned int n_refines);
void setup_system();
void assemble_system();
void solve_and_update(double newton_step);
bool run_newton_method(double tol,
unsigned int max_iter,
double newton_step = 1.0);
void output_solution(std::string filename);
private:
static constexpr int dim = 1;
dealii::Triangulation<dim> tria;
dealii::FE_Q<dim> fe;
dealii::DoFHandler<dim> dof_handler;
dealii::AffineConstraints<double> constraints;
dealii::SparsityPattern sparsity_pattern;
dealii::SparseMatrix<double> system_matrix;
dealii::Vector<double> solution;
dealii::Vector<double> system_rhs;
double eps;
};
#endif
#include "LiquidCrystalSystems/DzyaloshinskiiSystem.hpp"
#include <string>
#include <ostream>
int main(int ac, char* av[])
{
double eps = 0.5;
unsigned int degree = 1;
unsigned int n_refines = 8;
double tol = 1e-8;
unsigned int max_iter = 1;
std::string initial_filename("initial_dzyaloshinskii_solution.vtu");
std::string filename("dzyaloshinskii_solution_");
filename += std::to_string(eps);
filename += std::string(".vtu");
DzyaloshinskiiSystem dzyaloshinskii_system(eps, degree);
dzyaloshinskii_system.make_grid(n_refines);
dzyaloshinskii_system.setup_system();
dzyaloshinskii_system.output_solution(initial_filename);
bool converged = dzyaloshinskii_system.run_newton_method(tol, max_iter);
dzyaloshinskii_system.output_solution(filename);
std::cout << converged << "\n";
return 0;
}
#include "DzyaloshinskiiSystem.hpp"
#include <deal.II/base/types.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/point.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/base/table.h>
#include <deal.II/numerics/vector_tools.h>
#include <cmath>
#include <fstream>
#include <limits>
DzyaloshinskiiSystem::DzyaloshinskiiSystem(double eps_, unsigned int degree)
: fe(degree)
, dof_handler(tria)
, eps(eps_)
{}
void DzyaloshinskiiSystem::make_grid(unsigned int n_refines)
{
double left_endpoint = 0.0;
double right_endpoint = M_PI;
dealii::GridGenerator::hyper_cube(tria, left_endpoint, right_endpoint);
tria.refine_global(n_refines);
}
void DzyaloshinskiiSystem::setup_system()
{
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
{
dealii::Table<2, double> polynomial_exponents(2, dim);
std::vector<double> polynomial_coefficients(2);
polynomial_exponents[0][0] = 1.0;
polynomial_exponents[1][0] = 1.0;
polynomial_coefficients[0] = 0;
polynomial_coefficients[1] = 0.5;
dealii::Functions::Polynomial<dim>
initial_configuration(polynomial_exponents,
polynomial_coefficients);
constraints.clear();
dealii::DoFTools::make_hanging_node_constraints(dof_handler,
constraints);
dealii::VectorTools::
interpolate_boundary_values(dof_handler,
0,
initial_configuration,
constraints);
constraints.close();
dealii::VectorTools::project(dof_handler,
constraints,
dealii::QGauss<dim>(fe.degree + 1),
initial_configuration,
solution);
}
constraints.clear();
dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);
dealii::VectorTools::
interpolate_boundary_values(dof_handler,
0,
dealii::Functions::ZeroFunction<dim>(1),
constraints);
constraints.close();
dealii::DynamicSparsityPattern dsp(dof_handler.n_dofs());
dealii::DoFTools::
make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
void DzyaloshinskiiSystem::assemble_system()
{
dealii::QGauss<dim> quadrature_formula(fe.degree + 1);
dealii::FEValues<dim> fe_values(fe,
quadrature_formula,
dealii::update_values |
dealii::update_gradients |
dealii::update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
dealii::FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
dealii::Vector<double> cell_rhs(dofs_per_cell);
std::vector<dealii::types::global_dof_index>
local_dof_indices(dofs_per_cell);
std::vector<double> phi(n_q_points);
std::vector<dealii::Tensor<1, dim>> dphi(n_q_points);
std::vector<dealii::Point<dim>> quad_points(n_q_points);
double theta;
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
fe_values.get_function_values(solution, phi);
fe_values.get_function_gradients(solution, dphi);
quad_points = fe_values.get_quadrature_points();
for (const auto q : fe_values.quadrature_point_indices())
{
theta = quad_points[q](0);
for (const auto i : fe_values.dof_indices())
for (const auto j : fe_values.dof_indices())
cell_matrix(i, j) -=
(
( (1 - eps * std::cos(2 * (phi[q] - theta)))
* fe_values.shape_grad(i, q)
* fe_values.shape_grad(j, q) )
+
( dphi[q][0] * 2 * eps * std::sin(2 * (phi[q] - theta))
* fe_values.shape_grad(i, q)[0]
* fe_values.shape_value(j, q) )
+
( (4 * eps * std::cos(2 * (phi[q] - theta))
* dphi[q][0] * (dphi[q][0] - 1)
+
(2 * dphi[q][0] - dphi[q][0] * dphi[q][0])
* 2 * eps * std::cos(2 * (phi[q] - theta)) )
* fe_values.shape_value(i, q)
* fe_values.shape_value(j, q) )
)
* fe_values.JxW(q);
for (const auto i : fe_values.dof_indices())
cell_rhs(i) +=
(
( fe_values.shape_grad(i, q)[0]
* (1 - eps * std::cos(2 * (phi[q] - theta)))
* dphi[q][0] )
+
( fe_values.shape_value(i, q)
* dphi[q][0] * dphi[q][0]
* eps * std::sin(2 * (phi[q] - theta)) )
)
* fe_values.JxW(q);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
void DzyaloshinskiiSystem::solve_and_update(double newton_step)
{
dealii::SparseDirectUMFPACK solver;
solver.solve(system_matrix, system_rhs);
constraints.distribute(system_rhs);
dealii::DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(system_rhs, "dzyaloshinskii_solution");
data_out.build_patches();
std::ofstream output(std::string("dzyaloshinskii_update.vtu"));
data_out.write_vtu(output);
solution.add(newton_step, system_rhs);
}
bool DzyaloshinskiiSystem::
run_newton_method(double tol, unsigned int max_iters, double newton_step)
{
double res = std::numeric_limits<double>::max();
unsigned int iters = 0;
while (res >= tol && iters < max_iters)
{
assemble_system();
res = system_rhs.l2_norm();
solve_and_update(newton_step);
++iters;
}
return (res < tol);
}
void DzyaloshinskiiSystem::output_solution(std::string filename)
{
dealii::DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "dzyaloshinskii_solution");
data_out.build_patches();
std::ofstream output(filename);
data_out.write_vtu(output);
}
