Palabos_permeability解读

生成多孔介质几何体

\examples\tutorial\permeability 中提供了一个matlab函数 createDAT.m，以及一系列切面图。可以创建 matlab 文件 create.m 调用 createDAT.m 函数，在 create.m 中给定图像数、图像路径、图像基本名、三维结构名，将 \examples\tutorial\permeability\twoSpheres 中的切面图合成三维结构，本次模拟中，相关参数为 (48， ‘twoSpheres/’, ‘twoSpheres’, ‘twoSpheres.dat’ )。 可以生成一根沿x方向的管子，流体从中间流过，文件保存在 twoSphere.dat 中。同时对不同的区域赋予不同的数值：将固体区域设定为2，液体设定为0，固液之间的边界设定为1。

相关头文件

#include <cmath>
#include <cstdlib>
#include <vector>

#include "palabos3D.h"
#include "palabos3D.hh"

using namespace plb;

typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor

读取几何体

将 twoSpheres.dat 中的文件保存到多块结构 MultiScalarField3D 中。
其中利用了 generateMultiScalarField<int>(geometry, sliceBox)函数。根据 multiBlockGenerator3D.hh 文件中的说明：Generate a new multi-scalar-field with same distribution as original block, but intersected with a given domain. No data is being copied, the new scalar-field is initialized to zero values. The new boundingBox is either equal to the original one, or equal to the intersection, depending on the parameter crop. 但经过实际测试，得到的结果数据结构与 geometry 相同，但大小与 sliceBox 相同。

void readGeometry(std::string fNameIn, std::string fNameOut, MultiScalarField3D<int> &geometry)
{
   
    const plint nx = geometry.getNx();
    const plint ny = geometry.getNy();
    const plint nz = geometry.getNz();

    Box3D sliceBox(0, 0, 0, ny - 1, 0, nz - 1);  //在入口处生成一个截面
   // 创建一个与 geometry 类型相同，与 sliceBox 大小相同的新 multi-scalar-field,命名为 slice
    std::unique_ptr<MultiScalarField3D<int> > slice =
        generateMultiScalarField<int>(geometry, sliceBox);  

    plb_ifstream geometryFile(fNameIn.c_str());  //读取 twoSpheres.dat 的数据
    for (plint iX = 0; iX < nx - 1; ++iX) {
     // 利用循环提取 geometryFile 中的数据
        if (!geometryFile.is_open()) {
     // 检验文件是否有效
            pcout << "Error: could not open geometry file " << fNameIn << std::endl;
            exit(EXIT_FAILURE);
        }
        // 将 1*Ny*Nz 大小的 geometry 数据导入 slice，并利用下面的 copy 函数导入 geometry 中。
        geometryFile >> *slice;  
        copy(*slice, slice->getBoundingBox(), geometry, Box3D(iX, iX,