#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/lac/affine_constraints.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 <fstream>
#include <iostream>
using namespace dealii;
//class Boundary_values : public Function<2>
//{
//public:
//virtual double value (const Point <2> &p,
//const unsigned int component = 0) const override;
//};
//double Boundary_values::value(const Point<2> & p, const unsigned int )
const
//{
//double s=0.0;
//s=p(0)-p(0)+900;
//return s;
//}
class Right_side : public Function<2>
{
public:
virtual double value (const Point <2> & p,
const unsigned int component = 0) const override;
};
double Right_side::value(const Point<2> & p, const unsigned int ) const
{
double return_value=0.0;
return_value = std::cos(2*dealii::numbers::PI*(8*p(0)+1*p(1)));
return return_value;
// Code here for logic of "l" of governing equation
}
class fea
{
public:
fea();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<2> mesh;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
//AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> K_matrix;
Vector<double> d_vector;
Vector<double> F_vector;
};
fea::fea()
: fe(1)
, dof_handler(mesh)
{}
void fea::make_grid()
{
const Point<2> p1(0,0);
const Point<2> p2(27,45);
GridGenerator::hyper_rectangle(mesh, p1, p2);
mesh.refine_global(3);
//code here
}
void fea::setup_system()
{
dof_handler.distribute_dofs(fe);
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
std::cout <<"Number of degrees of freedom:" <<dof_handler.n_boundary_dofs()
<< std::endl;
K_matrix.reinit(sparsity_pattern);
F_vector.reinit(dof_handler.n_dofs());
d_vector.reinit(dof_handler.n_dofs());
//code here
}
void fea::assemble_system()
{
QGauss<2> quadrature_formula(fe.degree+1);
Right_side right_hand_side;
FEValues<2> fe_values(fe,quadrature_formula,update_values |
update_gradients | update_JxW_values);
const unsigned int dofs_per_cell =fe.n_dofs_per_cell();
FullMatrix<double> local_stiffness_matrix(dofs_per_cell,dofs_per_cell);
Vector<double> local_rhs_vector(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);
local_stiffness_matrix=0;
local_rhs_vector=0;
for (const unsigned int q_index :fe_values.quadrature_point_indices())
for(const unsigned int a : fe_values.dof_indices())
{
for(const unsigned int b : fe_values.dof_indices())
local_stiffness_matrix(a,b) += (fe_values.shape_grad(a,
q_index)*fe_values.shape_grad(b, q_index)*fe_values.JxW(q_index));
const auto & x_q= fe_values.quadrature_point(q_index);
local_rhs_vector(a) += (fe_values.shape_value(a,
q_index)*right_hand_side.value(x_q)*fe_values.JxW(q_index));
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
{
for (const unsigned int j : fe_values.dof_indices())
K_matrix.add(local_dof_indices[i],local_dof_indices[j],local_stiffness_matrix(i,j));
F_vector(local_dof_indices[i])+=local_rhs_vector(i);
}
}
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler, 0,
ConstantFunction<2> (9000), boundary_values);
MatrixTools::apply_boundary_values(boundary_values, K_matrix, d_vector,
F_vector);
//code here
}
void fea::solve()
{
SolverControl solver_control(1000,1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(K_matrix,d_vector,F_vector,PreconditionIdentity());
//code here
}
void fea::output_results() const
{
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(d_vector, "solution");
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
//code here
}
void fea::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
Please help to rectify the error from the following code
int main()
{
deallog.depth_console(2);
fea laplace_problem;
laplace_problem.run();
return 0;
}
[image: image.png]
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/CAEgfrC5TQQW51Y2q3XKabO0z22RgeHaz5dOo%3D42gZ5SuROY1zg%40mail.gmail.com.
