Hello all,
After I studied and understood the tutorials of step-8 and step-49. I
tried to use step-49 to import external grids to to calculate the force
condition and stress analysis,
So I modified the statement of step8 grid generation,
while failed...
Did I make a mistake or are there other ways to finish the force analysis
from external grids?
Attach my code...
Thanks!
--
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/39d0cd44-2231-4c10-a996-62de35a87c39n%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.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/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.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/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.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 <fstream>
#include <iostream>
#include <map>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,
std::vector<Tensor<1, dim>> & values)
{
Assert(values.size() == points.size(),
ExcDimensionMismatch(values.size(), points.size()));
Assert(dim >= 2, ExcNotImplemented());
Point<dim> point_1;
point_1(0) = 1.0;
point_1(1) = 1.0;
for (unsigned int point_n = 0; point_n < points.size();
++point_n)
{
if (((points[point_n] - point_1).norm_square() < 0.2 *
0.2))
values[point_n][1] = -1.0;
else
values[point_n][1] = 0.0;
}
}
template <int dim>
ElasticProblem<dim>::ElasticProblem()
: dof_handler(triangulation)
, fe(FE_Q<dim>(1), dim)
{}
template <int dim>
void print_mesh_info(const Triangulation<dim>& triangulation,
const std::string& filename)
{
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells()
<< std::endl;
{
std::map<types::boundary_id, unsigned int>
boundary_count;
for (const auto& face :
triangulation.active_face_iterators())
if (face->at_boundary())
boundary_count[face->boundary_id()]++;
std::cout << " boundary indicators: ";
for (const std::pair<const types::boundary_id, unsigned
int>& pair :
boundary_count)
{
std::cout << pair.first << "(" << pair.second
<< " times) ";
}
std::cout << std::endl;
}
void grid_1()
{
Triangulation<2> triangulation;
GridIn<2> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("example.msh");
gridin.read_msh(f);
print_mesh_info(triangulation, "grid-1.vtu");
}
template <int dim>
void ElasticProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
constraints);
for (auto &face : triangulation.active_face_iterators())
if (std::fabs(face->center()(1) - (-1.0)) < 1e-12)
face->set_boundary_id(1);
VectorTools::interpolate_boundary_values(dof_handler,
1,
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);
system_matrix.reinit(sparsity_pattern);
}
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = 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);
for (const auto &cell : dof_handler.active_cell_iterators())
{
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 (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);
}
}
template <int dim>
void ElasticProblem<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
constraints.distribute(solution);
}
template <int dim>
void ElasticProblem<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.emplace_back("displacement");
break;
case 2:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
break;
case 3:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
solution_names.emplace_back("z_displacement");
break;
default:
Assert(false, ExcNotImplemented());
}
data_out.add_data_vector(solution, solution_names);
data_out.build_patches();
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
template <int dim>
void ElasticProblem<dim>::run()
{
make_grid();
setup_system();
std::cout << " Number of degrees of freedom: " <<
dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve();
output_results();
}
} // namespace Step8
}
int main()
{
try
{
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;
}
