MoM矩量法（一）：剖分数据的处理

最新推荐文章于 2023-03-02 09:53:07 发布

置顶 clamtian

最新推荐文章于 2023-03-02 09:53:07 发布

阅读量3.1k

点赞数 7

分类专栏： MoM矩量法计算电磁学文章标签：电磁学

本文链接：https://blog.csdn.net/luckytsw/article/details/107545684

版权

计算电磁学同时被 2 个专栏收录

4 篇文章 3 订阅

订阅专栏

MoM矩量法

3 篇文章 5 订阅

订阅专栏

MoM矩量法（一）：剖分数据的处理

写在开头：沉寂了一年多，终于又要开始写博客了。这一年多来经历了很多事情，我也成长了不少，从Java到机器学习再到计算电磁学，方向的转变还是挺大的。虽然方向在变，但之前写代码留下的功底作用还是非常大的。好了，废话不多说，开始进入主题吧。
矩量法的计算理论还有公式推导很多专门的书上都有，这里不再赘述。需要注意的是，要看懂书上的推导过程还是需要比较深的物理和数学功底的，所以一开始看的比较吃力的话，不要担心，慢慢来。
这里采用的公式来自于聂在平主编的《目标与环境电磁散射特性建模——理论、方法与实现（应用篇）》，所有数据均采用国际标准单位。公式的推导就不再多说了，直接开始coding吧！

完整代码已贴在文章最后。文章和代码思路借鉴了https://blog.csdn.net/u014411646/article/details/99689460，谢谢原作者！

一、剖分数据

初始剖分数据包括两个文件，分别是包含点元素信息的文件pdata.txt和包含三角单元信息的文件tdata.txt，先一起来看一看初始数据。

1.1点元素数据文件pdata.txt

在这里插入图片描述
横向为坐标值，纵向为点元素序号。例如，上图中显示第4个点的坐标为(-3.7033783e-02, 6.4144393e-02, 4.4386259e-01)。所以pdata.txt的行个数就代表点元素的个数。

1.2三角单元数据文件tdata.txt

在这里插入图片描述
横向为点元素的序号，纵向为三角单元序号。例如，从上图中可以知道第3个三角单元的三个顶点分别为1,4和5。则tdata.txt的行个数即为三角单元的个数。

上面两组数据的意义是非常明显的，例如，从第二组数据中我们可以得到第2个三角单元的三个顶点分别是1,3和4，而从第一组数据又可以得到第1号点的坐标为（0.0000000e+00 0.0000000e+00 4.5000000e-01），第3号点的坐标为（3.7033783e-02 6.4144393e-02 4.4386259e-01），第4号点的坐标为（-3.7033783e-02 6.4144393e-02 4.4386259e-01）。所以，第2个三角单元三个顶点坐标分别为（0.0000000e+00 0.0000000e+00 4.5000000e-01）、（3.7033783e-02 6.4144393e-02 4.4386259e-01）和（-3.7033783e-02 6.4144393e-02 4.4386259e-01）。

二、初始数据的处理

初始数据的处理其实有很多种方法，但核心都是为了之后计算阻抗矩和激励矩阵的简便。
接下来对一些核心代码进行讲解：

mat Edge_list;              //[端点1， 端点2， 对点1， 对点2， 左三角单元， 右三角单元]
mat TriCenter;
vec EdgeLength;
vec TriArea;
int TriTotal, EdgeTotal;
mat pdata, tdata;

TriTotal为三角单元的个数，EdgeTotal为边的个数；
Edge_list是一个维度为（EdgeTotal，6）的二维矩阵，其中储存的是剖分数据中边的信息，行代表第几条边，列从左到右依次代表[第一个端点，第二个端点，与边相对的第一个端点，与边相对的第二个端点，左边的三角单元，右边的三角单元]；
TriCenter是（TriTotal，3）的二维矩阵，储存的是三角单元的重心坐标；
EdgeLength是长度为EdgeTotal向量，储存的是对应边的长度；
TriArea是长度为TriTotal向量，储存的是对应三角单元的面积；
pdata, tdata为原始的pdata.txt和tdata.txt初始数据。

下面贴上代码：

