tmp update

README.md

@@ -0,0 +1,12 @@
+### Compile
+
+```shell
+g++ main.cpp --std=c++11 -o a.exe
+```
+
+### Run
+
+```shell
+./a.exe
+```
+

cube.geo

@@ -0,0 +1,24 @@
+//+
+Point(1) = {40, 40, -20, 5};
+//+
+Point(2) = {60, 40, -20, 5};
+//+
+Point(3) = {60, 60, -20, 5};
+//+
+Point(4) = {40, 60, -20, 5};
+//+
+Line(1) = {4, 1};
+//+
+Line(2) = {1, 2};
+//+
+Line(3) = {2, 3};
+//+
+Line(4) = {3, 4};
+//+
+Curve Loop(1) = {1, 2, 3, 4};
+//+
+Plane Surface(1) = {1};
+//+
+Extrude {0, 0, -40} {
+  Surface{1}; 
+}

main.cpp

@@ -0,0 +1,42 @@
+#include "tetra_vz.h"
+#include "iostream"
+
+int main(int argc, char *argv[])
+{
+	// Declare a tetrahedron
+	point3d p[4];
+	p[0].x = 20; p[0].y = 50; p[0].z = -5;
+	p[1].x = 80; p[1].y = 20; p[1].z = -10;
+	p[2].x = 50; p[2].y = 80; p[2].z = -15;
+	p[3].x = 50; p[3].y = 50; p[3].z = -40;
+
+	tetrahedron t;
+	t.rho = 1.0;
+	for (int i = 0; i < 4; i++)
+	{
+		t.vec_ptr[i] = &p[i];
+	}
+	t.initialize_tensors();
+
+	// Declare the observation stations
+	double xmin = 0, dx = 2.5; // 41
+	double ymin = 0, dy = 2.5; // 41
+	point3d obs[1681]; // 41*41
+	for (int i = 0; i < 41; i++)
+	{
+		for (int j = 0; j < 41; j++)
+		{
+			obs[i*41+j].x = xmin + dx*j;
+			obs[i*41+j].y = ymin + dy*i;
+			obs[i*41+j].z = 0.0;
+		}
+	}
+
+	// Calculate the gravity values
+	for (int i = 0; i < 1681; i++)
+	{
+		std::cout << obs[i].x << " " << obs[i].y << " ";
+		std::cout << tetra_vz(&t, &obs[i]) << std::endl;
+	}
+	return 0;
+}

main_msh.cpp

@@ -0,0 +1,138 @@
+#include "tetra_vz.h"
+#include "string"
+#include "sstream"
+#include "iostream"
+#include "fstream"
+
+std::stringstream str2ss(std::string in)
+{
+	std::stringstream tmp_ss;
+	tmp_ss.clear();
+	tmp_ss.str(in);
+	return tmp_ss;
+}
+
+int main(int argc, char *argv[])
+{
+	int pnum, tnum; // Size of the vertexes and tetrahedrons
+	point3d *modelp = nullptr; // vertex class
+	tetrahedron *modelt = nullptr; // tetrahedron class
+
+	// file object
+	std::ifstream infile("cube.msh");
+
+	int tmp_i;
+	int tmp_tnum;
+	std::string line_str;
+	std::stringstream line_ss;
+	while (getline(infile, line_str))
+	{
+		if (line_str == "$Nodes") // Read vertexes
+		{
+			// Get vertex size firstly
+			getline(infile, line_str);
+			line_ss = str2ss(line_str);
+			line_ss >> pnum;
+			// Allocate memories
+			modelp = new point3d [pnum];
+			for (int i = 0; i < pnum; i++)
+			{
+				getline(infile, line_str);
+				line_ss = str2ss(line_str);
+				line_ss >> tmp_i >> modelp[i].x >> modelp[i].y >> modelp[i].z;
+			}
+		}
+		else if (line_str == "$Elements") // Read tetrahedral indexes
+		{
+			// Get tetrahedron's sizes
+			getline(infile, line_str);
+			line_ss = str2ss(line_str);
+			line_ss >> tmp_tnum;
+			tnum = tmp_tnum;
+
+			for (int i = 0; i < tmp_tnum; i++)
+			{
+				getline(infile, line_str);
+				line_ss = str2ss(line_str);
+				line_ss >> tmp_i >> tmp_i;
+				if (tmp_i != 4) // The tag for tetrahedron is 4
+				{
+					tnum--;
+				}
+				else break;
+			}
+
+			// Allocate memories
+			modelt = new tetrahedron [tnum];
+
+			line_ss = str2ss(line_str);
+
+			for (int j = 0; j < 5; j++) // Skip the first 5 numbers
+				line_ss >> tmp_i;
+			for (int j = 0; j < 4; j++)
+			{
+				line_ss >> tmp_i;
+				modelt[0].vec_ptr[j] = &modelp[tmp_i-1];
+			}
+
+			for (int i = 1; i < tnum; i++)
+			{
+				getline(infile, line_str);
+				line_ss = str2ss(line_str);
+
+				for (int j = 0; j < 5; j++) // Skip the first 5 numbers
+					line_ss >> tmp_i;
+				for (int j = 0; j < 4; j++)
+				{
+					line_ss >> tmp_i;
+					modelt[i].vec_ptr[j] = &modelp[tmp_i-1];
+				}
+			}
+		}
+	}
+
+	infile.close();
+
+	// Output for hunman validation
+	std::clog << "node number = " << pnum << std::endl;
+	std::clog << "tetra number = " << tnum << std::endl;
+
+	// Set tetrahedron's density and initialize it
+	for (int i = 0; i < tnum; i++)
+	{
+		modelt[i].rho = 1.0;
+		modelt[i].initialize_tensors();
+	}
+
+	// Declare observation stations
+	double xmin = 0, dx = 2.5; // 41
+	double ymin = 0, dy = 2.5; // 41
+	point3d obs[1681]; // 41*41
+	for (int i = 0; i < 41; i++)
+	{
+		for (int j = 0; j < 41; j++)
+		{
+			obs[i*41+j].x = xmin + dx*j;
+			obs[i*41+j].y = ymin + dy*i;
+			obs[i*41+j].z = 0.0;
+		}
+	}
+
+	// Calculate the gravity value
+	double sum;
+	for (int i = 0; i < 1681; i++)
+	{
+		std::cout << obs[i].x << " " << obs[i].y << " ";
+
+		sum = 0.0;
+		for (int t = 0; t < tnum; t++)
+		{
+			sum += tetra_vz(&modelt[t], &obs[i]);
+		}
+		std::cout << sum << std::endl;
+	}
+
+	if (modelp != nullptr) delete[] modelp;
+	if (modelt != nullptr) delete[] modelt;
+	return 0;
+}

