非线性有限元
#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/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.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 <cmath>
namespace Parameters
{
using namespace dealii;
struct allparameters
{
/*******************************************************************/
/*parameter function*/
allparameters(const std::string& input_file);
static void declare_parameters(ParameterHandler& prm);
void parse_parameters(ParameterHandler& prm);
/*******************************************************************/
/*Define FESystem*/
unsigned int poly_degree;
unsigned int quad_order;
/*******************************************************************/
/*Input elastic tensor*/
double C11;
double C12;
double C13;
double C14;
double C15;
double C16;
double C22;
double C23;
double C24;
double C25;
double C26;
double C33;
double C34;
double C35;
double C36;
double C44;
double C45;
double C46;
double C55;
double C56;
double C66;
/*******************************************************************/
/*solve parameters*/
double err;
unsigned int iteration_number;
/*******************************************************************/
/*Geometry parameters*/
double length_left;
double length_right;
unsigned int refine_number;
/*******************************************************************/
/*engstrain*/
double F11;
double F12;
double F13;
double F21;
double F22;
double F23;
double F31;
double F32;
double F33;
};
/****************************************************************************************************/
allparameters::allparameters(const std::string& input_file)
{
ParameterHandler prm;
declare_parameters(prm);
prm.parse_input(input_file);
parse_parameters(prm);
}
/*****************************************************************************************************/
void allparameters::declare_parameters(ParameterHandler& prm)
{
prm.enter_subsection("Finite strain PFM");
{
prm.declare_entry("Polynomial degree", "1", Patterns::Integer(0), "Displacement system polynomial order");
prm.declare_entry("Quadrature order", "3", Patterns::Integer(0), "Gauss quadrature order");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
prm.declare_entry("C11", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C12", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C13", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C14", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C15", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C16", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C22", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C23", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C24", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C25", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C26", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C33", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C34", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C35", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C36", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C44", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C45", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C46", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C55", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C56", "0", Patterns::Double(0.), "Elastic constant");
prm.declare_entry("C66", "0", Patterns::Double(0.), "Elastic constant");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
prm.declare_entry("err", "1e-12", Patterns::Double(), "solve error");
prm.declare_entry("iteration number", "2000", Patterns::Integer(), "iteration number");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
prm.declare_entry("length left", "0", Patterns::Double(), "Gemetry parameters");
prm.declare_entry("length right", "0", Patterns::Double(), "Gemetry parameters");
prm.declare_entry("refine number", "5", Patterns::Integer(), "Gemetry parameters");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
prm.declare_entry("F11", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F12", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F13", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F21", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F22", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F23", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F31", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F32", "0", Patterns::Double(), "enegstrain");
prm.declare_entry("F33", "0", Patterns::Double(), "enegstrain");
}
prm.leave_subsection();
}
void allparameters::parse_parameters(ParameterHandler& prm)
{
prm.enter_subsection("Finite strain PFM");
{
poly_degree = prm.get_integer("Polynomial degree");
quad_order = prm.get_integer("Quadrature order");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
C11 = prm.get_double("C11");
C12 = prm.get_double("C12");
C13 = prm.get_double("C13");
C14 = prm.get_double("C14");
C15 = prm.get_double("C15");
C16 = prm.get_double("C16");
C22 = prm.get_double("C22");
C23 = prm.get_double("C23");
C24 = prm.get_double("C24");
C25 = prm.get_double("C25");
C26 = prm.get_double("C26");
C33 = prm.get_double("C33");
C34 = prm.get_double("C34");
C35 = prm.get_double("C35");
C36 = prm.get_double("C36");
C44 = prm.get_double("C44");
C45 = prm.get_double("C45");
C46 = prm.get_double("C46");
C55 = prm.get_double("C55");
C56 = prm.get_double("C56");
C66 = prm.get_double("C66");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
err = prm.get_double("err");
iteration_number = prm.get_integer("iteration number");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
length_left = prm.get_double("length left");
length_right = prm.get_double("length right");
refine_number = prm.get_integer("refine number");
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
F11 = prm.get_double("F11");
F12 = prm.get_double("F12");
F13 = prm.get_double("F13");
F21 = prm.get_double("F21");
F22 = prm.get_double("F22");
F23 = prm.get_double("F23");
F31 = prm.get_double("F31");
F32 = prm.get_double("F32");
F33 = prm.get_double("F33");
}
prm.leave_subsection();
}
struct PointHistory
{
Vector<double> order_parameters;
Vector<double> stress_vector;
Vector<double> strain_vector;
PointHistory();
};
PointHistory::PointHistory()
{
order_parameters.reinit(12);
stress_vector.reinit(6);
strain_vector.reinit(6);
}
struct deformation
{
FullMatrix<double> Fi;
FullMatrix<double> Ei;
Vector<double> engstrain_vector;
FullMatrix<double> Fi_inv;
FullMatrix<double> Fi_inv_T;
FullMatrix<double> Fi_T;
FullMatrix<double> Fi_T_Fi;
FullMatrix<double> elastic_matrix;
FullMatrix<double> I;
double Ji;
deformation();
};
deformation::deformation()
{
Fi.reinit(3,3);
Ei.reinit(3, 3);
engstrain_vector.reinit(6);
Fi_inv.reinit(3, 3);
Fi_inv_T.reinit(3, 3);
Fi_T.reinit(3, 3);
Fi_T_Fi.reinit(3, 3);
elastic_matrix.reinit(6,6);
I.reinit(3, 3);
}
}
#include "Finite_strain_PFM.h"
namespace Finite_strain_PFM
{
using namespace dealii;
template <int dim>
class Finite_strain_problem
{
public:
Finite_strain_problem(const std::string& input_file);
void run();
private:
void create_grid();
void setup_system();
void parameter_setup();
void setup_quadrature_point_history();
void update_quadrature_point_history();
void assemble_small_strain_system();
void assemble_Finite_strain_system();
void solve_linear(int judge);//judge用来判断求解哪个线性方程组
void solve_nolinear();
void output_small_strain_displacement();
void output_Finite_strain_displacement(int number_of_iteration);
void output_small_strain_elastic_field();
void output_Finite_strain_elastic_field(int number_of_iteration);
std::vector<Parameters::PointHistory> quadrature_point_history;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> small_strain_stiff_matrix;
SparseMatrix<double> Finite_strain_tangent_stiff_matrix;
Parameters::allparameters parameters;
Vector<double> small_strain_system_solution;//小应变解
Vector<double> Finite_strain_system_solution;//每步结束后的解
Vector<double> Finite_strain_system_increment_solution;//增量解
Vector<double> small_strain_system_rhs;//小应变的右端项
Vector<double> Finite_strain_system_rhs;//使用Netwon_Raphson迭代所需的右端项
Parameters::deformation deformation_variable;
};
template <int dim>
Finite_strain_problem<dim>::Finite_strain_problem(const std::string& input_file)
: parameters(input_file)
, dof_handler(triangulation)
, fe(FE_Q<dim>(1), dim)
{}
template <int dim>
void Finite_strain_problem<dim>::parameter_setup()
{
clock_t start, end;
start=clock();
std::cout<<"start to set up parameters"<<std::endl;
deformation_variable.elastic_matrix(0, 0) = parameters.C11;
deformation_variable.elastic_matrix(0, 1) = parameters.C12;
deformation_variable.elastic_matrix(0, 2) = parameters.C13;
deformation_variable.elastic_matrix(0, 3) = parameters.C14;
deformation_variable.elastic_matrix(0, 4) = parameters.C15;
deformation_variable.elastic_matrix(0, 5) = parameters.C16;
deformation_variable.elastic_matrix(1, 0) = parameters.C12;
deformation_variable.elastic_matrix(1, 1) = parameters.C22;
deformation_variable.elastic_matrix(1, 2) = parameters.C23;
deformation_variable.elastic_matrix(1, 3) = parameters.C24;
deformation_variable.elastic_matrix(1, 4) = parameters.C25;
deformation_variable.elastic_matrix(1, 5) = parameters.C26;
deformation_variable.elastic_matrix(2, 0) = parameters.C13;
deformation_variable.elastic_matrix(2, 1) = parameters.C23;
deformation_variable.elastic_matrix(2, 2) = parameters.C33;
deformation_variable.elastic_matrix(2, 3) = parameters.C34;
deformation_variable.elastic_matrix(2, 4) = parameters.C35;
deformation_variable.elastic_matrix(2, 5) = parameters.C36;
deformation_variable.elastic_matrix(3, 0) = parameters.C14;
deformation_variable.elastic_matrix(3, 1) = parameters.C24;
deformation_variable.elastic_matrix(3, 2) = parameters.C34;
deformation_variable.elastic_matrix(3, 3) = parameters.C44;
deformation_variable.elastic_matrix(3, 4) = parameters.C45;
deformation_variable.elastic_matrix(3, 5) = parameters.C46;
deformation_variable.elastic_matrix(4, 0) = parameters.C15;
deformation_variable.elastic_matrix(4, 1) = parameters.C25;
deformation_variable.elastic_matrix(4, 2) = parameters.C35;
deformation_variable.elastic_matrix(4, 3) = parameters.C45;
deformation_variable.elastic_matrix(4, 4) = parameters.C55;
deformation_variable.elastic_matrix(4, 5) = parameters.C56;
deformation_variable.elastic_matrix(5, 0) = parameters.C16;
deformation_variable.elastic_matrix(5, 1) = parameters.C26;
deformation_variable.elastic_matrix(5, 2) = parameters.C36;
deformation_variable.elastic_matrix(5, 3) = parameters.C46;
deformation_variable.elastic_matrix(5, 4) = parameters.C56;
deformation_variable.elastic_matrix(5, 5) = parameters.C66;
/*for(int i=0;i<6;i++)
{
for (int j=0;j<6;j++)
std::cout<<deformation_variable.elastic_matrix(i,j)<<" ";
std::cout<<std::endl;
}*/
deformation_variable.Fi(0, 0) = parameters.F11;
deformation_variable.Fi(0, 1) = parameters.F12;
deformation_variable.Fi(0, 2) = parameters.F13;
deformation_variable.Fi(1, 0) = parameters.F21;
deformation_variable.Fi(1, 1) = parameters.F22;
deformation_variable.Fi(1, 2) = parameters.F23;
deformation_variable.Fi(2, 0) = parameters.F31;
deformation_variable.Fi(2, 1) = parameters.F32;
deformation_variable.Fi(2, 2) = parameters.F33;
/*for(int i=0;i<3;i++)
{
for (int j=0;j<3;j++)
std::cout<<deformation_variable.Fi(i,j)<<" ";
std::cout<<std::endl;
}*/
deformation_variable.I = 0;
deformation_variable.Fi_inv.invert(deformation_variable.Fi);
deformation_variable.Ji = deformation_variable.Fi.determinant();
for (int i = 0; i < 3; i++)
{
deformation_variable.I(i, i) = 1;
for (int j = 0; j < 3; j++)
{
deformation_variable.Fi_T(i, j) = deformation_variable.Fi(j, i);
deformation_variable.Fi_inv_T(i, j) = deformation_variable.Fi_inv(j, i);
}
}
deformation_variable.Fi_T.mmult(deformation_variable.Fi_T_Fi, deformation_variable.Fi);
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
deformation_variable.Ei(i, j) = deformation_variable.Fi_T_Fi(i, j) / 2 - deformation_variable.I(i, j) / 2;
deformation_variable.engstrain_vector(0) = deformation_variable.Ei(0, 0);
deformation_variable.engstrain_vector(1) = deformation_variable.Ei(1, 1);
deformation_variable.engstrain_vector(2) = deformation_variable.Ei(2, 2);
deformation_variable.engstrain_vector(3) = 2 * deformation_variable.Ei(1, 2);
deformation_variable.engstrain_vector(4) = 2 * deformation_variable.Ei(0, 2);
deformation_variable.engstrain_vector(5) = 2 * deformation_variable.Ei(0, 1);
end=clock();
std::cout<<"the time of set parameters "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
/***************/
}
template <int dim>
void Finite_strain_problem<dim>::create_grid()
{
clock_t start, end;
start=clock();
std::cout<<"start to creat grid"<<std::endl;
GridGenerator::hyper_cube(triangulation, parameters.length_left, parameters.length_right);
triangulation.refine_global(parameters.refine_number);
end=clock();
std::cout<<"the time of creat grid "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
//下面可增加施加边界条件的程序。
}
/******************************************设置系统，初始化参数*******************************/
template <int dim>
void Finite_strain_problem<dim>::setup_system()
{
clock_t start, end;
start=clock();
std::cout<<"start to set up system"<<std::endl;
dof_handler.distribute_dofs(fe);
small_strain_system_solution.reinit(dof_handler.n_dofs());
Finite_strain_system_solution.reinit(dof_handler.n_dofs());
Finite_strain_system_increment_solution.reinit(dof_handler.n_dofs());
small_strain_system_rhs.reinit(dof_handler.n_dofs());
Finite_strain_system_rhs.reinit(dof_handler.n_dofs());
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, false);
sparsity_pattern.copy_from(dsp);
small_strain_stiff_matrix.reinit(sparsity_pattern);
Finite_strain_tangent_stiff_matrix.reinit(sparsity_pattern);
end=clock();
std::cout<<"the time of set up system "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::setup_quadrature_point_history()
{
clock_t start, end;
start=clock();
std::cout<<"start to set up quadrature point history"<<std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
std::vector<Parameters::PointHistory> tmp;
quadrature_point_history.swap(tmp);
const unsigned int n_q_points = quadrature_formula.size();
quadrature_point_history.resize(triangulation.n_active_cells() * n_q_points);
long int the_number_of_quadrature_point = 0;//表示总的高斯点数量。
for (auto& cell : triangulation.active_cell_iterators())
{
fe_values.reinit(cell);
std::vector<Point<dim>> points = fe_values.get_quadrature_points();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
double radius = points[q_point][0] * points[q_point][0] + points[q_point][1] * points[q_point][1] + points[q_point][2] * points[q_point][2];
if (radius <= 1)
{
quadrature_point_history[the_number_of_quadrature_point].order_parameters(0) = 1;
}
the_number_of_quadrature_point++;
}
}
std::cout<<"the_number_of_quadrature_point= "<<the_number_of_quadrature_point<<std::endl;
end=clock();
std::cout<<"the time of set up quadrature point history "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::update_quadrature_point_history()
{
clock_t start, end;
start=clock();
std::cout<<"update quadrature point history"<<std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Vector<double> nodle_displace(dofs_per_cell);
long int the_number_of_quadrature_point = 0;
for (const auto& cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell->get_dof_indices(local_dof_indices);
nodle_displace = 0;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
nodle_displace(i) = Finite_strain_system_solution(local_dof_indices[i]);
FullMatrix<double> nodle_displace_x_y_z(3, 8);
for (int i = 0; i < 3; i++)
for (unsigned int j = 0; j < dofs_per_cell / 3; j++)
nodle_displace_x_y_z(i, j) = nodle_displace(3 * j + i);
//std::vector<Point<dim> > points = fe_values.get_quadrature_points();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
FullMatrix<double> shape_function_grad(3, dofs_per_cell / 3);
FullMatrix<double> F(3, 3);
F = 0;
for (int i = 0; i < 3; i++)
for (unsigned int j = 0; j < dofs_per_cell / 3; j++)
shape_function_grad(i, j) = fe_values.shape_grad(3 * j , q_point)[i];
FullMatrix<double> d_u_d_x(3, 3);
shape_function_grad.mTmult(d_u_d_x, nodle_displace_x_y_z);
F.add(1, deformation_variable.I, 1, d_u_d_x);
/*******************************/
FullMatrix<double> Fi_inv(3,3);
Fi_inv=0;
if(quadrature_point_history[the_number_of_quadrature_point].order_parameters(0)==1)
{
//Fi_inv_T=deformation_variable.Fi_inv_T;
Fi_inv=deformation_variable.Fi_inv;
}else{
//Fi_inv_T=deformation_variable.I;
Fi_inv=deformation_variable.I;
}
FullMatrix<double> Fe(3, 3);
F.mmult(Fe,Fi_inv);
FullMatrix<double> tao_Ee(3, 3);
FullMatrix<double> E(3, 3);
tao_Ee=0;
E=0;
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
for (int k = 0; k < 3; k++)
{
tao_Ee(i, j) += Fe(k, i) * Fe(k, j) ;
E(i, j) += F(k, i) * F(k, j);
}
tao_Ee(i, j) = tao_Ee(i, j)/2 - deformation_variable.I(i, j) / 2;
E(i, j)=E(i, j)/2 - deformation_variable.I(i, j) / 2;
}
}
Vector<double> tao_Ee_vector(6);
tao_Ee_vector=0;
Vector<double> E_vector(6);
E_vector=0;
tao_Ee_vector(0) = tao_Ee(0, 0);
tao_Ee_vector(1) = tao_Ee(1, 1);
tao_Ee_vector(2) = tao_Ee(2, 2);
tao_Ee_vector(3) = 2 * tao_Ee(1, 2);
tao_Ee_vector(4) = 2 * tao_Ee(0, 2);
tao_Ee_vector(5) = 2 * tao_Ee(0, 1);
/****************************/
E_vector(0) = E(0, 0);
E_vector(1) = E(1, 1);
E_vector(2) = E(2, 2);
E_vector(3) = 2 * E(1, 2);
E_vector(4) = 2 * E(0, 2);
E_vector(5) = 2 * E(0, 1);
/****************************/
for (int k = 0; k < 6; k++)
quadrature_point_history[the_number_of_quadrature_point].strain_vector(k) = E_vector(k);
for (int i = 0; i < 6; i++)
{
long double sum = 0;
for (int j = 0; j < 6; j++)
sum += deformation_variable.elastic_matrix(i, j) * tao_Ee_vector(j);
quadrature_point_history[the_number_of_quadrature_point].stress_vector(i) = sum;
}
the_number_of_quadrature_point++;
}
}
end=clock();
std::cout<<"the time of update quadrature point history "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::assemble_small_strain_system()
{
clock_t start, end;
start=clock();
std::cout<<"start to assemble small strain system"<<std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.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);
FullMatrix<double> B(6, dofs_per_cell);
long int the_number_of_quadrature_point = 0;//用来编号高斯点。
for (const auto& cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0;
cell_rhs = 0;
B = 0;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (unsigned int ii = 1; ii <= dofs_per_cell/3; ii++)
{
B(0, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
B(1, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(2, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(3, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(3, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(4, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(4, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
B(5, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(5, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
}
/*************************************下面计算B^T*C***********************************/
FullMatrix<double> temp1(dofs_per_cell, 6);
temp1=0;
B.Tmmult(temp1, deformation_variable.elastic_matrix);
FullMatrix<double> temp2(dofs_per_cell, dofs_per_cell);
temp2=0;
temp1.mmult(temp2, B);
cell_matrix.add(fe_values.JxW(q_point), temp2);
Vector<double> temp3(dofs_per_cell);
temp3=0;
temp1.vmult(temp3, deformation_variable.engstrain_vector);
cell_rhs.add(fe_values.JxW(q_point) * quadrature_point_history[the_number_of_quadrature_point].order_parameters(0), temp3);
the_number_of_quadrature_point++;
}
/*******************************************************************************/
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, small_strain_stiff_matrix, small_strain_system_rhs);
}
end=clock();
std::cout<<"the time of assemble small strain stiff matrix "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::output_small_strain_displacement()
{
clock_t start, end;
start=clock();
std::cout<<"start to output small strain displacement"<<std::endl;
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
std::vector<std::string> solution_names;
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
solution_names.emplace_back("z_displacement");
data_out.add_data_vector(small_strain_system_solution, solution_names);
data_out.build_patches();
std::ofstream output("small_strain_displacement.vtk");
data_out.write_vtk(output);
end=clock();
std::cout<<"the time of output small strain displacement "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::output_small_strain_elastic_field()
{
clock_t start, end;
start=clock();
std::cout<<"start to output small strain elastic field"<<std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Vector<double> nodle_displace(dofs_per_cell);
FullMatrix<double> B(6, dofs_per_cell);
FullMatrix<double> quadrature_strain(n_q_points, 6);
FullMatrix<double> quadrature_stress(n_q_points, 6);
FullMatrix<double> H(n_q_points, 8);
FullMatrix<double> H_T_multi_H(8, 8);
FullMatrix<double> H_T_multi_stress(8, 6);
FullMatrix<double> H_T_multi_strain(8, 6);
Vector<double> stress_xx_rhs(8);
Vector<double> stress_yy_rhs(8);
Vector<double> stress_zz_rhs(8);
Vector<double> stress_yz_rhs(8);
Vector<double> stress_xz_rhs(8);
Vector<double> stress_xy_rhs(8);
/*************************************************************/
Vector<double> strain_xx_rhs(8);
Vector<double> strain_yy_rhs(8);
Vector<double> strain_zz_rhs(8);
Vector<double> strain_yz_rhs(8);
Vector<double> strain_xz_rhs(8);
Vector<double> strain_xy_rhs(8);
/*************************************************************/
Vector<double> solution_stress_xx(8);
Vector<double> solution_stress_yy(8);
Vector<double> solution_stress_zz(8);
Vector<double> solution_stress_yz(8);
Vector<double> solution_stress_xz(8);
Vector<double> solution_stress_xy(8);
/*************************************************************/
Vector<double> solution_strain_xx(8);
Vector<double> solution_strain_yy(8);
Vector<double> solution_strain_zz(8);
Vector<double> solution_strain_yz(8);
Vector<double> solution_strain_xz(8);
Vector<double> solution_strain_xy(8);
/*************************************************************/
//这里定义数组，用来保存单元节点的应力和应变。
Vector<double> normal_stress(dof_handler.n_dofs());
normal_stress = 0;
Vector<double> normal_strain(dof_handler.n_dofs());
normal_strain = 0;
Vector<double> shear_stress_p(dof_handler.n_dofs());
shear_stress_p = 0;
Vector<double> shear_strain_p(dof_handler.n_dofs());
shear_strain_p = 0;
long long int the_number_of_quadrature_point = 0;
for (const auto& cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell->get_dof_indices(local_dof_indices);
nodle_displace = 0;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
nodle_displace(i) = small_strain_system_solution(local_dof_indices[i]);
B = 0;
quadrature_strain = 0;
quadrature_stress = 0;
H = 0;
std::vector<Point<dim> > points = fe_values.get_quadrature_points();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (int k = 0; k < 8; k++)
H(q_point, k) = fe_values.shape_value(3 * k, q_point); //即H矩阵。
for (int ii = 1; ii <= 8; ii++)
{
B(0, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
B(1, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(2, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(3, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(3, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(4, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[2];
B(4, 3 * ii - 1) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
B(5, 3 * ii - 3) = fe_values.shape_grad(3 * ii - 1, q_point)[1];
B(5, 3 * ii - 2) = fe_values.shape_grad(3 * ii - 1, q_point)[0];
}
for (int k = 0; k < 6; k++)
{
long double sum = 0;
for (unsigned int j = 0; j < dofs_per_cell;j++)
sum += B(k, j) * nodle_displace(j);
quadrature_strain(q_point, k) = sum;
}
//上面获得高斯点处的应变值
for (int i = 0; i < 6; i++)
{
long double sum = 0;
for (int j = 0; j < 6; j++)
sum += deformation_variable.elastic_matrix(i, j)*(quadrature_strain(q_point, j) - quadrature_point_history[the_number_of_quadrature_point].order_parameters(0) * deformation_variable.engstrain_vector(j));
quadrature_stress(q_point, i) = sum;
}
//上面获得高斯点处的应力值
the_number_of_quadrature_point++;
}
H_T_multi_H = 0;
H.Tmmult(H_T_multi_H, H);
H_T_multi_stress = 0;
H.Tmmult(H_T_multi_stress, quadrature_stress);
H_T_multi_strain = 0;
H.Tmmult(H_T_multi_strain, quadrature_strain);
/***************************************************/
strain_xx_rhs = 0;
strain_yy_rhs = 0;
strain_zz_rhs = 0;
strain_yz_rhs = 0;
strain_xz_rhs = 0;
strain_xy_rhs = 0;
/***************************************************/
stress_xx_rhs = 0;
stress_yy_rhs = 0;
stress_zz_rhs = 0;
stress_yz_rhs = 0;
stress_xz_rhs = 0;
stress_xy_rhs = 0;
/***************************************************/
for (int i = 0; i < 8; i++)
{
strain_xx_rhs(i) = H_T_multi_strain(i, 0);
strain_yy_rhs(i) = H_T_multi_strain(i, 1);
strain_zz_rhs(i) = H_T_multi_strain(i, 2);
strain_yz_rhs(i) = H_T_multi_strain(i, 3);
strain_xz_rhs(i) = H_T_multi_strain(i, 4);
strain_xy_rhs(i) = H_T_multi_strain(i, 5);
/******************************************/
stress_xx_rhs(i) = H_T_multi_stress(i, 0);
stress_yy_rhs(i) = H_T_multi_stress(i, 1);
stress_zz_rhs(i) = H_T_multi_stress(i, 2);
stress_yz_rhs(i) = H_T_multi_stress(i, 3);
stress_xz_rhs(i) = H_T_multi_stress(i, 4);
stress_xy_rhs(i) = H_T_multi_stress(i, 5);
}
/*****************************************************/
solution_strain_xx = 0;
solution_strain_yy = 0;
solution_strain_zz = 0;
solution_strain_yz = 0;
solution_strain_xz = 0;
solution_strain_xy = 0;
/*********************************************************/
solution_stress_xx = 0;
solution_stress_yy = 0;
solution_stress_zz = 0;
solution_stress_yz = 0;
solution_stress_xz = 0;
solution_stress_xy = 0;
/*****************************************************************************************/
SolverControl solver_control(1000, 1e-16);
SolverCG<> solver(solver_control);
solver.solve(H_T_multi_H, solution_strain_xx, strain_xx_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_yy, strain_yy_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_zz, strain_zz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_yz, strain_yz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_xz, strain_xz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_xy, strain_xy_rhs, PreconditionIdentity());
/*****************************************************************************************/
solver.solve(H_T_multi_H, solution_stress_xx, stress_xx_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_yy, stress_yy_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_zz, stress_zz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_yz, stress_yz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_xz, stress_xz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_xy, stress_xy_rhs, PreconditionIdentity());
/*****************************************************************************************/
for (int i = 0; i < 8; i++)
{
int a = local_dof_indices[3 * i];
int b = local_dof_indices[3 * i + 1];
int c = local_dof_indices[3 * i + 2];
/*****************************************************************************************/
normal_stress(a) += solution_stress_xx(i);
normal_stress(b) += solution_stress_yy(i);
normal_stress(c) += solution_stress_zz(i);
shear_stress_p(a) += solution_stress_yz(i);
shear_stress_p(b) += solution_stress_xz(i);
shear_stress_p(c) += solution_stress_xy(i);
/*****************************************************************************************/
normal_strain(a) += solution_strain_xx(i);
normal_strain(b) += solution_strain_yy(i);
normal_strain(c) += solution_strain_zz(i);
shear_strain_p(a) += solution_stress_yz(i);
shear_strain_p(b) += solution_strain_xz(i);
shear_strain_p(c) += solution_strain_xy(i);
}
}
//上面单元循环就此结束
for (unsigned int ii = 0; ii < dof_handler. n_dofs(); ii++)
{
normal_strain(ii) = normal_strain(ii) / 8;
normal_stress(ii) = normal_stress(ii) / 8;
shear_stress_p(ii) = shear_stress_p(ii) / 8;
shear_strain_p(ii) = shear_strain_p(ii) / 8;
}
//上面一个循环的目的是将这些应力应变放到和位移一样的序列上去，用于输出
/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/
//下面正式输出弹性场，应力和应变场
std::string strain = "small_strain_normal_strain_";
std::string stress = "small_strain_normal_stress_";
std::string shear_stress = "small_strain_shear_stress_";
std::string shear_strain = "small_strain_shear_strain_";
strain += ".vtk";
stress += ".vtk";
shear_stress += ".vtk";
shear_strain += ".vtk";
std::ofstream output_strain(strain.c_str());
std::ofstream output_stress(stress.c_str());
std::ofstream output_shear_stress(shear_stress.c_str());
std::ofstream output_shear_strain(shear_strain.c_str());
DataOut<dim> strain_out, stress_out, shear_stress_out, shear_strain_out;
strain_out.attach_dof_handler(dof_handler);
stress_out.attach_dof_handler(dof_handler);
shear_stress_out.attach_dof_handler(dof_handler);
shear_strain_out.attach_dof_handler(dof_handler);
std::vector<std::string> strain_name, stress_name, shear_stress_name, shear_strain_name;
strain_name.push_back("strain_xx");
strain_name.push_back("strain_yy");
strain_name.push_back("strain_zz");
stress_name.push_back("stress_xx");
stress_name.push_back("stress_yy");
stress_name.push_back("stress_zz");
shear_stress_name.push_back("stress_yz");
shear_stress_name.push_back("stress_xz");
shear_stress_name.push_back("stress_xy");
shear_strain_name.push_back("strain_yz");
shear_strain_name.push_back("strain_xz");
shear_strain_name.push_back("strain_xy");
strain_out.add_data_vector(normal_strain, strain_name);
stress_out.add_data_vector(normal_stress, stress_name);
shear_stress_out.add_data_vector(shear_stress_p, shear_stress_name);
shear_strain_out.add_data_vector(shear_strain_p, shear_strain_name);
strain_out.build_patches();
stress_out.build_patches();
shear_stress_out.build_patches();
shear_strain_out.build_patches();
strain_out.write_vtk(output_strain);
stress_out.write_vtk(output_stress);
shear_stress_out.write_vtk(output_shear_stress);
shear_strain_out.write_vtk(output_shear_strain);
end=clock();
std::cout<<"the time of output small strain elastic field "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::assemble_Finite_strain_system()
{
clock_t start, end;
start=clock();
std::cout<<"start to assemble Finite strain system"<<std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> cell_rhs(dofs_per_cell);
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> tangent_stiff(dofs_per_cell, dofs_per_cell);
FullMatrix<double> geometric_stiff(dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Vector<double> nodle_displace(dofs_per_cell);
FullMatrix<double> B(6, dofs_per_cell);
long int the_number_of_quadrature_point = 0;
for (const auto& cell : dof_handler.active_cell_iterators())
{
cell_rhs = 0;
cell_matrix = 0;
tangent_stiff=0;
geometric_stiff=0;
fe_values.reinit(cell);
cell->get_dof_indices(local_dof_indices);
nodle_displace = 0;
for (unsigned int i = 0; i < dofs_per_cell; ++i)
nodle_displace(i) = Finite_strain_system_solution(local_dof_indices[i]);
FullMatrix<double> nodle_displace_x_y_z(3, 8);
for (int i = 0; i < 3; i++)
for (unsigned int j = 0; j < dofs_per_cell / 3; j++)
nodle_displace_x_y_z(i, j) = nodle_displace(3 * j + i);
B = 0;
// std::vector<Point<dim> > points = fe_values.get_quadrature_points();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
FullMatrix<double> shape_function_grad(3, dofs_per_cell / 3);
shape_function_grad=0;
FullMatrix<double> F(3, 3);
F = 0;
for (int i = 0; i < 3; i++)
for (unsigned int j = 0; j < dofs_per_cell/3; j++)
shape_function_grad(i, j) = fe_values.shape_grad(3 * j, q_point)[i];
FullMatrix<double> d_u_d_x(3, 3);
shape_function_grad.mTmult(d_u_d_x, nodle_displace_x_y_z);
F.add(1, deformation_variable.I, 1, d_u_d_x);
//FullMatrix<double> Fi_inv_T(3,3);
FullMatrix<double> Fi_inv(3,3);
double Ji;
if(quadrature_point_history[the_number_of_quadrature_point].order_parameters(0)==1)
{
//Fi_inv_T=deformation_variable.Fi_inv_T;
Fi_inv=deformation_variable.Fi_inv;
Ji=deformation_variable.Ji;
}else{
//Fi_inv_T=deformation_variable.I;
Fi_inv=deformation_variable.I;
Ji=1;
}
for (unsigned int a = 0; a < dofs_per_cell / 3; a++)
{
for (int l = 0; l < 3; l++)
{
for (int m = 0; m < 3; m++)
{
for(int k=0;k<3;k++)
{
B(0, 3 * a+k) += Fi_inv(l, 0)* Fi_inv(m, 0)*(shape_function_grad(l, a)*F(k,m)/2+shape_function_grad(m, a)*F(k,l)/2);
B(1, 3 * a+k) += Fi_inv(l, 1)* Fi_inv(m, 1)*(shape_function_grad(l, a)*F(k,m)/2+shape_function_grad(m, a)*F(k,l)/2);
B(2, 3 * a+k) += Fi_inv(l, 2)* Fi_inv(m, 2)*(shape_function_grad(l, a)*F(k,m)/2+shape_function_grad(m, a)*F(k,l)/2);
B(3, 3 * a+k) += Fi_inv(l, 1)* Fi_inv(m, 2)*(shape_function_grad(l, a)*F(k,m)+shape_function_grad(m, a)*F(k,l));
B(4, 3 * a+k) += Fi_inv(l, 0)* Fi_inv(m, 2)*(shape_function_grad(l, a)*F(k,m)+shape_function_grad(m, a)*F(k,l));
B(5, 3 * a+k) += Fi_inv(l, 0)* Fi_inv(m, 1)*(shape_function_grad(l, a)*F(k,m)+shape_function_grad(m, a)*F(k,l));
}
}
}
}
Vector<double> temp(dofs_per_cell);
temp=0;
B.Tvmult(temp,quadrature_point_history[the_number_of_quadrature_point].stress_vector);
cell_rhs.add(fe_values.JxW(q_point)*Ji, temp);
FullMatrix<double> temp1(dofs_per_cell,6);
temp1=0;
B.Tmmult(temp1,deformation_variable.elastic_matrix);
FullMatrix<double> temp2(dofs_per_cell,dofs_per_cell);
temp2=0;
temp1.mmult(temp2,B);
tangent_stiff.add(fe_values.JxW(q_point)*Ji,temp2);
FullMatrix<double> temp3(8,8);
temp3=0;
FullMatrix<double> temp4(dofs_per_cell,dofs_per_cell);
temp4=0;
FullMatrix<double> tao_S(3,3);
tao_S(0,0)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(0);
tao_S(0,1)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(5);
tao_S(0,2)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(4);
tao_S(1,0)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(5);
tao_S(1,1)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(1);
tao_S(1,2)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(3);
tao_S(2,0)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(4);
tao_S(2,1)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(3);
tao_S(2,2)=quadrature_point_history[the_number_of_quadrature_point].stress_vector(2);
for (int a=0;a<8;a++)
for (int b=0;b<8;b++)
for (int p=0;p<3;p++)
for (int q=0;q<3;q++)
for (int i=0;i<3;i++)
for (int j=0;j<3;j++)
temp3(a,b)+=Fi_inv(i,p)*Fi_inv(j,q)*tao_S(p,q)*shape_function_grad(i,a)*shape_function_grad(j,b);
for (int a=0;a<8;a++)
for (int b=0;b<8;b++)
for(int i=0;i<3;i++)
temp4(3*a+i,3*b+i)=temp3(a,b);
geometric_stiff.add(fe_values.JxW(q_point)*Ji,temp4);
the_number_of_quadrature_point++;
}
cell_matrix.add(1,tangent_stiff,1,geometric_stiff);
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, Finite_strain_tangent_stiff_matrix, Finite_strain_system_rhs);
}
for (unsigned int i = 0; i < dof_handler.n_dofs(); i++)
Finite_strain_system_rhs(i) = -Finite_strain_system_rhs(i);
end=clock();
std::cout<<"the time of assemble Finite strain system "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::output_Finite_strain_displacement(int number_of_iteration)
{
clock_t start, end;
start=clock();
std::cout<<"start to output Finite strain displacement"<<std::endl;
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
std::vector<std::string> solution_names;
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
solution_names.emplace_back("z_displacement");
data_out.add_data_vector(Finite_strain_system_solution, solution_names);
data_out.build_patches();
std::ofstream output("Finite_strain_dispancement_"+std::to_string(number_of_iteration)+".vtk");
data_out.write_vtk(output);
end=clock();
std::cout<<"the time of output Finite_ train displacement "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::output_Finite_strain_elastic_field(int number_of_iteration)
{
clock_t start, end;
start = clock();
std::cout << "start to output Finite strain elastic field" << std::endl;
QGauss<dim> quadrature_formula(parameters.quad_order);
FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
FullMatrix<double> quadrature_strain(n_q_points, 6);
FullMatrix<double> quadrature_stress(n_q_points, 6);
FullMatrix<double> H(n_q_points, 8);
FullMatrix<double> H_T_multi_H(8, 8);
FullMatrix<double> H_T_multi_stress(8, 6);
FullMatrix<double> H_T_multi_strain(8, 6);
Vector<double> stress_xx_rhs(8);
Vector<double> stress_yy_rhs(8);
Vector<double> stress_zz_rhs(8);
Vector<double> stress_yz_rhs(8);
Vector<double> stress_xz_rhs(8);
Vector<double> stress_xy_rhs(8);
/*************************************************************/
Vector<double> strain_xx_rhs(8);
Vector<double> strain_yy_rhs(8);
Vector<double> strain_zz_rhs(8);
Vector<double> strain_yz_rhs(8);
Vector<double> strain_xz_rhs(8);
Vector<double> strain_xy_rhs(8);
/*************************************************************/
Vector<double> solution_stress_xx(8);
Vector<double> solution_stress_yy(8);
Vector<double> solution_stress_zz(8);
Vector<double> solution_stress_yz(8);
Vector<double> solution_stress_xz(8);
Vector<double> solution_stress_xy(8);
/*************************************************************/
Vector<double> solution_strain_xx(8);
Vector<double> solution_strain_yy(8);
Vector<double> solution_strain_zz(8);
Vector<double> solution_strain_yz(8);
Vector<double> solution_strain_xz(8);
Vector<double> solution_strain_xy(8);
/*************************************************************/
//这里定义数组，用来保存单元节点的应力和应变。
Vector<double> normal_stress(dof_handler.n_dofs());
normal_stress = 0;
Vector<double> normal_strain(dof_handler.n_dofs());
normal_strain = 0;
Vector<double> shear_stress_p(dof_handler.n_dofs());
shear_stress_p = 0;
Vector<double> shear_strain_p(dof_handler.n_dofs());
shear_strain_p = 0;
long int the_number_of_quadrature_point = 0;
for (const auto& cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
quadrature_strain = 0;
quadrature_stress = 0;
H = 0;
cell->get_dof_indices(local_dof_indices);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (int k = 0; k < 8; k++)
H(q_point, k) = fe_values.shape_value(3 * k, q_point); //即H矩阵。
for (int i = 0; i < 6; i++)
{
quadrature_stress(q_point, i) = quadrature_point_history[the_number_of_quadrature_point].stress_vector(i);
quadrature_strain(q_point, i) = quadrature_point_history[the_number_of_quadrature_point].strain_vector(i);
}
the_number_of_quadrature_point++;
}
H_T_multi_H = 0;
H.Tmmult(H_T_multi_H, H);
H_T_multi_stress = 0;
H.Tmmult(H_T_multi_stress, quadrature_stress);
H_T_multi_strain = 0;
H.Tmmult(H_T_multi_strain, quadrature_strain);
/***************************************************/
strain_xx_rhs = 0;
strain_yy_rhs = 0;
strain_zz_rhs = 0;
strain_yz_rhs = 0;
strain_xz_rhs = 0;
strain_xy_rhs = 0;
/***************************************************/
stress_xx_rhs = 0;
stress_yy_rhs = 0;
stress_zz_rhs = 0;
stress_yz_rhs = 0;
stress_xz_rhs = 0;
stress_xy_rhs = 0;
/***************************************************/
for (int i = 0; i < 8; i++)
{
strain_xx_rhs(i) = H_T_multi_strain(i, 0);
strain_yy_rhs(i) = H_T_multi_strain(i, 1);
strain_zz_rhs(i) = H_T_multi_strain(i, 2);
strain_yz_rhs(i) = H_T_multi_strain(i, 3);
strain_xz_rhs(i) = H_T_multi_strain(i, 4);
strain_xy_rhs(i) = H_T_multi_strain(i, 5);
/******************************************/
stress_xx_rhs(i) = H_T_multi_stress(i, 0);
stress_yy_rhs(i) = H_T_multi_stress(i, 1);
stress_zz_rhs(i) = H_T_multi_stress(i, 2);
stress_yz_rhs(i) = H_T_multi_stress(i, 3);
stress_xz_rhs(i) = H_T_multi_stress(i, 4);
stress_xy_rhs(i) = H_T_multi_stress(i, 5);
}
/*****************************************************/
solution_strain_xx = 0;
solution_strain_yy = 0;
solution_strain_zz = 0;
solution_strain_yz = 0;
solution_strain_xz = 0;
solution_strain_xy = 0;
/*********************************************************/
solution_stress_xx = 0;
solution_stress_yy = 0;
solution_stress_zz = 0;
solution_stress_yz = 0;
solution_stress_xz = 0;
solution_stress_xy = 0;
/*****************************************************************************************/
SolverControl solver_control(1000, 1e-16);
SolverCG<> solver(solver_control);
solver.solve(H_T_multi_H, solution_strain_xx, strain_xx_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_yy, strain_yy_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_zz, strain_zz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_yz, strain_yz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_xz, strain_xz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_strain_xy, strain_xy_rhs, PreconditionIdentity());
/*****************************************************************************************/
solver.solve(H_T_multi_H, solution_stress_xx, stress_xx_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_yy, stress_yy_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_zz, stress_zz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_yz, stress_yz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_xz, stress_xz_rhs, PreconditionIdentity());
solver.solve(H_T_multi_H, solution_stress_xy, stress_xy_rhs, PreconditionIdentity());
/*****************************************************************************************/
for (int i = 0; i < 8; i++)
{
int a = local_dof_indices[3 * i];
int b = local_dof_indices[3 * i + 1];
int c = local_dof_indices[3 * i + 2];
/*****************************************************************************************/
normal_stress(a) += solution_stress_xx(i);
normal_stress(b) += solution_stress_yy(i);
normal_stress(c) += solution_stress_zz(i);
shear_stress_p(a) += solution_stress_yz(i);
shear_stress_p(b) += solution_stress_xz(i);
shear_stress_p(c) += solution_stress_xy(i);
/*****************************************************************************************/
normal_strain(a) += solution_strain_xx(i);
normal_strain(b) += solution_strain_yy(i);
normal_strain(c) += solution_strain_zz(i);
shear_strain_p(a) += solution_stress_yz(i);
shear_strain_p(b) += solution_strain_xz(i);
shear_strain_p(c) += solution_strain_xy(i);
}
}
//上面单元循环就此结束
for (unsigned int ii = 0; ii < dof_handler.n_dofs(); ii++)
{
normal_strain(ii) = normal_strain(ii) / 8;
normal_stress(ii) = normal_stress(ii) / 8;
shear_stress_p(ii) = shear_stress_p(ii) / 8;
shear_strain_p(ii) = shear_strain_p(ii) / 8;
}
//上面一个循环的目的是将这些应力应变放到和位移一样的序列上去，用于输出
/***********************************************************************************/
/***********************************************************************************/
/***********************************************************************************/
//下面正式输出弹性场，应力和应变场
std::string strain = "Finite_strain_normal_strain_";
std::string stress = "Finite_strain_normal_stress_";
std::string shear_stress = "Finite_strain_shear_stress_";
std::string shear_strain = "Finite_strain_shear_strain_";
strain +=std::to_string(number_of_iteration)+".vtk";
stress +=std::to_string(number_of_iteration) + ".vtk";
shear_stress +=std::to_string(number_of_iteration) + ".vtk";
shear_strain +=std::to_string(number_of_iteration) + ".vtk";
std::ofstream output_strain(strain.c_str());
std::ofstream output_stress(stress.c_str());
std::ofstream output_shear_stress(shear_stress.c_str());
std::ofstream output_shear_strain(shear_strain.c_str());
DataOut<dim> strain_out, stress_out, shear_stress_out, shear_strain_out;
strain_out.attach_dof_handler(dof_handler);
stress_out.attach_dof_handler(dof_handler);
shear_stress_out.attach_dof_handler(dof_handler);
shear_strain_out.attach_dof_handler(dof_handler);
std::vector<std::string> strain_name, stress_name, shear_stress_name, shear_strain_name;
strain_name.push_back("strain_xx");
strain_name.push_back("strain_yy");
strain_name.push_back("strain_zz");
stress_name.push_back("stress_xx");
stress_name.push_back("stress_yy");
stress_name.push_back("stress_zz");
shear_stress_name.push_back("stress_yz");
shear_stress_name.push_back("stress_xz");
shear_stress_name.push_back("stress_xy");
shear_strain_name.push_back("strain_yz");
shear_strain_name.push_back("strain_xz");
shear_strain_name.push_back("strain_xy");
strain_out.add_data_vector(normal_strain, strain_name);
stress_out.add_data_vector(normal_stress, stress_name);
shear_stress_out.add_data_vector(shear_stress_p, shear_stress_name);
shear_strain_out.add_data_vector(shear_strain_p, shear_strain_name);
strain_out.build_patches();
stress_out.build_patches();
shear_stress_out.build_patches();
shear_strain_out.build_patches();
strain_out.write_vtk(output_strain);
stress_out.write_vtk(output_stress);
shear_stress_out.write_vtk(output_shear_stress);
shear_strain_out.write_vtk(output_shear_strain);
end = clock();
std::cout << "the time of output Finite strain elastic field " << double(end - start) / CLOCKS_PER_SEC << "s" << std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::solve_linear(int judge)
{
clock_t start, end;
start=clock();
std::cout<<"start to solve linear"<<std::endl;
SolverControl solver_control(parameters.iteration_number, parameters.err);
SolverCG<Vector<double>> cg(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
if (judge == 0)
{
preconditioner.initialize(small_strain_stiff_matrix, 1.2);
cg.solve(small_strain_stiff_matrix, small_strain_system_solution, small_strain_system_rhs, preconditioner);
constraints.distribute(small_strain_system_rhs);
}
if (judge == 1)
{
Finite_strain_system_increment_solution=0;
preconditioner.initialize(Finite_strain_tangent_stiff_matrix, 1.2);
cg.solve(Finite_strain_tangent_stiff_matrix, Finite_strain_system_increment_solution, Finite_strain_system_rhs, preconditioner);
constraints.distribute(Finite_strain_system_increment_solution);
}
end=clock();
std::cout<<"the time of solve linear "<<double(end-start)/CLOCKS_PER_SEC<<"s"<<std::endl;
}
template <int dim>
void Finite_strain_problem<dim>::run()
{
parameter_setup();
create_grid();
setup_system();
setup_quadrature_point_history();
assemble_small_strain_system();
solve_linear(0);
output_small_strain_displacement();
output_small_strain_elastic_field();
solve_nolinear();
}
template <int dim>
void Finite_strain_problem<dim>::solve_nolinear()
{
for (unsigned int i = 0; i < dof_handler.n_dofs(); i++)
Finite_strain_system_solution(i) = small_strain_system_solution(i);
update_quadrature_point_history();
int number_of_iteration = 0;
long double solve_err = 100;
while((number_of_iteration == 0) || (solve_err > 1e-12))
{
Finite_strain_system_rhs=0;
Finite_strain_tangent_stiff_matrix=0;
assemble_Finite_strain_system();
solve_err = Finite_strain_system_rhs.l2_norm();
std::cout<<"solve_err="<<solve_err<<std::endl;
solve_linear(1);
for (unsigned int i = 0; i < dof_handler.n_dofs(); i++)
Finite_strain_system_solution(i)+= Finite_strain_system_increment_solution(i);
update_quadrature_point_history();
output_Finite_strain_displacement(number_of_iteration);
output_Finite_strain_elastic_field(number_of_iteration);
number_of_iteration++;
}
}
}
int main()
{
try
{
Finite_strain_PFM::Finite_strain_problem<3> finite_strain_problem_3d("parameters.prm");
finite_strain_problem_3d.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;
}
# Listing of Parameters
# ---------------------
#*******************YangPing***************************#
subsection Finite strain PFM
#############################################################
set Polynomial degree = 1
set Quadrature order = 2
#############################################################
set C11=3.5e+11
set C12=1.5e+11
set C13=1.5e+11
set C14=0
set C15=0
set C16=0
set C22=3.5e+11
set C23=1.5e+11
set C24=0
set C25=0
set C26=0
set C33=3.5e+11
set C34=0
set C35=0
set C36=0
set C44=1e+11
set C45=0
set C46=0
set C55=1e+11
set C56=0
set C66=1e+11
#############################################################
set err=1e-12
set iteration number=2000
#############################################################
set length left=-5
set length right=5
set refine number=6
#############################################################
set F11=0.93
set F12=0
set F13=0
set F21=0
set F22=0.93
set F23=0
set F31=0
set F32=0
set F33=0.93
#############################################################
end
扫一扫
4323
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
