Hello,
I am trying to solve a 2D elasticity problem (step 8) by considering a
cantilever beam which is fixed at the right end and on the left end there
is traction.
But I am facing the following issues:
1. While considering only body forces my code is not working (no traction
case).
2. When I am neglecting the body forces and considering only traction then
code is working but solutions are not correct (no body force case).
Please see the attached code for both issues and advise.
--
Thanks and Regards
Deepika Kushwah
PhD Scholar, School of Mechanical Science,
IIT Goa
--
**************************************************************************
This e-mail is for the sole use of the intended recipient(s) and may
contain confidential and privileged information. If you are not the
intended recipient, please contact the sender by reply e-mail and destroy
all copies and the original message. Any unauthorized review, use,
disclosure, dissemination, forwarding, printing or copying of this email
is strictly prohibited and appropriate legal action will be taken.
************************************************************************************************
--
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/CAKTdrppatXH%3Dio%3DpETZbwsJm-%2BKT856fjR0KPckQ9dEemfD2kQ%40mail.gmail.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/tensor.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/solver_cg.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
//#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <fstream>
#include <iostream>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem ();
~ElasticProblem ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
//ConstraintMatrix hanging_node_constraints;
//QGauss<dim> quadrature_formula; //Quadrature
//QGauss<dim-1> face_quadrature_formula; //Face Quadrature
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
std::vector<std::string> nodal_solution_names;
int ncx, ncy;
};
template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,
std::vector<Tensor<1, dim>> & values)
{
AssertDimension(values.size(), points.size());
Assert(dim >= 2, ExcNotImplemented());
for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
values[point_n][0] = 0;
values[point_n][1] = 10;
}
}
/*
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
8*/
template <int dim>
ElasticProblem<dim>::ElasticProblem ()
:
dof_handler (triangulation),
fe (FE_Q<dim>(1), dim)
// quadrature_formula(quadRule),
// face_quadrature_formula(quadRule)
{}
/*
template <int dim>
ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
}
*/
template <int dim>
void ElasticProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
constraints.clear();
//hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
/*
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(dim),
constraints);
*/
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
//hanging_node_constraints.close ();
/*
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
constraints.condense (sparsity_pattern);
sparsity_pattern.compress();
*/
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
void ElasticProblem<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
//QGauss<dim-1> face_quadrature_formula(2);
// RightHandSide<dim> right_hand_side;
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
/*
FEFaceValues<dim> 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.n_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);
std::vector<double> lambda_values (n_q_points);
std::vector<double> mu_values (n_q_points);
// const RightHandSide<dim> right_hand_side;
//std::vector<double> rhs_values(n_q_points);
Functions::ConstantFunction<dim> lambda(1.), mu(1.);
std::vector<Tensor<1, dim>> rhs_values(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
//double valE=1.0;
//double valnu=0.3;
//double vallambda=(valE*valnu)/((1+valnu)*(1-2*valnu));
//double vallambda=(valE*valnu)/((1+valnu)*(1-valnu));
//double valmu=valE/(2*(1+valnu));
// ConstantFunction<dim> lambda(vallambda), mu(valmu);
// RightHandSide<dim> right_hand_side;
// std::vector<Vector<double> > rhs_values (n_q_points,
// Vector<double>(dim));
std::map<types::global_dof_index,double> boundary_values;
/*
typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active(),
endc = dof_handler.end();
int counter=0;
for (; cell!=endc; ++cell)
{
counter++;
cell_matrix = 0;
cell_rhs = 0;
//double x = 1;
fe_values.reinit (cell);
*/
lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
mu.value_list (fe_values.get_quadrature_points(), mu_values);
right_hand_side(fe_values.get_quadrature_points(), rhs_values);
/*
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int
component_i = fe.system_to_component_index(i).first;
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const unsigned int
component_j = fe.system_to_component_index(j).first;
for (unsigned int q_point=0; q_point<n_q_points;
++q_point)
{
cell_matrix(i,j)
+=
(
(fe_values.shape_grad(i,q_point)[component_i] *
fe_values.shape_grad(j,q_point)[component_j] *
lambda_values[q_point])
+
(fe_values.shape_grad(i,q_point)[component_j] *
fe_values.shape_grad(j,q_point)[component_i] *
mu_values[q_point])
+
((component_i == component_j) ?
(fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point) *
mu_values[q_point]) :
0)
)
*
fe_values.JxW(q_point);
}
}
}
*/
for (const unsigned int i : fe_values.dof_indices())
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (const unsigned int j : fe_values.dof_indices())
{
const unsigned int component_j =
fe.system_to_component_index(j).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
{
cell_matrix(i, j) +=
(
(fe_values.shape_grad(i, q_point)[component_i] *
fe_values.shape_grad(j, q_point)[component_j] *
lambda_values[q_point])
+
(fe_values.shape_grad(i, q_point)[component_j] *
fe_values.shape_grad(j, q_point)[component_i] *
mu_values[q_point])
+
((component_i == component_j) ?
(fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point) *
mu_values[q_point]) :
0)
) *
fe_values.JxW(q_point);
}
}
}
for (const unsigned int i : fe_values.dof_indices())
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (const unsigned int q_point :
fe_values.quadrature_point_indices())
cell_rhs(i) += fe_values.shape_value(i, q_point) *
rhs_values[q_point][component_i] *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
/*
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
cell_rhs(i)=0;
const unsigned int
component_i = fe.system_to_component_index(i).first;
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += fe_values.shape_value(i,q_point) *
rhs_values[q_point](component_i) *
fe_values.JxW(q_point);
}
*/
/*
for (const auto &face : cell->face_iterators())
if (face->at_boundary() && (face->boundary_id() == 1))
{
fe_face_values.reinit(cell, face);
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));
const auto &x_q = fe_values.quadrature_point(n_face_q_points );
for (unsigned int i = 0; i < dofs_per_cell; ++i)
// double x = fe_face_values.quadrature_point(q)[0];
cell_rhs(i) +=
(fe_face_values.shape_value(i, q_point) * // phi_i(x_q)
//right_hand_side.value(x_q) * // f(x_q)
10.0e4* // g(x_q)
fe_face_values.JxW(q_point)); // dx
}
}
*/
// assemble cell_rhs, which is 0 except at cell with (0,ncy), there it
is -1 at the corresponding vertex (not using the RightHandSide-Class)
/*
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);
}
*/
//rhs is 0 except at cell with (0,ncy), there it is -1 at the
corresponding vertex (not using the RightHandSide-Class)
for (unsigned int f=0; f < GeometryInfo<dim>::faces_per_cell; ++f){
// find faces at right side of rectangle and set Dirichlet
boundary in y-Direction
if ((cell->face(f)->center()[0] == 15)){
for (unsigned int v=0; v <
GeometryInfo<dim>::vertices_per_face; ++v) {
boundary_values[cell->face(f)->vertex_dof_index(v, 1)] = 0;
}
}
}
for (unsigned int f=0; f < GeometryInfo<dim>::faces_per_cell; ++f){
// find faces at right side of rectangle and set Dirichlet
boundary in x-Direction
if ((cell->face(f)->center()[0] == 15)){
for (unsigned int v=0; v <
GeometryInfo<dim>::vertices_per_face; ++v) {
boundary_values[cell->face(f)->vertex_dof_index(v, 0)] = 0;
}
}
}
// A different possibility would be to set the boundary_indicator and
fill the boundary_values list by applying the function
// "interpolate_boundary_values"
(cell->face(f)->set_boundary_indicator(1))
}
constraints.condense (system_matrix);
constraints.condense (system_rhs);
std::map<types::global_dof_index,double> boundary_values;
/*
VectorTools::interpolate_boundary_values (dof_handler,
0,
ConstantFunction<dim>(2),
boundary_values);
*/
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim>
void ElasticProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
constraints.distribute (solution);
}
template <int dim>
void ElasticProblem<dim>::refine_grid ()
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(dof_handler,
QGauss<dim - 1>(fe.degree + 1),
{},
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
0.3,
0.03);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void ElasticProblem<dim>::output_results (const unsigned int cycle) const
{
std::string filename = "solution-";
filename += ('0' + cycle);
Assert (cycle < 1, ExcInternalError());
filename += ".vtk";
std::ofstream output (filename.c_str());
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.push_back ("displacement");
break;
case 2:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
break;
case 3:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
solution_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}
data_out.add_data_vector (solution, solution_names);
data_out.build_patches ();
//std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
data_out.write_vtk (output);
}
/*
std::vector<std::string> solution_names (dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names,
DataOut<dim>::type_dof_data, data_component_interpretation);
data_out.build_patches ();
data_out.write_vtk (output);
}
*/
template <int dim>
void ElasticProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<1; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
// We generate a rectangular grid with 60 cells in the x-direction
(ncx) and 40 cells in the y-direction (ncy)
ncx=15, ncy=5;
// The grid generator needs the lower left corner and the upper
right corner of the rectangle
const Point<2> point_1 (0,0);
const Point<2> point_2 (ncx,ncy);
std::vector<unsigned int> subdivisions;
subdivisions.push_back (ncx);
subdivisions.push_back (ncy);
GridGenerator::subdivided_hyper_rectangle (triangulation,
subdivisions, point_1, point_2);
triangulation.refine_global(1);
for (const auto &cell : triangulation.cell_iterators())
for (const auto &face : cell->face_iterators())
{
const auto center = face->center();
if ((std::fabs(center(0) - (0)) < 1e-12) )//||
// (std::fabs(center(1) - (0)) < 1e-12))
face->set_boundary_id(1);
}
}
else
refine_grid ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
setup_system ();
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
output_results (cycle);
}
}
}
int main ()
{
try
{
dealii::deallog.depth_console (0);
Step8::ElasticProblem<2> elastic_problem_2d;
elastic_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
/* ---------------------------------------------------------------------
* This test case is based on the dealii tutorial step-8. It calculates a
* Cantilever beam poblem
* BCs are : right end is fixed
: left end is subjected to force or traction
* ---------------------------------------------------------------------
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/tensor.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/solver_cg.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
//#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <fstream>
#include <iostream>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem ();
~ElasticProblem ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
//ConstraintMatrix hanging_node_constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
int ncx, ncy;
};
/*
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override;
};
*/
template <int dim>
ElasticProblem<dim>::ElasticProblem ()
:
dof_handler (triangulation),
fe (FE_Q<dim>(1), dim)
{}
template <int dim>
ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
}
template <int dim>
void ElasticProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
constraints.clear();
//hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
constraints.close();
//hanging_node_constraints.close ();
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
constraints.condense (sparsity_pattern);
sparsity_pattern.compress();
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
void ElasticProblem<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
QGauss<dim-1> face_quadrature_formula(2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> 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);
std::vector<double> lambda_values (n_q_points);
std::vector<double> mu_values (n_q_points);
Functions::ConstantFunction<dim> lambda(1.), mu(1.);
std::vector<Tensor<1, dim>> rhs_values(n_q_points);
//RightHandSide<dim> right_hand_side;
std::map<types::global_dof_index,double> boundary_values;
typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active(),
endc = dof_handler.end();
int counter=0;
for (; cell!=endc; ++cell)
{
counter++;
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
mu.value_list (fe_values.get_quadrature_points(), mu_values);
// right_hand_side(fe_values.get_quadrature_points(), rhs_values);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int
component_i = fe.system_to_component_index(i).first;
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const unsigned int
component_j = fe.system_to_component_index(j).first;
for (unsigned int q_point=0; q_point<n_q_points;
++q_point)
{
cell_matrix(i,j)
+=
(
(fe_values.shape_grad(i,q_point)[component_i] *
fe_values.shape_grad(j,q_point)[component_j] *
lambda_values[q_point])
+
(fe_values.shape_grad(i,q_point)[component_j] *
fe_values.shape_grad(j,q_point)[component_i] *
mu_values[q_point])
+
((component_i == component_j) ?
(fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point) *
mu_values[q_point]) :
0)
)
*
fe_values.JxW(q_point);
}
}
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
cell_rhs(i)=0;
const unsigned int
component_i = fe.system_to_component_index(i).first;
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += fe_values.shape_value(i,q_point) *0*
// rhs_values[q_point](component_i) *
fe_values.JxW(q_point);
}
// for (const auto &face : cell->face_iterators())
// if (face->at_boundary() && (face->boundary_id() == 1))
// {
//fe_face_values.reinit(cell, face);
// for (unsigned int q_point = 0; q_point < n_face_q_points;
// ++q_point)
//{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary() == true
&& cell->face(f)->boundary_id() == 1)
{
fe_face_values.reinit (cell, f);
for (unsigned int fq_point=0; fq_point<n_face_q_points;
++fq_point)
{
const Tensor<1, dim> &N =
fe_face_values.normal_vector(fq_point);
static const double T = 50e3; // arbitraty magnitude of
traction
const Tensor<1,dim> traction = (T * N);
for (int i = 0; i<dofs_per_cell; ++i)
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
cell_rhs(i) += traction[component_i]
* fe_face_values.shape_value(i,fq_point)
* fe_face_values.JxW(fq_point);
}
}
}
//const auto &x_q = fe_values.quadrature_point(n_face_q_points );
/*
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(fe_face_values.shape_value(i, q_point) * // phi_i(x_q)
//right_hand_side.value(x_q) * // f(x_q)
10.e3*
// g(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);
}
// displacement boundary condition for right end
for (unsigned int f=0; f < GeometryInfo<dim>::faces_per_cell; ++f){
// find faces at right side of rectangle and set Dirichlet
boundary in x-Direction
if ((cell->face(f)->center()[0] == 15)){
for (unsigned int v=0; v <
GeometryInfo<dim>::vertices_per_face; ++v) {
boundary_values[cell->face(f)->vertex_dof_index(v, 0)] = 0;
}
}
}
for (unsigned int f=0; f < GeometryInfo<dim>::faces_per_cell; ++f){
// find faces at right side of rectangle and set Dirichlet
boundary in y-Direction
if ((cell->face(f)->center()[0] == 15)){
for (unsigned int v=0; v <
GeometryInfo<dim>::vertices_per_face; ++v) {
boundary_values[cell->face(f)->vertex_dof_index(v, 1)] = 0;
}
}
}
// A different possibility would be to set the boundary_indicator and
fill the boundary_values list by applying the function
// "interpolate_boundary_values"
(cell->face(f)->set_boundary_indicator(1))
}
constraints.condense (system_matrix);
constraints.condense (system_rhs);
//std::map<types::global_dof_index,double> boundary_values;
/*
VectorTools::interpolate_boundary_values (dof_handler,
0,
ConstantFunction<dim>(2),
boundary_values);
*/
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim>
void ElasticProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
constraints.distribute (solution);
}
template <int dim>
void ElasticProblem<dim>::refine_grid ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
//typename FunctionMap<dim>::type neumann_boundary;
/* KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(2),
//neumann_boundary,
solution,
estimated_error_per_cell);
*/
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.03);
triangulation.execute_coarsening_and_refinement ();
}
template <int dim>
void ElasticProblem<dim>::output_results (const unsigned int cycle) const
{
std::string filename = "solution-";
filename += ('0' + cycle);
Assert (cycle < 1, ExcInternalError());
filename += ".vtk";
std::ofstream output (filename.c_str());
/*DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.push_back ("displacement");
break;
case 2:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
break;
case 3:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
solution_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}
data_out.add_data_vector (solution, solution_names);
data_out.build_patches ();
data_out.write_vtk (output);*/
std::vector<std::string> solution_names (dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names,
DataOut<dim>::type_dof_data, data_component_interpretation);
data_out.build_patches ();
data_out.write_vtk (output);
}
template <int dim>
void ElasticProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<1; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
// We generate a rectangular grid with 15 cells in the x-direction
(ncx) and 5 cells in the y-direction (ncy)
ncx=15, ncy=5;
// The grid generator needs the lower left corner and the upper
right corner of the rectangle
const Point<2> point_1 (0,0);
const Point<2> point_2 (ncx,ncy);
std::vector<unsigned int> subdivisions;
subdivisions.push_back (ncx);
subdivisions.push_back (ncy);
GridGenerator::subdivided_hyper_rectangle (triangulation,
subdivisions, point_1, point_2,true);
for (const auto &cell : triangulation.cell_iterators())
for (const auto &face : cell->face_iterators())
{
const auto center = face->center();
if ((std::fabs(center(0) - (0)) < 1e-12) ||
(std::fabs(center(1) - (0)) < 1e-12))
face->set_boundary_id(1);
}
}
else
refine_grid ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
setup_system ();
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
output_results (cycle);
}
}
}
int main ()
{
try
{
dealii::deallog.depth_console (0);
Step8::ElasticProblem<2> elastic_problem_2d;
elastic_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
