MoM(二)

一、计算单元面积

1、原理

利用向量的思想，

其中，与分别为向量AB与AC在空间直角坐标系下的坐标表达，即：向量构成三角形面积等于向量构成平行四边形面积的一半

即三角型的面积是四边形面积的一半

2、重难点

写代码的时候，发现函数不能返回数组，所以在这记一下。

1、使用指针

int* fun(int x1[3], int x2[3])
{
    int* x3 = new int[3];

    x3[0] = x1[0] + x2[0];
    x3[1] = x1[1] + x2[1] ;
    x3[2] = x1[2] + x2[2] ;

    return x3;
}

int main()
{
   int x1[3] = { 1,2,3 };
   int x2[3] = { 5,6,7 };
   int* x3 = new int[3];
   x3 = fun(x1, x2);
   for (int i = 0; i < 3; i++) {
       cout << x3[i] << '\t';
   }
   delete[] x3;   //必不可少
}
  

2、使用 static 

int* fun(int x1[3], int x2[3])
{
    //int* x3 = new int[3];
    static int x3[3];
    x3[0] = x1[0] + x2[0];
    x3[1] = x1[1] + x2[1] ;
    x3[2] = x1[2] + x2[2] ;

    return x3;
}

int main()
{
   int x1[3] = { 1,2,3 };
   int x2[3] = { 5,6,7 };
   //int* x3 = new int[3];
   int *x3;
   x3 = fun(x1, x2);
   for (int i = 0; i < 3; i++) {
       cout << x3[i] << '\t';
   }
}

3、代码

/*向量叉乘运算*/
double* cross(double x1[3], double x2[3])
{
	double* x3 = new double[3];
	x3[0] = x1[1] * x2[2] - x1[2] * x2[1];
	x3[1] = x1[2] * x2[0] - x1[0] * x2[2];
	x3[2] = x1[0] * x2[1] - x1[1] * x2[0];
	return x3;
}
/*计算每个三角片面的面积*/
void area()
{
	vec area(elem_num);
	mat center(elem_num, 3);
	double vector1[3], vector2[3];   //两向量
	double * vector3 = new double[3]; //叉乘结果
	for (int i = 0; i < elem_num; i++) {
		for (int j = 0; j < 3; j++) {
			vector1[j] = element_list.slice(i)(0, j) - element_list.slice(i)(1, j);
			vector2[j] = element_list.slice(i)(2, j) - element_list.slice(i)(1, j);
		}
		vector3 = cross(vector1, vector2); 
		area(i) = 0.5 * sqrt(pow(vector3[0], 2) + pow(vector3[1], 2) + pow(vector3[2], 2));
		/*center(i, 0) = (1.0 / 3.0) * (element_list.slice(i)(0, 0) + element_list.slice(i)(1, 0) + element_list.slice(i)(2, 0));
		center(i, 1) = (1.0 / 3.0) * (element_list.slice(i)(0, 1) + element_list.slice(i)(1, 1) + element_list.slice(i)(2, 1));
		center(i, 2) = (1.0 / 3.0) * (element_list.slice(i)(0, 2) + element_list.slice(i)(1, 2) + element_list.slice(i)(2, 2));*/
		for (size_t k = 0; k < 3; k++) {
			center(i, k) = (1.0 / 3.0) * (element_list.slice(i)(0, k) + element_list.slice(i)(1, k) + element_list.slice(i)(2, k));
		}
	}
	delete[] vector3;
	Area = area;
	Center = center;
}

二、求公共边的长度

/*寻找组成公共边的两个点 以及 非公共边的点*/
void commonEdge() {
	mat commonEdge_p(rwg_num, 2);
	mat commonEdge_n(rwg_num, 2);
	vec vertex_pos(rwg_num);
	vec vertex_nega(rwg_num);
	for (int i = 0; i < rwg_num; i++) {
		int t = 0;
		for (int m = 0; m < 3; m++) {
			for (int n = 0; n < 3; n++) {
				//slice 选层，即选 element    (basis_list 里面是配对了的两个 element)
			// (m,0)  m:组成element的第m个点   0:x坐标   1：y坐标  2：z坐标
				if (element_list.slice(basis_list(i, 0))(m, 0) == element_list.slice(basis_list(i, 1))(n, 0) &&
					element_list.slice(basis_list(i, 0))(m, 1) == element_list.slice(basis_list(i, 1))(n, 1) &&
					element_list.slice(basis_list(i, 0))(m, 2) == element_list.slice(basis_list(i, 1))(n, 2))
				{
					commonEdge_p(i, t) = m;
					commonEdge_n(i, t) = n;
					t++;
				}
				//if (t > 2) break;
			}
			if (t > 2) {
				printf("相同点数大于2，退出");
				break;
			}
		}
		for (int k = 0; k < 3; k++) {
			if (k != commonEdge_p(i,0) && k != commonEdge_p(i,1)) {
				vertex_pos(i) = k;
			}
			if (k != commonEdge_n(i, 0) && k != commonEdge_n(i, 1)) {
				vertex_nega(i) = k;
			}
		}

		//if (t > 2) break;
	}
	CommonEdgeP = commonEdge_p;
	CommonEdgeN = commonEdge_n;
	Vertex_pos = vertex_pos;
	Vertex_naga = vertex_nega;
}
void edgelength()
{
	commonEdge();
	vec edgeLength(rwg_num);
	rowvec vector1(1, 3), vector2(1, 3), vector3(1, 3);
	for (int i = 0; i < rwg_num; i++) {
		// 取出 被定义为正的element 的编号
		int NoPos = basis_list(i, 0);
		vector1 = element_list.slice(NoPos).row(CommonEdgeP(i, 0)) - element_list.slice(NoPos).row(CommonEdgeP(i, 1));
		edgeLength(i) = sqrt(as_scalar(vector1 * strans(vector1)));
	}
	EdgeLength = edgeLength;
}