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

最新推荐文章于 2025-07-06 09:43:28 发布
原创 最新推荐文章于 2025-07-06 09:43:28 发布 · 644 阅读 · 1 ·
CC 4.0 BY-SA版权
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

本文介绍了一个三维弹性问题的数值求解方法，采用deal-ii库实现有限元方法，并通过自适应网格细化提高求解精度。文章详细展示了从网格生成到后处理的全过程，包括应力应变的计算与输出。

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;
最低0.47元/天 解锁文章
使用Deal.II过程
monkey3的专栏
04-16 6070
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
开源有限元dealii学习——step-1
weixin_42312251的博客
05-02 4201
#简介 dealii是一款开源的求解偏微分方程的有限元软件，它优点： 1.多种单元类型 2.可以自适应网格 3.文档和范例齐全 4.与其它库有良好的接口 5.使用c++编写 #安装 dealii可以在官网上或者github上下载源码，解压后进入源文件目录安装。 mkdir build cd build cmake -DCMAKE_INSTALL_PREFIX=/path/to/install...
MMTensor 1.0：使用弹性、织物和 Eshelby 张量进行微机械计算的教程。-matlab开发
05-30
MMTensor 是一组用于微机械计算的已记录的 MATLAB 脚本。 它包含使用矩阵符号的张量进行基本计算、计算和分析织物张量以及（有限）包含建模的方法。 该软件包应被视为微机械建模的“入门”指南。 我们采取了认真的措施以有据可查，易于理解且仍然有效的形式实施这些方法。 该包主要提供基本功能。 我们实现全面或详尽无遗的目标从来都不是我们的目标。
Digimat-MF：平均场均匀化——（二）Eshelby问题和单夹杂问题
mansin_csa的专栏
09-12 6183
本文介绍Digimat-MF中的Eshelby问题和单夹杂问题。
推导Eshelby夹杂解析解的maple程序
py753159123456的博客
03-18 752
计算内点的弹性场 #计算内点: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链接PETSc过程记录
摄心神，致虚极，守静笃
11-23 2247
deal.II链接PETSc过程记录 Posted on 2016-09-29   |   In computational material science   |   暂无评论 2016-11-16 Update p4est与deal.II的链接注意事项： p4est安装时需要开启mpi，即configure时加上–enable-mpi选项。链接时增添：
Eshelby夹杂理论的matlab程序
05-26
在matlab中，可以通过有限元方法（FEM）来实现Eshelby夹杂理论的计算。 以下是一个简单的matlab程序，可以计算一个球形夹杂在一个无限大的均匀弹性介质中的应力场和应变场： ```matlab % Eshelby夹杂理论的matlab...
Eshelby夹杂理论
05-26
Eshelby夹杂理论是一种用于研究材料中夹杂物（如裂纹、孔洞、切口等）的力学行为的理论。该理论由英国物理学家John D. Eshelby于1957年提出。 Eshelby夹杂理论的基本思想是，将夹杂物看作一个孤立的弹性体，其内部...
基于MATLAB的椭球面与Eshelby张量模拟工具包
因此，借助MEX技术将最耗时的部分用C语言实现，既保留了MATLAB强大的数据可视化与交互优势，又获得了接近原生C程序的执行效率。 此外，从标签信息可以看出，该软件包涵盖了多个跨学科知识点：包括MATLAB编程环境的...
2.Windows系统下利用Ubuntu安装deal.II的可视化界面
lrj688的博客
12-12 1002
上一篇文章介绍了如何在Windows系统下利用Ubuntu安装deal.II，其实如果对实时可视化要求不高的话，可以输出文件本地查看，不用进行可视化界面的设置，但有了可视化界面显然要方便一些。若只想进行可视化操作，可以先参考专栏上一篇文章进行Ubuntu下载以及WSL设置，成功后进行如下操作。
有限元开源代码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
1.Windows系统下利用Ubuntu安装Deal.ii
lrj688的博客
11-23 3416
Windows系统目前支持Visio Studio直接安装Deal.ii，但官网也提到正处于试验阶段，另一方面有些库只支持Linux系统，所以还是要考虑在Windows中开辟出一个Linux来。可以选择虚拟机，但Windows目前也可以使用WSL（Windows Subsystem for Linux）建立Linux子系统，非常方便。本文将讲解如何在WSL中使用Ubuntu安装dealii。
Deal.ii安装
ppzzgg666的博客
02-27 1148
需要安装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，安装过程看官网就可以，非常简单。
有限元deal.ii库在wsl子系统的安装编译
qq_45241855的博客
07-12 906
编译命令：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。
dealii-step-3
weixin_42312251的博客
05-02 583
算例3介绍了泊松方程在dealii中的求解 泊松方程为边界为0，右端项非零,可以转化为AU=F： 使用类Step3开始程序的运行： class Step3 { public: Step3(); void run(); 使用一些私有对象完成对应的功能： private: void make_grid(); void setup_system(); void assemble_s...
DEALII适配器在OpenFOAM7与precice2.1.1中的应用
最新发布
weixin_42466723的博客
07-06 709
DEALII（Distributed and Efficient Algorithms for Finite Element Analysis and Computational Fluid Dynamics）是一个高度模块化且灵活的C++库，主要用于有限元分析和计算流体动力学。它广泛应用于物理、工程和其他科学领域。DEALII的一个重要特征是其在网格处理和离散化方面的能力，它支持多种类型的单元和几何，提供了一套完备的工具集以构建和模拟复杂的物理现象。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值