/***********************************************************/
/*         该函数主要用来对输入的初始数据进行处理            */
/*                                                         */
/*                                                         */
/***********************************************************/
/*                                                         */
/*                 Author：Tianshaowenwen                  */
/*                                                         */
/***********************************************************/
#include "rwg1.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "armadillo"

using namespace arma;
using namespace std;


extern string p_filename, t_filename;


mat Edge_list;              //[端点1， 端点2， 对点1， 对点2， 左三角单元， 右三角单元]
mat TriCenter;
vec EdgeLength;
vec TriArea;
int TriTotal, EdgeTotal;
mat pdata, tdata;

void rwg1_main(){
	pdata.load(p_filename, raw_ascii);
	tdata.load(t_filename, raw_ascii);
	TriTotal = tdata.n_rows;
	edgeNumber();
	edgeElement();
	area();
	Edge_list.save("Edge_list.txt", raw_ascii);
	TriCenter.save("TriCenter.txt", raw_ascii);
	EdgeLength.save("EdgeLength.txt", raw_ascii);
	TriArea.save("TriArea.txt", raw_ascii);
	cout << "\n" << "TriTotal= " << TriTotal << endl;
	cout << "\n" << "EdgeTotal= " << EdgeTotal << endl;
}
/***********************************************************/
/*              Function：void edgeNumber()                */
/***********************************************************/
/* 1.统计边的个数                                          */
/***********************************************************/
void edgeNumber(){
	int edgetotal = 0;
	for (int i = 0;i < TriTotal;i++) {
		int N[3] = { 0 }, M[3] = { 0 };
		N[0] = tdata(i, 0) - 1;
		N[1] = tdata(i, 1) - 1;
		N[2] = tdata(i, 2) - 1;
		for (int j = i + 1;j < TriTotal;j++) {
			M[0] = tdata(j, 0) - 1;
			M[1] = tdata(j, 1) - 1;
			M[2] = tdata(j, 2) - 1;
			int sameNode = 0;
			for(int i1 = 0;i1 < 3;i1++) {
				for (int i2 = 0;i2 < 3;i2++) {
					if (N[i1] == M[i2]) {
						sameNode += 1;
					}
				}
			}
			if (sameNode == 2) {
				edgetotal += 1;
			}
		}
	}
	EdgeTotal = edgetotal;
}