tetra_vz.h

@@ -0,0 +1,244 @@
+#ifndef _TETRA_VZ_H
+#define _TETRA_VZ_H
+
+#include "cmath"
+
+#define ZERO 1e-20
+#define G0 6.67408e-11
+
+/**
+ * Define the point structure in 3D space
+ */
+
+struct point3d
+{
+	double x, y, z;
+	double module() const; // module
+	point3d normal() const; // normal vector
+};
+
+double point3d::module() const
+{
+	return sqrt(x*x + y*y + z*z);
+}
+
+point3d point3d::normal() const
+{
+	point3d nor_v;
+	double length = module();
+	if (length <= ZERO)
+	{
+		nor_v.x = nor_v.y = nor_v.z = 0.0;
+	}
+	else
+	{
+		nor_v.x = x/length;
+		nor_v.y = y/length;
+		nor_v.z = z/length;
+	}
+	return nor_v;
+}
+
+point3d operator+(const point3d &a, const point3d &b) // Add
+{
+	point3d m;
+	m.x = a.x+b.x;
+	m.y = a.y+b.y;
+	m.z = a.z+b.z;
+	return m;
+}
+
+point3d operator-(const point3d &a, const point3d &b) // Minus
+{
+	point3d m;
+	m.x = a.x-b.x;
+	m.y = a.y-b.y;
+	m.z = a.z-b.z;
+	return m;
+}
+
+point3d operator*(double scale, const point3d &b) // Multiply
+{
+	point3d m;
+	m.x = scale*b.x;
+	m.y = scale*b.y;
+	m.z = scale*b.z;
+	return m;
+}
+
+double dot(const point3d &a, const point3d &b) // Dot
+{
+	return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+double distance(const point3d &a, const point3d &b) // Distance
+{
+	return (a - b).module();
+}
+
+point3d cross(const point3d &a, const point3d &b) // Cross
+{
+	point3d v;
+	v.x = a.y*b.z-a.z*b.y;
+	v.y = a.z*b.x-a.x*b.z;
+	v.z = a.x*b.y-a.y*b.x;
+	if (fabs(v.x) <= ZERO) v.x = 0;
+	if (fabs(v.y) <= ZERO) v.y = 0;
+	if (fabs(v.z) <= ZERO) v.z = 0;
+	return v;
+}
+
+/**
+ * Define the tensor structure
+ */
+
+struct tensor
+{
+	double val[3][3];
+};
+
+tensor kron(const point3d &a, const point3d &b) // tensor product
+{
+	tensor t;
+	t.val[0][0]=a.x*b.x; t.val[0][1]=a.x*b.y; t.val[0][2]=a.x*b.z;
+	t.val[1][0]=a.y*b.x; t.val[1][1]=a.y*b.y; t.val[1][2]=a.y*b.z;
+	t.val[2][0]=a.z*b.x; t.val[2][1]=a.z*b.y; t.val[2][2]=a.z*b.z;
+	return t;
+}
+
+/**
+ * Define the tetrahedron structure
+ * 
+ * regular type of tetrahedron
+ *          3
+ *         /|\
+ *        / | \    y
+ *   z   /  |  \  /
+ *   |  /   1   \
+ *   | /  /   \  \
+ *   |/ /       \ \
+ *   0-------------2----> x
+ * 
+ * triangle list (anti-clockwise)
+ * 0 1 2
+ * 0 2 3
+ * 0 3 1
+ * 1 3 2
+ * 
+ * inverse type of tetrahedron
+ *          3
+ *         /|\
+ *        / | \    y
+ *   z   /  |  \  /
+ *   |  /   2   \
+ *   | /  /   \  \
+ *   |/ /       \ \
+ *   0-------------1----> x
+ * 
+ * triangle list (anti-clockwise)
+ * 0 1 3
+ * 0 2 1
+ * 0 3 2
+ * 1 2 3
+ */
+static int regular_order[12] = {0,1,2,0,2,3,0,3,1,1,3,2};
+static int inverse_order[12] = {0,1,3,0,2,1,0,3,2,1,2,3};
+
+struct tetrahedron
+{
+	bool regular_type;
+	point3d *vec_ptr[4];
+	tensor F[4];
+	tensor E[12];
+	double rho;
+	double edglen[12];
+	void initialize_tensors();
+};
+
+void tetrahedron::initialize_tensors()
+{
+	point3d v1, v2, v3, nf, ne;
+
+	v1 = *vec_ptr[1] - *vec_ptr[0];
+	v2 = *vec_ptr[2] - *vec_ptr[0];
+	v3 = *vec_ptr[3] - *vec_ptr[0];
+
+	int *triangle_order;
+	if (dot(v3, cross(v1, v2)) <= 0.0)
+	{
+		regular_type = true;
+		triangle_order = regular_order;
+	}
+	else
+	{
+		regular_type = false;
+		triangle_order = inverse_order;
+	}
+
+	for (int i = 0; i < 4; i++)
+	{
+		v1 = *vec_ptr[triangle_order[1+i*3]] - *vec_ptr[triangle_order[i*3]];
+		v2 = *vec_ptr[triangle_order[2+i*3]] - *vec_ptr[triangle_order[i*3]];
+		nf = cross(v1, v2).normal();
+		F[i] = kron(nf, nf);
+
+		for (int j = 0; j < 3; j++)
+		{
+			edglen[j+i*3] = distance(*vec_ptr[triangle_order[j+i*3]], *vec_ptr[triangle_order[(j+1)%3+i*3]]);
+
+			v3 = *vec_ptr[triangle_order[(j+1)%3+i*3]] - *vec_ptr[triangle_order[j+i*3]];
+			ne = cross(v3, nf).normal();
+			E[j+i*3] = kron(nf, ne);
+		}
+	}
+	return;
+}
+
+double tetra_vz(tetrahedron *ele_ptr, point3d *op_ptr)
+{
+	int f,e;
+	double Le,wf;
+	double dv1,dv2;
+	double face_sum,edge_sum;
+	point3d re;
+	point3d r_ijk[3];
+	double L_ijk[3];
+
+	int *triangle_order;
+	if (ele_ptr->regular_type)
+		triangle_order = regular_order;
+	else
+		triangle_order = inverse_order;
+
+	face_sum = edge_sum = 0.0;
+	for (f = 0; f < 4; f++)
+	{
+		r_ijk[0] = *ele_ptr->vec_ptr[triangle_order[3*f]] - *op_ptr;
+		r_ijk[1] = *ele_ptr->vec_ptr[triangle_order[3*f+1]] - *op_ptr;
+		r_ijk[2] = *ele_ptr->vec_ptr[triangle_order[3*f+2]] - *op_ptr;
+
+		L_ijk[0] = r_ijk[0].module();
+		L_ijk[1] = r_ijk[1].module();
+		L_ijk[2] = r_ijk[2].module();
+
+		wf =2*atan2(dot(r_ijk[0], cross(r_ijk[1],r_ijk[2])),
+					L_ijk[0]*L_ijk[1]*L_ijk[2] + L_ijk[0]*dot(r_ijk[1],r_ijk[2]) + 
+					L_ijk[1]*dot(r_ijk[2],r_ijk[0]) + L_ijk[2]*dot(r_ijk[0],r_ijk[1]));
+
+		face_sum += (ele_ptr->F[f].val[2][0]*r_ijk[0].x + ele_ptr->F[f].val[2][1]*r_ijk[0].y + ele_ptr->F[f].val[2][2]*r_ijk[0].z) * wf;
+
+		for (e = 0; e < 3; e++)
+		{
+			dv1 = distance(*ele_ptr->vec_ptr[triangle_order[e+3*f]], *op_ptr);
+			dv2 = distance(*ele_ptr->vec_ptr[triangle_order[(e+1)%3+3*f]], *op_ptr);
+
+			re = 0.5*(*ele_ptr->vec_ptr[triangle_order[e+3*f]] + *ele_ptr->vec_ptr[triangle_order[(e+1)%3+3*f]]) - *op_ptr;
+			Le = log((dv1+dv2+ele_ptr->edglen[e+3*f])/(dv1+dv2-ele_ptr->edglen[e+3*f]));
+
+			edge_sum += (ele_ptr->E[e+3*f].val[2][0]*re.x + ele_ptr->E[e+3*f].val[2][1]*re.y + ele_ptr->E[e+3*f].val[2][2]*re.z) * Le;
+		}
+	}
+	return 1.0e+8*G0*(ele_ptr->rho)*(edge_sum - face_sum);
+}
+
+#endif // _TETRA_VZ_H