deal-II软件计算Eshelby夹杂程序，以及matlab处理数据程序

最新推荐文章于 2024-07-12 15:51:28 发布
最新推荐文章于 2024-07-12 15:51:28 发布
阅读量484 收藏 1
点赞数
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/py753159123456/article/details/114987537
版权

deal-ii程序

#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/constraint_matrix.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 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;
    ConstraintMatrix     hanging_node_constraints;
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
  };
  /
  /
  /
  /
  /
  //此处添加右端项 ，如果是本征应变问题，同样在这里添加。 
  /
  /
  /
  /
  /
  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);
    hanging_node_constraints.clear ();
    DoFTools::make_hanging_node_constraints (dof_handler,
                                             hanging_node_constraints);
    hanging_node_constraints.close ();
    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern(dof_handler,
                                    dsp,
                                    hanging_node_constraints,
                                    /*keep_constrained_dofs = */ true);
    sparsity_pattern.copy_from (dsp);
    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);//选取高斯积分点 
    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);//获取全局节点编号。 
    /*std::vector<double>     lambda_values (n_q_points);
    std::vector<double>     mu_values (n_q_points);
    ConstantFunction<dim> lambda(1.), mu(1.);*///此处是设定弹性张量，一般为常数，故不用此方法。
	double nu=0.3,mu=100000000000;
    double lambda=2*mu*nu/(1-2*nu); 
    double temperature_epsilon[6]={0.000001,0.000001,0.000001,0,0,0};
    double sum=temperature_epsilon[0]+temperature_epsilon[1]+temperature_epsilon[2];
    std::vector<Tensor<1, dim> > rhs_values (n_q_points);//将高斯积分点传入右手边计算函数，积分得到右手边的高斯积分点处的值。 
    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();//单元迭代器。 
    for (; cell!=endc; ++cell)
      {
        cell_matrix = 0;//单元矩阵赋值为0. 
        cell_rhs = 0;//单元右边项赋值为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);//此步骤待考虑。
        std::vector<Point<dim> > points=fe_values.get_quadrature_points();
        //int number_quadrature_points=points.size();
        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)
                        +
                        (fe_values.shape_grad(i,q_point)[component_j] *
                         fe_values.shape_grad(j,q_point)[component_i] *
                         mu)
                        +
                        ((component_i == component_j) ?
                         (fe_values.shape_grad(i,q_point) *
                          fe_values.shape_grad(j,q_point) *
                          mu)  :
                         0)
                      )
                      *
                      fe_values.JxW(q_point);
                  }
              }
          }
        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 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(points[q_point][0]<1&&points[q_point][0]>-1)     
           // std::cout<<"radius="<<points[q_point][0]<<std::endl;
            if(radius<=1){
            cell_rhs(i) +=((2*mu*fe_values.shape_grad(i,q_point)[component_i]*temperature_epsilon[component_i])
			    +lambda*sum*fe_values.shape_grad(i,q_point)[component_i])*fe_values.JxW(q_point);
			   }                             
          }
          }
        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);
          }
      }
    hanging_node_constraints.condense (system_matrix);
    hanging_node_constraints.condense (system_rhs);
    std::map<types::global_dof_index,double> boundary_values;
    VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
                                              ZeroFunction<dim>(dim),
                                              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-16);
    SolverCG<>              cg (solver_control);
    PreconditionSSOR<> preconditioner;
    preconditioner.initialize(system_matrix, 1.2);
    cg.solve (system_matrix, solution, system_rhs,
              preconditioner);
    hanging_node_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>(2),
                                        typename FunctionMap<dim>::type(),
                                        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 < 10, 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);*/
    // 
    // 
    // 
    // 
    // 
    // 
    // 
    QGauss<dim>  quadrature_formula(2);
    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();
    double nu=0.3,mu=100000000000;
    double lambda=2*mu*nu/(1-2*nu);
    //	double E=mu*2*(1+nu); 
    double temperature_epsilon[6]={0.000001,0.000001,0.000001,0,0,0};//本征应变。
    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
                                                   endc = dof_handler.end();                                               
    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
    //
    FullMatrix<double> C(6,6); 
    C=0;
    for (int i=0;i<3;i++){
    	C(i,i)=2*mu+lambda;
	}
    for (int i=3;i<6;i++){
    	C(i,i)=mu;
	}
	C(0,1)=lambda;
	C(0,2)=lambda;
	C(1,0)=lambda;
	C(1,2)=lambda;
	C(2,1)=lambda;
	C(2,0)=lambda;  
	// define elactic tensor
	//std::vector<Point<dim> > points=fe_values.get_quadrature_points();//此过程获得高斯积分点。
        int a=dof_handler.n_dofs()/3;
	Vector<double>  dof_solution_strain_xx (a);// define nodle strain_xx
	dof_solution_strain_xx=0;
	Vector<double>  dof_solution_stress_xx (a);// define nodle stress_xx
	dof_solution_stress_xx=0; 
	Vector<double>  dof_solution_strain_yy (a);// define nodle strain_yy
	dof_solution_strain_yy=0;	
	Vector<double>  dof_solution_stress_yy (a);// define nodle stress_yy
	dof_solution_stress_yy=0; 
	Vector<double>  dof_solution_strain_zz (a);// define nodle strain_zz
	dof_solution_strain_zz=0;
	Vector<double>  dof_solution_stress_zz (a);// define nodle stress_zz
	dof_solution_stress_zz=0; 
	Vector<double>  dof_solution_stress_yz (a);// define nodle stress_zz
	dof_solution_stress_yz=0; 
	Vector<double>  dof_solution_stress_xz (a);// define nodle stress_zz
	dof_solution_stress_xz=0; 
	Vector<double>  dof_solution_stress_xy (a);// define nodle stress_zz
	dof_solution_stress_xy=0; 
	Vector<double>      nodle_displace(dofs_per_cell);
	FullMatrix<double>  quadrature_shape_value(n_q_points,8);//定义高斯点的形函数值。 
	FullMatrix<double>  quadrature_shape_value_transpose(8,n_q_points);//定义高斯点的形函数值的转置。 
	FullMatrix<double>  quadrature_strain(n_q_points,6);//定义所有高斯点处的应变。 for one cell
	FullMatrix<double>  B(6,dofs_per_cell);//定义B矩阵。define B matrix
	FullMatrix<double>  quadrature_stress(n_q_points,6);
	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> 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> 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>  solution_strain_xx(8);
	Vector<double>  solution_stress_xx(8);
	Vector<double>  solution_strain_yy(8);
	Vector<double>  solution_stress_yy(8);
	Vector<double>  solution_strain_zz(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>  all_dof_solution_strain(dof_handler.n_dofs());// define nodle strain_xx
	all_dof_solution_strain=0;//define strain_x strain_y strain_z
	
	Vector<double>  all_dof_solution_stress(dof_handler.n_dofs());// define nodle strain_xx
	all_dof_solution_stress=0;//define stress_x stress_y stress_z
	
	Vector<double>  all_dof_solution_shear_stress(dof_handler.n_dofs());// define nodle strain_xx
	all_dof_solution_shear_stress=0;//define stress_x stress_y stress_z
	//creat file;
	    std::string stress_data = "stress_data-";
	    std::string strain_data = "strain_data-";
	    stress_data += ('0' + cycle);
	    strain_data += ('0' + cycle);
	    Assert (cycle < 10, ExcInternalError());
	    stress_data += ".txt";
	    strain_data += ".txt";
	    std::ofstream stress(stress_data.c_str());
	    std::ofstream strain(strain_data.c_str());
	   // strain<<"x      "<<"y       "<<"z       "<<"strain_xx       "<<"strain_yy       "<<"strain_zz       "<<"strain_yz       "<<"strain_xz       "<<"strain_xy       "<<std::endl;
	   // stress<<"x      "<<"y       "<<"z       "<<"stress_xx       "<<"stress_yy       "<<"stress_zz       "<<"stress_yz       "<<"stress_xz       "<<"stress_xy       "<<std::endl;
	 for (; cell!=endc; ++cell)
      {
	     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)=solution(local_dof_indices[i]);
          }
		quadrature_shape_value=0;//定义高斯点的形函数值。 
		quadrature_shape_value_transpose=0;//定义高斯点的形函数值的转置。 
		quadrature_strain=0;//定义所有高斯点处的应变。 
		B=0;//定义B矩阵。
		quadrature_stress=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++){
		  	   quadrature_shape_value(q_point,k)=fe_values.shape_value(3*k,q_point);	//即H矩阵。
			   quadrature_shape_value_transpose[k][q_point]=fe_values.shape_value(3*k,q_point);//transpose 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];
					  }
					  //上面求得了高斯点处的B矩阵，下面将其乘以节点位移。便可以获得高斯点处的应变。
					  for (int k=0;k<6;k++){
					  	 double sum=0;
					  	for(int j=0;j<dofs_per_cell;j++){
					  		sum+=B(k,j)*nodle_displace(j);
						  }
					  	quadrature_strain(q_point,k)=sum;
					  } 
					  if(fabs(points[q_point][0]-0.066039)<=0.000002){
					  strain<<points[q_point][0]<<"       "<<points[q_point][1]<<"       "<<points[q_point][2]<<"       "<<quadrature_strain(q_point,0)<<"       "<<quadrature_strain(q_point,1)<<"       "<<quadrature_strain(q_point,2)<<"       "<<quadrature_strain(q_point,3)<<"       "<<quadrature_strain(q_point,4)<<"       "<<quadrature_strain(q_point,5)<<std::endl;
					  }//上面求得了高斯点处的应变。下面求高斯点处的应力
					 	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){
					 	  for (int i=0;i<6;i++){
					  		double sum=0;
					  		for (int j=0;j<6;j++){
					  			sum+=C(i,j)*(quadrature_strain(q_point,j)-temperature_epsilon[j]);  
					  		}
							 quadrature_stress(q_point,i)=sum;
						    }	
						 }else{
						   for (int i=0;i<6;i++){
					  		  double sum=0;
					  		  for (int j=0;j<6;j++){
					  			sum+=C(i,j)*(quadrature_strain(q_point,j));
							  }
							 quadrature_stress(q_point,i)=sum;
						    }	
						 }
					 	if(fabs(points[q_point][0]-0.066039)<=0.000002){
					 	stress<<points[q_point][0]<<"       "<<points[q_point][1]<<"       "<<points[q_point][2]<<"       "<<quadrature_stress(q_point,0)<<"       "<<quadrature_stress(q_point,1)<<"       "<<quadrature_stress(q_point,2)<<"       "<<quadrature_stress(q_point,3)<<"       "<<quadrature_stress(q_point,4)<<"       "<<quadrature_stress(q_point,5)<<std::endl;
					 	}
		    //上面成功获得所有高斯点处的应力和应变，下面开始还原节点处的应力应变。
			//先建立节点处的应力应变还原值。
        }
      }
        /*H_T_multi_H=0;
         for (int i=0;i<8;i++){
         	for(int j=0;j<8;j++){
         		double sum=0;
         	  for (int k=0;k<n_q_points;k++){
         	  	sum+=quadrature_shape_value_transpose(i,k)*quadrature_shape_value(k,j);
			   }
			   H_T_multi_H(i,j)=sum;	
			 }
		 }
		H_T_multi_stress=0;//此过程非常重要，获得了六个应力与H的转置的乘积。 
		 for (int i=0;i<8;i++){
         	for(int j=0;j<6;j++){
         		double sum=0;
         	  for (int k=0;k<n_q_points;k++){
         	  	sum+=quadrature_shape_value_transpose(i,k)*quadrature_stress(k,j);
			   }
			   H_T_multi_stress(i,j)=sum;	
			 }
		 }
		H_T_multi_strain=0;//此过程非常重要，获得了六个应变与H的转置的乘积。 
		 for (int i=0;i<8;i++){
         	for(int j=0;j<6;j++){
         		double sum=0;
         	  for (int k=0;k<n_q_points;k++){
         	  	sum+=quadrature_shape_value_transpose(i,k)*quadrature_strain(k,j);
			   }
			   H_T_multi_strain(i,j)=sum;	
			 }
		 }
		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_stress_xx=0;
		solution_strain_yy=0;
		solution_stress_yy=0;
		solution_strain_zz=0;
		solution_stress_zz=0;
		solution_stress_yz=0;
		solution_stress_xz=0;
		solution_stress_xy=0;
		
		 
        SolverControl           solver_control (1000, 1e-10);
        SolverCG<>              solver (solver_control);
       // std::cout<<"correct"<<std::endl;
        solver.solve (H_T_multi_H, solution_strain_xx, strain_xx_rhs,PreconditionIdentity());
        solver.solve (H_T_multi_H, solution_stress_xx, stress_xx_rhs,PreconditionIdentity());
        solver.solve (H_T_multi_H, solution_strain_yy, strain_yy_rhs,PreconditionIdentity());
        
        solver.solve (H_T_multi_H, solution_stress_yy, stress_yy_rhs,PreconditionIdentity());
        solver.solve (H_T_multi_H, solution_strain_zz, strain_zz_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 (unsigned int i=0;i<8; ++i)
          { int a=local_dof_indices[3*i]/3;
          	dof_solution_strain_xx(a)=solution_strain_xx(i);//此步还原有大问题. 
          	dof_solution_stress_xx(a)=solution_stress_xx(i);
          	
          	dof_solution_strain_yy(a)=solution_strain_yy(i);//此步还原有大问题. 
          	dof_solution_stress_yy(a)=solution_stress_yy(i);
          	
          	dof_solution_strain_zz(a)=solution_strain_zz(i);//此步还原有大问题. 
          	dof_solution_stress_zz(a)=solution_stress_zz(i);
          	
          	dof_solution_stress_yz(a)=solution_stress_yz(i);
          	dof_solution_stress_xz(a)=solution_stress_xz(i);
          	dof_solution_stress_xy(a)=solution_stress_xy(i);
          }
      // std::cout<<"a="<<a;
      } 
	 // get out cell loop.
	 	  int b=dof_handler.n_dofs()/3;
	 	  for(int i=0;i<b;i++){
	 	 //recover normal_strain
	 	  all_dof_solution_strain(3*i)=dof_solution_strain_xx(i);
		  all_dof_solution_strain(3*i+1)=dof_solution_strain_yy(i);
		  all_dof_solution_strain(3*i+2)=dof_solution_strain_zz(i);
		  //recover normal_stress
		  all_dof_solution_stress(3*i)=dof_solution_stress_xx(i);
		  all_dof_solution_stress(3*i+1)=dof_solution_stress_yy(i);
		  all_dof_solution_stress(3*i+2)=dof_solution_stress_zz(i);
		 //recover shear_stress
		  all_dof_solution_shear_stress(3*i)=dof_solution_stress_yz(i);
		  all_dof_solution_shear_stress(3*i+1)=dof_solution_stress_xz(i);
		  all_dof_solution_shear_stress(3*i+2)=dof_solution_stress_xy(i);
		  
	 	  }
	 	  
	 
	 
	 

    std::string strain = "strain-";
    std::string stress = "stress-";
    std::string shear_stress = "shear_stress-";
    strain += ('0' + cycle);
    stress += ('0' + cycle);
    shear_stress+=('0' + cycle);
    Assert (cycle < 10, ExcInternalError());
    strain += ".vtk";
    stress += ".vtk";
    shear_stress+=".vtk";
    
    std::ofstream output_strain (strain.c_str());
    std::ofstream output_stress (stress.c_str());
    std::ofstream output_shear_stress (shear_stress.c_str());
    
    DataOut<dim> strain_out,stress_out,shear_stress_out;
    
    strain_out.attach_dof_handler(dof_handler);
    stress_out.attach_dof_handler(dof_handler);
    shear_stress_out.attach_dof_handler(dof_handler);
    
    std::vector<std::string> strain_name,stress_name, shear_stress_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");
    
    strain_out.add_data_vector (all_dof_solution_strain, strain_name);
    stress_out.add_data_vector (all_dof_solution_stress, stress_name);
    shear_stress_out.add_data_vector (all_dof_solution_shear_stress, shear_stress_name);
    
    strain_out.build_patches ();
    stress_out.build_patches ();
    shear_stress_out.build_patches ();
    
    strain_out.write_vtk (output_strain);
    stress_out.write_vtk (output_stress);
    shear_stress_out.write_vtk (output_shear_stress);*/
}
  /
  /
  /
  /
  /
  //下面划分网格并且控制自适应划分。 
  template <int dim>
  void ElasticProblem<dim>::run ()
  {
    for (unsigned int cycle=0; cycle<1; ++cycle)
      {
        std::cout << "Cycle " << cycle << ':' << std::endl;
        if (cycle == 0)
          {
            GridGenerator::hyper_cube (triangulation, -10,10);
            triangulation.refine_global (6);
          }
        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
    {
      Step8::ElasticProblem<3> elastic_problem_3d;
      elastic_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;
}

matlab处理数据程序
数值解部分

clear 
clc
%这里是数值解的数据
stress_data=load("stress_data-0.txt");
strain_data=load("strain_data-0.txt");
stress_data=stress_data(:,2:end);
strain_data=strain_data(:,2:end);
a_strain=unique(strain_data(:,1));
n=length(a_strain);
data_epsilon={};
for i=1:n
     p=find(a_strain(i)==strain_data(:,1));
     for j=1:length(p)
     data_epsilon{1}(i,j)=strain_data(p(j),3);%epsilon_xx
     data_epsilon{2}(i,j)=strain_data(p(j),4);%epsilon_yy
     data_epsilon{3}(i,j)=strain_data(p(j),5);%epsilon_zz
     data_epsilon{4}(i,j)=strain_data(p(j),6);%epsilon_yz
     data_epsilon{5}(i,j)=strain_data(p(j),7);%epsilon_xz
     data_epsilon{6}(i,j)=strain_data(p(j),8);%epsilon_xy
     end
end
a_stress=unique(stress_data(:,1));
data_sigma={};
for i=1:n
     p=find(a_stress(i)==stress_data(:,1));
     for j=1:length(p)
     data_sigma{1}(i,j)=stress_data(p(j),3);
     data_sigma{2}(i,j)=stress_data(p(j),4);
     data_sigma{3}(i,j)=stress_data(p(j),5);
     data_sigma{4}(i,j)=stress_data(p(j),6);
     data_sigma{5}(i,j)=stress_data(p(j),7);
     data_sigma{6}(i,j)=stress_data(p(j),8);
     end
end

解析解部分

%解析解的数据
clear
clc
nu=0.3;
mu=10^11;
E=mu*2*(1+nu);
a=1;
strain_inner=[-1/3000000*(1+nu)/(-1+nu), 0, 0;0, -1/3000000*(1+nu)/(-1+nu), 0;0, 0, -1/3000000*(1+nu)/(-1+nu)];
stress_inner=[1/750000*(nu+1)*mu/(-1+nu), 0, 0;0, 1/750000*(nu+1)*mu/(-1+nu), 0;0, 0, 1/750000*(nu+1)*mu/(-1+nu)];
x1=0.066039;
x2=[-9.93396000000000;-9.75354000000000;-9.62146000000000;-9.44104000000000;-9.30896000000000;-9.12854000000000;-8.99646000000000;-8.81604000000000;-8.68396000000000;-8.50354000000000;-8.37146000000000;-8.19104000000000;-8.05896000000000;-7.87854000000000;-7.74646000000000;-7.56604000000000;-7.43396000000000;-7.25354000000000;-7.12146000000000;-6.94104000000000;-6.80896000000000;-6.62854000000000;-6.49646000000000;-6.31604000000000;-6.18396000000000;-6.00354000000000;-5.87146000000000;-5.69104000000000;-5.55896000000000;-5.37854000000000;-5.24646000000000;-5.06604000000000;-4.93396000000000;-4.75354000000000;-4.62146000000000;-4.44104000000000;-4.30896000000000;-4.12854000000000;-3.99646000000000;-3.81604000000000;-3.68396000000000;-3.50354000000000;-3.37146000000000;-3.19104000000000;-3.05896000000000;-2.87854000000000;-2.74646000000000;-2.56604000000000;-2.43396000000000;-2.25354000000000;-2.12146000000000;-1.94104000000000;-1.80896000000000;-1.62854000000000;-1.49646000000000;-1.31604000000000;-1.18396000000000;-1.00354000000000;-0.871461000000000;-0.691039000000000;-0.558961000000000;-0.378539000000000;-0.246461000000000;-0.0660390000000000;0.0660390000000000;0.246461000000000;0.378539000000000;0.558961000000000;0.691039000000000;0.871461000000000;1.00354000000000;1.18396000000000;1.31604000000000;1.49646000000000;1.62854000000000;1.80896000000000;1.94104000000000;2.12146000000000;2.25354000000000;2.43396000000000;2.56604000000000;2.74646000000000;2.87854000000000;3.05896000000000;3.19104000000000;3.37146000000000;3.50354000000000;3.68396000000000;3.81604000000000;3.99646000000000;4.12854000000000;4.30896000000000;4.44104000000000;4.62146000000000;4.75354000000000;4.93396000000000;5.06604000000000;5.24646000000000;5.37854000000000;5.55896000000000;5.69104000000000;5.87146000000000;6.00354000000000;6.18396000000000;6.31604000000000;6.49646000000000;6.62854000000000;6.80896000000000;6.94104000000000;7.12146000000000;7.25354000000000;7.43396000000000;7.56604000000000;7.74646000000000;7.87854000000000;8.05896000000000;8.19104000000000;8.37146000000000;8.50354000000000;8.68396000000000;8.81604000000000;8.99646000000000;9.12854000000000;9.30896000000000;9.44104000000000;9.62146000000000;9.75354000000000;9.93396000000000]';%这个待得到。
x3=[-9.93396000000000;-9.75354000000000;-9.62146000000000;-9.44104000000000;-9.30896000000000;-9.12854000000000;-8.99646000000000;-8.81604000000000;-8.68396000000000;-8.50354000000000;-8.37146000000000;-8.19104000000000;-8.05896000000000;-7.87854000000000;-7.74646000000000;-7.56604000000000;-7.43396000000000;-7.25354000000000;-7.12146000000000;-6.94104000000000;-6.80896000000000;-6.62854000000000;-6.49646000000000;-6.31604000000000;-6.18396000000000;-6.00354000000000;-5.87146000000000;-5.69104000000000;-5.55896000000000;-5.37854000000000;-5.24646000000000;-5.06604000000000;-4.93396000000000;-4.75354000000000;-4.62146000000000;-4.44104000000000;-4.30896000000000;-4.12854000000000;-3.99646000000000;-3.81604000000000;-3.68396000000000;-3.50354000000000;-3.37146000000000;-3.19104000000000;-3.05896000000000;-2.87854000000000;-2.74646000000000;-2.56604000000000;-2.43396000000000;-2.25354000000000;-2.12146000000000;-1.94104000000000;-1.80896000000000;-1.62854000000000;-1.49646000000000;-1.31604000000000;-1.18396000000000;-1.00354000000000;-0.871461000000000;-0.691039000000000;-0.558961000000000;-0.378539000000000;-0.246461000000000;-0.0660390000000000;0.0660390000000000;0.246461000000000;0.378539000000000;0.558961000000000;0.691039000000000;0.871461000000000;1.00354000000000;1.18396000000000;1.31604000000000;1.49646000000000;1.62854000000000;1.80896000000000;1.94104000000000;2.12146000000000;2.25354000000000;2.43396000000000;2.56604000000000;2.74646000000000;2.87854000000000;3.05896000000000;3.19104000000000;3.37146000000000;3.50354000000000;3.68396000000000;3.81604000000000;3.99646000000000;4.12854000000000;4.30896000000000;4.44104000000000;4.62146000000000;4.75354000000000;4.93396000000000;5.06604000000000;5.24646000000000;5.37854000000000;5.55896000000000;5.69104000000000;5.87146000000000;6.00354000000000;6.18396000000000;6.31604000000000;6.49646000000000;6.62854000000000;6.80896000000000;6.94104000000000;7.12146000000000;7.25354000000000;7.43396000000000;7.56604000000000;7.74646000000000;7.87854000000000;8.05896000000000;8.19104000000000;8.37146000000000;8.50354000000000;8.68396000000000;8.81604000000000;8.99646000000000;9.12854000000000;9.30896000000000;9.44104000000000;9.62146000000000;9.75354000000000;9.93396000000000]';%这个待得到。
n=length(x3);
%下面得到的数据都是在x1=0.066039 ,x2(i)=0.066039的这根线上。
strain_xx=zeros(n,n);
strain_yy=zeros(n,n);
strain_zz=zeros(n,n);
strain_xy=zeros(n,n);
strain_xz=zeros(n,n);
strain_yz=zeros(n,n);
stress_xx=zeros(n,n);
stress_yy=zeros(n,n);
stress_zz=zeros(n,n);
stress_xy=zeros(n,n);
stress_xz=zeros(n,n);
stress_yz=zeros(n,n);
for i=1:n
    for j=1:n
   if((x1^2+x2(i)^2+x3(j)^2)<=1) 
    strain_xx(i,j)=strain_inner(1,1);
    strain_yy(i,j)=strain_inner(2,2);
    strain_zz(i,j)=strain_inner(3,3);
    strain_xy(i,j)=strain_inner(1,2);
    strain_xz(i,j)=strain_inner(1,3);
    strain_yz(i,j)=strain_inner(2,3);
    stress_xx(i,j)=stress_inner(1,1);
    stress_yy(i,j)=stress_inner(2,2);
    stress_zz(i,j)=stress_inner(3,3);
    stress_xy(i,j)=stress_inner(1,2);
    stress_xz(i,j)=stress_inner(1,3);
    stress_yz(i,j)=stress_inner(2,3);
   end
   if((x1^2+x2(i)^2+x3(j)^2)>1)
    strain_xx(i,j)=(10*a^3*((40927176184870261*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 - (129078017198436977*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 + (40927176184870261*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 - (40927176184870261*a^2*x1^2)/47223664828696452136960 + (40927176184870261*a^2*x2(i)^2)/94447329657392904273920 + (40927176184870261*a^2*x3(j)^2)/94447329657392904273920))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    strain_yy(i,j)=(10*a^3*((40927176184870261*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 - (129078017198436977*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 + (40927176184870261*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 + (40927176184870261*a^2*x1^2)/94447329657392904273920 - (40927176184870261*a^2*x2(i)^2)/47223664828696452136960 + (40927176184870261*a^2*x3(j)^2)/94447329657392904273920))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    strain_zz(i,j)=(10*a^3*((40927176184870261*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 + (40927176184870261*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 - (129078017198436977*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/94447329657392904273920 + (40927176184870261*a^2*x1^2)/94447329657392904273920 + (40927176184870261*a^2*x2(i)^2)/94447329657392904273920 - (40927176184870261*a^2*x3(j)^2)/47223664828696452136960))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    strain_xy(i,j)=-(10*a^3*x1*x2(i)*((12278152855461077*a^2)/9444732965739290427392 - (8500259669165361*a)/4722366482869645213696 + (8500259669165361*x1^2)/4722366482869645213696 + (8500259669165361*x2(i)^2)/4722366482869645213696 + (8500259669165361*x3(j)^2)/4722366482869645213696))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    strain_xz(i,j)=-(10*a^3*x1*x3(j)*((12278152855461077*a^2)/9444732965739290427392 - (8500259669165361*a)/4722366482869645213696 + (8500259669165361*x1^2)/4722366482869645213696 + (8500259669165361*x2(i)^2)/4722366482869645213696 + (8500259669165361*x3(j)^2)/4722366482869645213696))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    strain_yz(i,j)=-(10*a^3*x2(i)*x3(j)*((12278152855461077*a^2)/9444732965739290427392 - (8500259669165361*a)/4722366482869645213696 + (8500259669165361*x1^2)/4722366482869645213696 + (8500259669165361*x2(i)^2)/4722366482869645213696 + (8500259669165361*x3(j)^2)/4722366482869645213696))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_xx(i,j)=(5*a^3*((43042402838655623046875*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 - (1285123170468432173828125*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 + (43042402838655623046875*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 - (79935890986074728515625*a^2*x1^2)/1152921504606846976 + (79935890986074728515625*a^2*x2(i)^2)/2305843009213693952 + (79935890986074728515625*a^2*x3(j)^2)/2305843009213693952))/(2*((7*x1^2)/10 + (7*x2(i)^2)/10 + (7*x3(j)^2)/10)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_yy(i,j)=(5*a^3*((43042402838655623046875*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 - (1285123170468432173828125*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 + (43042402838655623046875*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 + (79935890986074728515625*a^2*x1^2)/2305843009213693952 - (79935890986074728515625*a^2*x2(i)^2)/1152921504606846976 + (79935890986074728515625*a^2*x3(j)^2)/2305843009213693952))/(2*((7*x1^2)/10 + (7*x2(i)^2)/10 + (7*x3(j)^2)/10)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_zz(i,j)=(5*a^3*((43042402838655623046875*x1^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 + (43042402838655623046875*x2(i)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 - (1285123170468432173828125*x3(j)^2*(x1^2 + x2(i)^2 + x3(j)^2 - a))/9223372036854775808 + (79935890986074728515625*a^2*x1^2)/2305843009213693952 + (79935890986074728515625*a^2*x2(i)^2)/2305843009213693952 - (79935890986074728515625*a^2*x3(j)^2)/1152921504606846976))/(2*((7*x1^2)/10 + (7*x2(i)^2)/10 + (7*x3(j)^2)/10)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_xy(i,j)=-(1000000000000*a^3*x1*x2(i)*((13*a^2)/5000000 - (9*a)/2500000 + (9*x1^2)/2500000 + (9*x2(i)^2)/2500000 + (9*x3(j)^2)/2500000))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_xz(i,j)=-(1000000000000*a^3*x1*x3(j)*((13*a^2)/5000000 - (9*a)/2500000 + (9*x1^2)/2500000 + (9*x2(i)^2)/2500000 + (9*x3(j)^2)/2500000))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
    stress_yz(i,j)=-(1000000000000*a^3*x2(i)*x3(j)*((13*a^2)/5000000 - (9*a)/2500000 + (9*x1^2)/2500000 + (9*x2(i)^2)/2500000 + (9*x3(j)^2)/2500000))/(7*(x1^2 + x2(i)^2 + x3(j)^2)*(a^2 - a + x1^2 + x2(i)^2 + x3(j)^2)^(5/2));
   end
    end
end
力学数学研究生
关注
使用Deal.II过程
monkey3的专栏
04-16 5632
deal.II 是一个有限元C++有限元库。 1、在官方网站http://dealii.org/download下载Virtual Machine Image。 2、安装vitualBox 3、在vitualBox启动虚拟机dealvm 4 、下载deal.ii 源代码，dealii-8.5.0.tar.gz。 5、展开dealii-8.5.0.tar.gz，  t
matlab做椭球面代码-CGpaper:用于CG纸的matlab软件包
05-26
Eshelby张量数值评估的最佳方案及其在MATLAB软件包中的实现，用于模拟粘性流体在缓慢流动中的运动孟梦萌，姜大志，露西X.卢计算机与地球科学，第96卷，2016年11月，第98-108页 使用此或衍生代码生成图形或其他内容的...
MMTensor 1.0：使用弹性、织物和 Eshelby 张量进行微机械计算的教程。-matlab开发
05-30
MMTensor 是一组用于微机械计算的已记录的 MATLAB 脚本。 它包含使用矩阵符号的张量进行基本计算、计算和分析织物张量以及（有限）包含建模的方法。 该软件包应被视为微机械建模的“入门”指南。 我们采取了认真的措施以有据可查，易于理解且仍然有效的形式实施这些方法。 该包主要提供基本功能。 我们实现全面或详尽无遗的目标从来都不是我们的目标。
Digimat-MF：平均场均匀化——（二）Eshelby问题和单夹杂问题
mansin_csa的专栏
09-12 5174
本文介绍Digimat-MF中的Eshelby问题和单夹杂问题。
有限元deal.ii库在wsl子系统的安装编译
最新发布
qq_45241855的博客
07-12 279
编译命令：cmake -DCMAKE_INSTALL_PREFIX=~/mycode/APP/dealii ~/mycode/dealii。命令分析：编译并确定具体的安装路径：~/mycode/APP/dealii。1、进入提前安装好的： Ubuntu22.04.3 LTS 子系统中。我们进入：APP/dealii/example/step-1。接着进行安装命令：make install。3、运行算例：make run。最后测试：make test。1、配置命令：cmake .2、编译命令：make。
推导Eshelby夹杂解析解的maple程序
py753159123456的博客
03-18 518
计算内点的弹性场 #计算内点:restart: > with(linalg): > assume(a>0): > a1:=a: > a2:=a: > a3:=a: > engstrain:=array(1..3,1..3,sparse): > for i from 1 to 3 do: > engstrain[i,i]:=10^(-6); > od: > deta:=array(1..3,1..3,sparse): > for i f
推荐开源项目：deal.II——高级有限元求解器库
gitblog_00097的博客
05-14 576
推荐开源项目：deal.II——高级有限元求解器库 项目地址:https://gitcode.com/dealii/dealii 项目介绍 deal.II是一个用C++编写的高级程序库，专注于通过自适应有限元素方法解决偏微分方程。它采用现代编程技巧，提供了一种面向对象的接口，使开发者可以轻松地处理复杂的数学数据结构和算法，为科学计算领域带来强大的工具。 项目技术分析 deal.II库的设计基于模块...
Eshelby夹杂理论的matlab程序
05-26
在matlab中，可以通过有限元方法（FEM）来实现Eshelby夹杂理论的计算。 以下是一个简单的matlab程序，可以计算一个球形夹杂在一个无限大的均匀弹性介质中的应力场和应变场： ```matlab % Eshelby夹杂理论的matlab...
Eshelby夹杂理论
05-26
Eshelby夹杂理论是一种用于研究材料中夹杂物（如裂纹、孔洞、切口等）的力学行为的理论。该理论由英国物理学家John D. Eshelby于1957年提出。 Eshelby夹杂理论的基本思想是，将夹杂物看作一个孤立的弹性体，其内部...
Matlab 函数用于模拟不连续纤维c 中流动诱导的纤维取向和机械性能_Matlab_代码_下载
07-03
通过平均场模型（Halpin-Tsai、稀释 Eshelby、Mori-Tanaka 和 Lielens/双夹杂物）的线性弹性刚度，用于单向或分布式纤维取向。 任何刚度模型的热应力和热膨胀张量。 使用经典层压板理论，其方向随厚度变化的复合板的...
deal.II链接PETSc过程记录
摄心神，致虚极，守静笃
11-23 1826
deal.II链接PETSc过程记录 Posted on 2016-09-29   |   In computational material science   |   暂无评论 2016-11-16 Update p4est与deal.II的链接注意事项： p4est安装时需要开启mpi，即configure时加上–enable-mpi选项。链接时增添：
有限元开源代码dealii_C编写
05-18
有限元开源代码，C编写的多物理场有限元dealii。目前活跃度较高。
Collected Works Of Eshelby(part 1)
02-03
固体力学的大牛人，提出了夹杂理论的框架和基本概念，以及演化微结构中的构形力，能量动量张量等光辉的思想。
dealii:deal.II有限元库的开发资料库
03-17
什么是交易II？ deal.II是一个C ++程序库，旨在使用自适应有限元解决偏微分方程的计算问题。 它使用最先进的编程技术，为您提供所需的复杂数据结构和算法的现代化接口。 对于急躁的人： 假设您已经将.tar.gz文件解压缩到目录/ path / to / dealii / sources中。 然后使用以下命令配置，编译和安装deal.II库： $ mkdir build $ cd build $ cmake -DCMAKE_INSTALL_PREFIX=/path/where/dealii/should/be/installed/to /path/to/dealii/sources $ make install (alternatively $ make -j<N> install) $ make test 要从存储库进行构建，请首先执行以下命令： $ git clone
开源有限元dealii学习——step-1
weixin_42312251的博客
05-02 3151
#简介 dealii是一款开源的求解偏微分方程的有限元软件，它优点： 1.多种单元类型 2.可以自适应网格 3.文档和范例齐全 4.与其它库有良好的接口 5.使用c++编写 #安装 dealii可以在官网上或者github上下载源码，解压后进入源文件目录安装。 mkdir build cd build cmake -DCMAKE_INSTALL_PREFIX=/path/to/install...
2.Windows系统下利用Ubuntu安装deal.II的可视化界面
lrj688的博客
12-12 683
上一篇文章介绍了如何在Windows系统下利用Ubuntu安装deal.II，其实如果对实时可视化要求不高的话，可以输出文件本地查看，不用进行可视化界面的设置，但有了可视化界面显然要方便一些。若只想进行可视化操作，可以先参考专栏上一篇文章进行Ubuntu下载以及WSL设置，成功后进行如下操作。
1.Windows系统下利用Ubuntu安装Deal.ii
lrj688的博客
11-23 2258
Windows系统目前支持Visio Studio直接安装Deal.ii，但官网也提到正处于试验阶段，另一方面有些库只支持Linux系统，所以还是要考虑在Windows中开辟出一个Linux来。可以选择虚拟机，但Windows目前也可以使用WSL（Windows Subsystem for Linux）建立Linux子系统，非常方便。本文将讲解如何在WSL中使用Ubuntu安装dealii。
Deal.ii安装
ppzzgg666的博客
02-27 458
需要安装hdf5支持 sudo apt-get install libhdf5-serial-dev（此路不通，还是自己下载编译hdf5比较方便https://www.jianshu.com/p/d0db51925797，注意要enable-parallel），然后cmake的时候-DHDF5_DIR=/path/to/hdf5。安装教程：https://www.bilibili.com/video/BV12x411d7Uv?需要使用测试功能，还需要安装numdiff，安装过程看官网就可以，非常简单。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值