Matlab编写comsol模型——创建不相交圆柱体

matlab创建Comsol中的圆柱体较为简单，但是在有些模型中需要创建多个随机分布的且不相交的圆柱体，这需要在matlab中进行相交检测的操作。

圆柱体属于凸多面体，对于凸多面体的检测通常可用GJK算法。那么我们可以先建立Comsol圆柱体模型的外包络多面体，然后在判断这些多面体是否相交。

创建与Comsol模型对应的多面体

利用Comsol创建不同取向的圆柱体的代码可参考以下方式：

model.component('comp1').geom('geom1').create('Cylinder_1', 'Cylinder');
model.component('comp1').geom('geom1').feature('Cylinder_1').set('pos', [0 0 0]);
model.component('comp1').geom('geom1').feature('Cylinder_1').set('r', r);
model.component('comp1').geom('geom1').feature('Cylinder_1').set('h', h);
geom1.feature('Cylinder_1').set('ax3'，['{',num2str(n(1)),',',num2str(n(2)),',',num2str(n(3)),'}']);

先在Comsol中创建一个圆柱体，在对他进行平移、拉伸和旋转等操作，但是要想在Matlab中使用GJK算法判断多个圆柱体之间是否相交，首先要创建圆柱体的外包络多面体，且必须与Comsol中平移、拉伸和旋转后的圆柱体对应。Comsol中平移的起点是圆柱体轴线的起点位置，高度是圆柱体轴线的长度，方向是圆柱体轴线的法向量。因此在利用matlab的cylinder函数建立圆柱体外包络几何体的信息时必须将圆柱体的方向旋转到Comsol中‘ax3’指定的方向。cylinder函数的初始方向为z轴方向，即[0 0 1]；需将其旋转至‘ax3’指定的方向，在Matlab中可利用rotationVectorToMatrix函数得到旋转矩阵，代码如下：

代码如下：

[Xr,Yr,Zr]=Rotate_Cylinder(X0*r_z*ErrorCoof,Y0*r_z*ErrorCoof,Z0*H_z*ErrorCoof,NormalVec_Cylinder);
Xn=Xr+x;Yn=Yr+y;Zn=Zr+z;
function [X,Y,Z]=Rotate_Cylinder(X0,Y0,Z0,ve)
    ve=ve/sqrt(ve*ve');
    vr=cross(ve,[0 0 1]);vr=vr/sqrt(vr*vr');vr=vr*acos([0 0 1]*ve');
    rm=rotationVectorToMatrix(vr);
    SurfNum=size(X0,1);
    P=[reshape(X0,SurfNum*2,1) reshape(Y0,SurfNum*2,1) reshape(Z0,SurfNum*2,1)]';
    P=rm*P;P=P';
    vD=[1,-ve(:,1)/ve(:,2),0];vD=vD/(sqrt(vD*vD'));
    vF=cross(vD,ve);vF=vF/(sqrt(vF*vF'));
    vx=rm*[1 0 0]';vx=vx';
    vee=ve*acos(vF*vx');
    rm=rotationVectorToMatrix(vee);
    P=rm*P';P=P';
    X=reshape(P(:,1),SurfNum,2);
    Y=reshape(P(:,2),SurfNum,2);
    Z=reshape(P(:,3),SurfNum,2);
end

判断相交

例如[Xn,Yn,Zn]表示新的圆柱体的外包络多面体，S2表示已有圆柱体的外包络多面体，判断已有圆柱体与其它圆柱体是否相交的代码如下：

FLAG=0;
    S1.XData=Xn;S1.YData=Yn;S1.ZData=Zn;
    for n_index=1:All_Cylinder_Index-1
        S2.XData=All_Cylinder(:,1:2,n_index);
        S2.YData=All_Cylinder(:,3:4,n_index);
        S2.ZData=All_Cylinder(:,5:6,n_index);
        if GJK(S1,S2,6)
           FLAG=1;
           break
        end
    end
    if FLAG==1
        '相交'
    end