/***********************************************************/
/*             Function：void edgeelement()                */
/***********************************************************/
/* 1.绑定edge两端节点和相对节点和左右三角单元                */
/* 2.求出edge长度                                           */
/*                                                         */
/***********************************************************/
void edgeElement(){
	vec edgeLength(EdgeTotal);
	mat edge_list(EdgeTotal, 6);
	int edgetotal = 0;
	for (int i = 0;i < TriTotal;i++) {
		int N[3] = { 0 }, M[3] = { 0 };
		N[0] = tdata(i, 0) - 1;
		N[1] = tdata(i, 1) - 1;
		N[2] = tdata(i, 2) - 1;
		for (int j = i + 1;j < TriTotal;j++) {
			M[0] = tdata(j, 0) - 1;
			M[1] = tdata(j, 1) - 1;
			M[2] = tdata(j, 2) - 1;
			int sameNode = 0;
			int same[2] = { 0 };
			for (int i1 = 0;i1 < 3;i1++) {
				for (int i2 = 0;i2 < 3;i2++) {
					if (N[i1] == M[i2]) {
						same[sameNode] = N[i1];
						sameNode += 1;
					}
				}
			}
			if (sameNode == 2) {
				edge_list(edgetotal, 0) = same[0];
				edge_list(edgetotal, 1) = same[1];
				for (int j = 0;j < 3;j++) {
					if (N[j] != same[0] && N[j] != same[1]) {
						edge_list(edgetotal, 2) = N[j];
					}
					if (M[j] != same[0] && M[j] != same[1]) {
						edge_list(edgetotal, 3) = M[j];
					}
				}
				edge_list(edgetotal, 4) = i;
				edge_list(edgetotal, 5) = j;

				double len = getDist(pdata.row(same[0]), pdata.row(same[1]));
				edgeLength(edgetotal) = len;
				edgetotal += 1;
			}
		}
	}
	Edge_list = edge_list;
	EdgeLength = edgeLength;
}
/***********************************************************/
/*                 Function：void area()                  */
/***********************************************************/
/* 1.计算三角单元的面积                                     */
/* 2.计算三角单元的重心                                     */
/***********************************************************/
void area() {
	vec triArea(TriTotal);
	mat triCenter(TriTotal, 3);
	float test = 0, test1 = 0;
	vec vector1(3), vector2(3);
	int N[3] = { 0 };
	for (int ia = 0;ia < TriTotal;ia++)
	{
		for (int i = 0;i < 3;i++)
			N[i] = tdata(ia, i) - 1;
		for (int j = 0;j < 3;j++)
		{
			vector1(j) = pdata(N[0], j) - pdata(N[1], j);
			vector2(j) = pdata(N[2], j) - pdata(N[1], j);
		}
		vec vector3 = cross(vector1, vector2);
		triArea(ia) = 0.5 * sqrt(vector3(0) * vector3(0) + vector3(1) * vector3(1) + vector3(2) * vector3(2));
		triCenter(ia, 0) = (1.0 / 3.0) * (pdata(N[0], 0) + pdata(N[1], 0) + pdata(N[2], 0));
		triCenter(ia, 1) = (1.0 / 3.0) * (pdata(N[0], 1) + pdata(N[1], 1) + pdata(N[2], 1));
		triCenter(ia, 2) = (1.0 / 3.0) * (pdata(N[0], 2) + pdata(N[1], 2) + pdata(N[2], 2));
	}
	TriArea = triArea;
	TriCenter = triCenter;
}
/***********************************************************/
/*                 Function：double getDist()              */
/***********************************************************/
/* 1.求两个向量的距离                                       */
/***********************************************************/
double getDist(rowvec r1, rowvec r2) {
	double dx = r1(0) - r2(0);
	double dy = r1(1) - r2(1);
	double dz = r1(2) - r2(2);
	double len = sqrt(dx * dx + dy * dy + dz * dz);
	return len;
}
/***********************************************************/
/*                 Function：rowvec cross()                */
/***********************************************************/
/* 1.求两个向量的差乘                                       */
/***********************************************************/
rowvec cross(rowvec r1,rowvec r2){
	rowvec res(3);
	res(0) = r1(1) * r2(2) - r1(2) * r2(1);
	res(1) = r1(2) * r2(0) - r1(0) * r2(2);
	res(2) = r1(0) * r2(1) - r1(1) * r2(0);
	return res;
}
/***********************************************************/
/*                 Function：double getDot()               */
/***********************************************************/
/* 1.求两个向量的点乘                                       */
/***********************************************************/
double getDot(rowvec r1, rowvec r2) {
	return r1(0) * r2(0) + r1(1) * r2(1) + r1(2) * r2(2);
}
/***********************************************************/
/*                 Function：vec getUnit()                 */
/***********************************************************/
/* 1.求单位向量                                            */
/***********************************************************/
rowvec getUnit(rowvec r1) {
	rowvec res(3);
	res(0) = 0.0;
	res(1) = 0.0;
	res(2) = 0.0;
	double len = sqrt(r1(0) * r1(0) + r1(1) * r1(1) + r1(2) * r1(2));
	if (len != 0.0) {
		res(0) = r1(0) / len;
		res(1) = r1(1) / len;
		res(2) = r1(2) / len;
	}
	return res;
}

今天写的这第一部分其实就是对初始数据的一些处理和需要用到的工具函数的创建。矩量法的重头其实还是阻抗矩阵的填充和矩阵方程的求解，关于阻抗矩阵的填充我会放在下次讲。

clamtian

关注

7
点赞
踩
19

收藏

觉得还不错? 一键收藏
1
评论
MoM矩量法（一）：剖分数据的处理

MoM矩量法（一）：剖分数据的处理写在开头：沉寂了一年多，终于又要开始写博客了。这一年多来经历了很多事情，我也成长了不少，从Java到机器学习再到计算电磁学，方向的转变还是挺大的。虽然方向在变，但之前写代码留下的功底作用还是非常大的。好了，废话不多说，开始进入主题吧。矩量法的计算理论还有公式推导很多专门的书上都有，这里不再赘述。需要注意的是，要看懂书上的推导过程还是需要比较深的物理和数学功底的，所以一开始看的比较吃力的话，不要担心，慢慢来。这里采用的公式来自于聂在平主编的《目标与环境电磁散射特性建模
复制链接

扫一扫