/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include <algorithm>
#include <string>
#include <iostream>
#include <cmath>
#include <map>
#include <vector>

#include "geometry2d.h"

// 计算点t是否在直线ab上。
bool gctl::geometry2d::collinear(const point2dc &a, const point2dc &b, 
	const point2dc &t, double cut_off)
{
	point2dc at = t - a;
	point2dc ab = b - a;
	double cos_arc = dot(at, ab)/(at.module()*ab.module());
	if (sqrt(1.0 - cos_arc*cos_arc)*at.module() <= cut_off)
		return true;
	else return false;
}

// 计算直线（射线）外一点在线上的投影点
bool gctl::geometry2d::dot_on_line(const point2dc &p1, const point2dc &p2, const point2dc &out_p, point2dc &on_p)
{
	if (p1 == p2) // 线上两点为同一点 线外点的投影位置也为此同一点 视为在射线正方向
	{
		on_p = p2; return true;
	}

	if (p1 == out_p) // 线外点与第一个点重合
	{
		on_p = p1; return true;
	}
	
	if (p2 == out_p) // 线外点与第二个点重合
	{
		on_p = p2;
		return true;
	}

	gctl::point2dc line = p2 - p1;
	double mod = dot(line, out_p - p1)/line.module();
	on_p = mod * line.normal() + p1;

	if (mod < 0.0) return false;
	return true;
}

// 求两条直线的交点。
void gctl::geometry2d::line_intersect(const point2dc &l1_p1, const point2dc &l1_p2, const point2dc &l2_p1, const point2dc &l2_p2, point2dc &out)
{
	// 如果两条直线平行则无交点
	if (fabs((l1_p2.x - l1_p1.x)*(l2_p2.y - l2_p1.y) - (l1_p2.y - l1_p1.y)*(l2_p2.x - l2_p1.x)) > GCTL_ZERO)
	{
		double sigma = ((l2_p2.x - l2_p1.x)*(l1_p1.y - l2_p1.y) - (l2_p2.y - l2_p1.y)*(l1_p1.x - l2_p1.x))
		/((l1_p2.x - l1_p1.x)*(l2_p2.y - l2_p1.y) - (l1_p2.y - l1_p1.y)*(l2_p2.x - l2_p1.x));
		out.x = l1_p1.x + sigma*(l1_p2.x - l1_p1.x);
		out.y = l1_p1.y + sigma*(l1_p2.y - l1_p1.y);
		return;
	}
}

// 求两条直线的交点。
gctl::point2dc gctl::geometry2d::line_intersect(const point2dc &l1_p1, const point2dc &l1_p2, 
	const point2dc &l2_p1, const point2dc &l2_p2)
{
	point2dc out;
	if (fabs((l1_p2.x - l1_p1.x)*(l2_p2.y - l2_p1.y) - (l1_p2.y - l1_p1.y)*(l2_p2.x - l2_p1.x)) <= GCTL_ZERO)
		return out;
	else
	{
		double sigma = ((l2_p2.x - l2_p1.x)*(l1_p1.y - l2_p1.y) - (l2_p2.y - l2_p1.y)*(l1_p1.x - l2_p1.x))
		/((l1_p2.x - l1_p1.x)*(l2_p2.y - l2_p1.y) - (l1_p2.y - l1_p1.y)*(l2_p2.x - l2_p1.x));
		out.x = l1_p1.x + sigma*(l1_p2.x - l1_p1.x);
		out.y = l1_p1.y + sigma*(l1_p2.y - l1_p1.y);
		return out;
	}
}

//计算三角形内一点到三条边的垂线讲三角形分割所得三个菱形区域的面积
void gctl::geometry2d::tri_split_area(const point2dc &t_p1, const point2dc &t_p2, 
	const point2dc &t_p3, const point2dc &in_p, double &a1, double &a2, double &a3)
{
	point2dc t[3] = {t_p1, t_p2, t_p3};
	double a[3] = {0.0, 0.0, 0.0};

	point2dc l1, l2;
	double cos_edge, sin_edge;
	for (int i = 0; i < 3; i++)
	{
		if (in_p != t[i] && in_p != t[(i+1)%3])
		{
			l1 = t[(i+1)%3] - t[i];
			l2 = in_p - t[i];
			cos_edge = dot(l1, l2)/l1.module();
			sin_edge = sqrt(l2.x*l2.x + l2.y*l2.y - cos_edge*cos_edge);
			a[i] += 0.5*cos_edge*sin_edge;

			l1 = t[i] - t[(i+1)%3];
			l2 = in_p - t[(i+1)%3];
			cos_edge = dot(l1, l2)/l1.module();
			sin_edge = sqrt(l2.x*l2.x + l2.y*l2.y - cos_edge*cos_edge);
			a[(i+1)%3] += 0.5*cos_edge*sin_edge;
		}
	}

	a1 = a[0]; a2 = a[1]; a3 = a[2];
	return;
}

//细分三角形面积比形式的三角形内插值函数
double gctl::geometry2d::tri_interp_area(const point2dc &p, const point2dc &p1, const point2dc &p2, 
	const point2dc &p3, double d1, double d2, double d3)
{
	point2dc pp1 = p1 - p;
	point2dc pp2 = p2 - p;
	point2dc pp3 = p3 - p;
	//三角形的面积等于叉乘的1/2 这里只计算比值 所以不需要乘以1/2
	double a1 = fabs(cross(pp2,pp3));
	double a2 = fabs(cross(pp3,pp1));
	double a3 = fabs(cross(pp1,pp2));
	return (d1*a1+d2*a2+d3*a3)/(a1+a2+a3);
}

// 距离反比形式的三角形内插值函数
double gctl::geometry2d::tri_interp_dist(const point2dc &p, const point2dc &p1, 
	const point2dc &p2, const point2dc &p3, double d1, double d2, double d3, double w)
{
	double w1 = 1.0/(pow((p - p1).module(), w) + GCTL_ZERO);
	double w2 = 1.0/(pow((p - p2).module(), w) + GCTL_ZERO);
	double w3 = 1.0/(pow((p - p3).module(), w) + GCTL_ZERO);
	return (d1*w1 + d2*w2 + d3*w3)/(w1 + w2 + w3);
}

// 利用平面方程计算二维坐标上三点数据所确定的方向导数
void gctl::geometry2d::tri_plane_gradient(const point2dc &p1, const point2dc &p2, 
	const point2dc &p3, const double &d1, const double &d2, const double &d3, 
	double &dx, double &dy)
{
	double A = (p2.y - p1.y)*(d3 - d1) - (d2 - d1)*(p3.y - p1.y);
	double B = (p3.x - p1.x)*(d2 - d1) - (d3 - d1)*(p2.x - p1.x);
	double C = (p2.x - p1.x)*(p3.y - p1.y) - (p3.x - p1.x)*(p2.y - p1.y);
	dx = -1.0*A/C;
	dy = -1.0*B/C;
	return;
}

// 判断一个点是否在三角形内（二维）
bool gctl::geometry2d::node_in_triangle(const point2dc &test_p, const point2dc &p1, 
	const point2dc &p2, const point2dc &p3, bool on_boundary)
{
	point2dc l1, l2;
	point2dc corner[3] = {p1, p2, p3};

	if (cross(corner[1] - corner[0], corner[2] - corner[0]) > 0)
	{
		if (on_boundary)
		{
			for (int i = 0; i < 3; i++)
			{
				l1 = corner[(i+1)%3] - corner[i];
				l2 = test_p - corner[i];

				if (cross(l1, l2) < 0)
					return false;
			}
		}
		else
		{
			for (int i = 0; i < 3; i++)
			{
				l1 = corner[(i+1)%3] - corner[i];
				l2 = test_p - corner[i];

				if (cross(l1, l2) <= 0)
					return false;
			}
		}
	}
	else if (cross(corner[1] - corner[0], corner[2] - corner[0]) < 0)
	{
		if (on_boundary)
		{
			for (int i = 0; i < 3; i++)
			{
				l1 = corner[(i+1)%3] - corner[i];
				l2 = test_p - corner[i];

				if (cross(l1, l2) > 0)
					return false;
			}
		}
		else
		{
			for (int i = 0; i < 3; i++)
			{
				l1 = corner[(i+1)%3] - corner[i];
				l2 = test_p - corner[i];

				if (cross(l1, l2) >= 0)
					return false;
			}
		}
	}
	else
	{
		std::string err_str = "Invalid triangular nodes. From gctl::geometry2d::node_in_triangle2d(...)";
		throw runtime_error(err_str);
	}

	return true;
}

bool gctl::geometry2d::node_in_polygon(const array<point2dc> &poly_ps, const point2dc &one_p, 
	bool on_boundary)
{
	int count, pnum;
	double tmp_x;
	double xmin,xmax,ymin,ymax;

	std::string err_str;
	if (poly_ps.size() <= 2)
	{
		err_str = "Invalid polygon nodes. From gctl::geometry2d::node_in_polygon(...)";
		throw runtime_error(err_str);
		return false;
	}
	else
	{
		pnum = poly_ps.size();
		//确定外接矩形
		xmin = GCTL_BDL_MAX; ymin = GCTL_BDL_MAX;
		xmax = GCTL_BDL_MIN; ymax = GCTL_BDL_MIN;
		for (int j = 0; j < pnum; j++)
		{
			xmin = GCTL_MIN(xmin, poly_ps[j].x);
			xmax = GCTL_MAX(xmax, poly_ps[j].x);
			ymin = GCTL_MIN(ymin, poly_ps[j].y);
			ymax = GCTL_MAX(ymax, poly_ps[j].y);
		}

		// 测试点在外接矩形外 直接返回假
		if (one_p.x < xmin || one_p.x > xmax || one_p.y < ymin || one_p.y > ymax)
		{
			return false;
		}
		else
		{
			count = 0;
			if (on_boundary)
			{
				for (int i = 0; i < pnum; i++)
				{
					// 包含点落在边上的情况
					if ((one_p.y >= poly_ps[i].y && one_p.y <= poly_ps[(i+1)%pnum].y) ||
						(one_p.y <= poly_ps[i].y && one_p.y >= poly_ps[(i+1)%pnum].y))
					{
						tmp_x = (poly_ps[(i+1)%pnum].x - poly_ps[i].x) 
							* (one_p.y - poly_ps[i].y)/(poly_ps[(i+1)%pnum].y - poly_ps[i].y) 
							+ poly_ps[i].x;

						if (one_p.x <= tmp_x)
							count++;
					}
				}
			}
			else
			{
				for (int i = 0; i < pnum; i++)
				{
					// 包含点落在边上的情况
					if ((one_p.y >= poly_ps[i].y && one_p.y < poly_ps[(i+1)%pnum].y) ||
						(one_p.y <= poly_ps[i].y && one_p.y > poly_ps[(i+1)%pnum].y))
					{
						tmp_x = (poly_ps[(i+1)%pnum].x - poly_ps[i].x) 
							* (one_p.y - poly_ps[i].y)/(poly_ps[(i+1)%pnum].y - poly_ps[i].y) 
							+ poly_ps[i].x;

						if (one_p.x < tmp_x)
							count++;
					}
				}
			}

			if (pow(-1,count) < 0) return true;
			else return false;
		}
	}
}

void gctl::geometry2d::get_common_edges(const array<triangle2d> &input_ele, array<edge2d> &out_edge)
{
	if (input_ele.empty())
	{
		std::string err_str = "The input mesh is empty. From gctl::geometry2d::get_common_edges(...)";
		throw runtime_error(err_str);
	}

	int ele_size = input_ele.size();
	// find edge neighbors
	std::vector<std::vector<int> > edge_n;
	edge_n.resize(ele_size);
	for (int i = 0; i < ele_size; i++)
	{
		edge_n[i].reserve(6); // 每个顶点的list最多6个
	}

	int small_id, big_id;
	for (int i = 0; i < ele_size; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			if (input_ele[i].neigh[j] != nullptr) // 邻居必须存在
			{
				// 这里我们总是按小-大排序序号，保证不会输入混乱
				small_id = std::min(input_ele[i].id, input_ele[i].neigh[j]->id);
				big_id   = std::max(input_ele[i].id, input_ele[i].neigh[j]->id);
				edge_n[small_id].push_back(big_id);
			}
		}
	}

	// 去除edge_n中的重复输入
	std::vector<int>::iterator pos;
	for (int i = 0; i < ele_size; i++)
	{
		if (!edge_n[i].empty()) // 很有可能某个边没有输入
		{
			std::sort(edge_n[i].begin(), edge_n[i].end()); //对顶点序列由小到大排序
			pos = std::unique(edge_n[i].begin(), edge_n[i].end()); //获取重复序列开始的位置
			edge_n[i].erase(pos, edge_n[i].end()); //删除重复点
		}
	}

	// 保存三角形的公共边的列表 首先统计公共边的数量
	int edge_num = 0;
	for (int i = 0; i < ele_size; i++)
		edge_num += edge_n[i].size();
	out_edge.resize(edge_num);

	int count = 0;
	for (int i = 0; i < ele_size; i++)
	{
		if (!edge_n[i].empty())
		{
			for (int j = 0; j < edge_n[i].size(); j++)
			{
				out_edge[count].id = count;
				out_edge[count].neigh[0] = input_ele.get(i);
				out_edge[count].neigh[1] = input_ele.get(edge_n[i][j]);
				count++;
			}
		}
	}

	// 最后找到边的两个顶点
	// 如果相邻的两个三角形中的其中一个三角形的一个顶点不等于另一个三角形的所有顶点，则这个三角形的
	// 另两个顶点连成的边即为公共边
	int tmp_int;
	for (int i = 0; i < edge_num; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			tmp_int = out_edge[i].neigh[0]->vert[j]->id;
			if (tmp_int != out_edge[i].neigh[1]->vert[0]->id &&
				tmp_int != out_edge[i].neigh[1]->vert[1]->id &&
				tmp_int != out_edge[i].neigh[1]->vert[2]->id)
			{
				out_edge[i].vert[0] = out_edge[i].neigh[0]->vert[(j+1)%3];
				out_edge[i].vert[1] = out_edge[i].neigh[0]->vert[(j+2)%3];
				break;
			}
		}
	}

	//清除向量
	destroy_vector(edge_n);
	return;
}

int gctl::geometry2d::sort_node_number(const array<triangle2d> &input_ele)
{
	if (input_ele.empty())
	{
		std::string err_str = "The input mesh is empty. From gctl::geometry2d::sort_node_number(...)";
		throw runtime_error(err_str);
		return 0;
	}

	std::vector<int> node_index;
	node_index.reserve(3*input_ele.size());
	for (int i = 0; i < input_ele.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			node_index.push_back(input_ele[i].vert[j]->id);
		}
	}

	std::vector<int>::iterator pos;
	std::sort(node_index.begin(), node_index.end()); //对顶点序列由小到大排序
	pos = std::unique(node_index.begin(), node_index.end()); //获取重复序列开始的位置
	node_index.erase(pos, node_index.end()); //删除重复点

	int out_num = node_index.size();
	destroy_vector(node_index);
	return out_num;
}

void gctl::geometry2d::sort_node_neighbor(const array<triangle2d> &input_ele, 
	std::vector<std::vector<triangle2d*> > &out_node_neigh, int node_num)
{
	// get node number
	if (node_num < 0)
	{
		throw std::runtime_error("Invalid node size. From gctl::geometry2d::sort_node_neighbor(...)");
	}

	out_node_neigh.resize(node_num);
	for (int i = 0; i < input_ele.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			out_node_neigh[input_ele[i].vert[j]->id].push_back(input_ele.get(i));
		}
	}
	return;
}

void gctl::geometry2d::sort_triangle_neighbor(array<triangle2d> &input_ele, const std::vector<std::vector<triangle2d*> > &node_neighs)
{
	bool double_break;
	int loop_size; //equal_num;
	for (int i = 0; i < node_neighs.size(); i++)
	{
		loop_size = node_neighs[i].size();
		for (int j = 0; j < loop_size-1; j++)
		{
			for (int k = j+1; k < loop_size; k++)
			{
				double_break = false;
				for (int p = 0; p < 3; p++)
				{
					for (int q = 0; q < 3; q++)
					{
						// 顶点必须是逆时针排序的
						if (node_neighs[i][j]->vert[p] == node_neighs[i][k]->vert[(q+1)%3] && 
							node_neighs[i][j]->vert[(p+1)%3] == node_neighs[i][k]->vert[q])
						{
							node_neighs[i][j]->neigh[p] = node_neighs[i][k];
							node_neighs[i][k]->neigh[q] = node_neighs[i][j];
							double_break = true;
							break;
						}
					}

					if (double_break) break;
				}

				/*
				equal_num = 0;
				for (int p = 0; p < 3; p++)
				{
					for (int q = 0; q < 3; q++)
					{
						if (node_neighs[i][j]->vert[p] == node_neighs[i][k]->vert[q])
						{
							equal_num++;
							break;
						}
					}
				}

				if (equal_num == 2)
				{
					// 相邻三角形可能被公共边的两个顶点重复添加
					if (node_neighs[i][k] != node_neighs[i][j]->neigh[0] &&
						node_neighs[i][k] != node_neighs[i][j]->neigh[1] &&
						node_neighs[i][k] != node_neighs[i][j]->neigh[2])
					{
						// loop afterward to free space for new input
						node_neighs[i][j]->neigh[2] = node_neighs[i][j]->neigh[1];
						node_neighs[i][j]->neigh[1] = node_neighs[i][j]->neigh[0];
						node_neighs[i][j]->neigh[0] = node_neighs[i][k];
					}

					if (node_neighs[i][j] != node_neighs[i][k]->neigh[0] &&
						node_neighs[i][j] != node_neighs[i][k]->neigh[1] &&
						node_neighs[i][j] != node_neighs[i][k]->neigh[2])
					{
						node_neighs[i][k]->neigh[2] = node_neighs[i][k]->neigh[1];
						node_neighs[i][k]->neigh[1] = node_neighs[i][k]->neigh[0];
						node_neighs[i][k]->neigh[0] = node_neighs[i][j];
					}
				}
				*/
			}
		}
	}

	destroy_vector(node_neighs);
	return;
}

// 使用直线切割二维三角网络
void gctl::geometry2d::cut_triangular_mesh_2d(const array<triangle2d> &in_eles, const point2dc &p1, 
	const point2dc &p2, array<vertex2dc> &out_nodes, array<triangle2d> &out_eles)
{
	bool created;
	point2dc tmp_p, tmp_p2;
	vertex2dc *tmp_vert, *tmp_vert2;
	triangle2d tmp_tri2d;
	std::vector<triangle2d> ele_index;
	std::vector<vertex2dc*> extra_nodes; // 注意这里是新生的顶点的指针向量

	int left_count; // 统计三角形顶点在直线左侧的个数
	for (int i = 0; i < in_eles.size(); ++i)
	{
		left_count = 0;
		for (int j = 0; j < 3; ++j)
		{
			// 注意这里包含等号 即点在直线上也等于在线的左侧
			// 此种设计会导致一点或两点在线上，其他点在右侧也会被记入在内
			if (cross(p2 - p1, *in_eles[i].vert[j] - p1) >= 0) left_count++;
		}

		if (left_count == 3) // 全在左侧 添加到输出列表
		{
			tmp_tri2d = in_eles[i];
			tmp_tri2d.id = ele_index.size();
			ele_index.push_back(tmp_tri2d);
		}

		if (left_count == 2) // 两个点在左侧或者线上 计算两个交点的位置添加到顶点列表 并新建两个三角形
		{
			for (int j = 0; j < 3; ++j)
			{
				if (cross(p2 - p1, *in_eles[i].vert[j] - p1) < 0) // 找到直线右侧的点 这里不包含等号 必须要完全在右侧
				{
					line_intersect(*in_eles[i].vert[j], *in_eles[i].vert[(j+1)%3], p1, p2, tmp_p);
					// 新建顶点 把顶点的指针交给 extra_nodes 保管
					// 检查交点是否已经存在
					created = false;
					tmp_p2.x = in_eles[i].vert[(j+1)%3]->x;
					tmp_p2.y = in_eles[i].vert[(j+1)%3]->y;
					if (!isequal(tmp_p, tmp_p2, 1e-6)) // 如果交点1等于顶点eles[i].vert[(j+1)%3] 不需要新建顶点以及三角形
					{
						// 检查交点是否已经存在
						for (int e = 0; e < extra_nodes.size(); ++e)
						{
							tmp_p2.x = extra_nodes[e]->x;
							tmp_p2.y = extra_nodes[e]->y;
							if (isequal(tmp_p, tmp_p2, 1e-6))
							{
								tmp_vert = extra_nodes[e];
								created = true;
								break;
							}
						}

						// 确定点不存在则新建
						if (!created)
						{
							tmp_vert = new vertex2dc;
							tmp_vert->id = 0; // 先统一把顶点索引设置为0 后面再改
							tmp_vert->x  = tmp_p.x;
							tmp_vert->y  = tmp_p.y;
							extra_nodes.push_back(tmp_vert);
						}

						tmp_tri2d.id = ele_index.size();
						tmp_tri2d.vert[0] = tmp_vert;
						tmp_tri2d.vert[1] = in_eles[i].vert[(j+1)%3];
						tmp_tri2d.vert[2] = in_eles[i].vert[(j+2)%3];
						ele_index.push_back(tmp_tri2d);
					}
					else tmp_vert = in_eles[i].vert[(j+1)%3]; // 注意即使不新建顶点也要对临时顶点赋值

					line_intersect(*in_eles[i].vert[j], *in_eles[i].vert[(j+2)%3], p1, p2, tmp_p);
					// 新建顶点 把顶点的指针交给 extra_nodes 保管
					created = false;
					tmp_p2.x = in_eles[i].vert[(j+2)%3]->x;
					tmp_p2.y = in_eles[i].vert[(j+2)%3]->y;
					if (!isequal(tmp_p, tmp_p2, 1e-6)) // 如果交点1等于顶点eles[i].vert[(j+2)%3] 不需要新建顶点以及三角形
					{
						for (int e = 0; e < extra_nodes.size(); ++e)
						{
							tmp_p2.x = extra_nodes[e]->x;
							tmp_p2.y = extra_nodes[e]->y;
							if (isequal(tmp_p, tmp_p2, 1e-6))
							{
								tmp_vert2 = extra_nodes[e];
								created = true;
								break;
							}
						}

						if (!created)
						{
							tmp_vert2 = new vertex2dc;
							tmp_vert2->id = 0; // 先统一把顶点索引设置为0 后面再改
							tmp_vert2->x  = tmp_p.x;
							tmp_vert2->y  = tmp_p.y;
							extra_nodes.push_back(tmp_vert2);
						}

						// 新建三角形
						tmp_tri2d.id = ele_index.size();
						tmp_tri2d.vert[0] = in_eles[i].vert[(j+2)%3];
						tmp_tri2d.vert[1] = tmp_vert2;
						tmp_tri2d.vert[2] = tmp_vert;
						ele_index.push_back(tmp_tri2d);
					}

					break; // 完成跳出
				}
			}
		}

		if (left_count == 1)
		{
			for (int j = 0; j < 3; ++j)
			{
				if (cross(p2 - p1, *in_eles[i].vert[j] - p1) > 0) // 找到直线左侧的点 注意这里不包含等号 因为我们需要点完全在线的左侧 不包括刚好在线上的情况
				{
					line_intersect(*in_eles[i].vert[j], *in_eles[i].vert[(j+1)%3], p1, p2, tmp_p);
					created = false;
					if (!created)
					{
						for (int e = 0; e < extra_nodes.size(); ++e)
						{
							tmp_p2.x = extra_nodes[e]->x;
							tmp_p2.y = extra_nodes[e]->y;
							if (isequal(tmp_p, tmp_p2, 1e-6))
							{
								tmp_vert = extra_nodes[e];
								created = true;
								break;
							}
						}
					}

					if (!created)
					{
						tmp_vert = new vertex2dc;
						tmp_vert->id = 0; // 先统一把顶点索引设置为0 后面再改
						tmp_vert->x  = tmp_p.x;
						tmp_vert->y  = tmp_p.y;
						extra_nodes.push_back(tmp_vert);
					}
					
					line_intersect(*in_eles[i].vert[j], *in_eles[i].vert[(j+2)%3], p1, p2, tmp_p);
					// 新建顶点 把顶点的指针交给 extra_nodes 保管
					created = false;
					if (!created)
					{
						for (int e = 0; e < extra_nodes.size(); ++e)
						{
							tmp_p2.x = extra_nodes[e]->x;
							tmp_p2.y = extra_nodes[e]->y;
							if (isequal(tmp_p, tmp_p2, 1e-6))
							{
								tmp_vert2 = extra_nodes[e];
								created = true;
								break;
							}
						}
					}

					if (!created)
					{
						tmp_vert2 = new vertex2dc;
						tmp_vert2->id = 0; // 先统一把顶点索引设置为0 后面再改
						tmp_vert2->x  = tmp_p.x;
						tmp_vert2->y  = tmp_p.y;
						extra_nodes.push_back(tmp_vert2);
					}

					// 新建三角形
					tmp_tri2d.id = ele_index.size();
					tmp_tri2d.vert[0] = in_eles[i].vert[j];
					tmp_tri2d.vert[1] = tmp_vert;
					tmp_tri2d.vert[2] = tmp_vert2;
					ele_index.push_back(tmp_tri2d);

					break; // 完成跳出
				}
			}
		}
	}

	//整理出新的顶点与元素列表
	std::vector<vertex2dc*> node_index;
	out_eles.resize(ele_index.size());
	for (int i = 0; i < ele_index.size(); ++i)
	{
		out_eles[i] = ele_index[i];
		out_eles[i].id = i;
		for (int j = 0; j < 3; ++j)
		{
			node_index.push_back(ele_index[i].vert[j]);
		}
	}

	std::vector<vertex2dc*>::iterator pos;
	std::sort(node_index.begin(), node_index.end()); //对顶点序列由小到大排序
	pos = std::unique(node_index.begin(), node_index.end()); //获取重复序列开始的位置
	node_index.erase(pos, node_index.end()); //删除重复点

	std::map<vertex2dc*, vertex2dc*> vert_map;
	out_nodes.resize(node_index.size());
	for (int i = 0; i < node_index.size(); ++i)
	{
		out_nodes[i] = *node_index[i];
		out_nodes[i].id = i;
		vert_map[node_index[i]] = out_nodes.get(i); // 建立新旧顶点的映射
	}

	for (int i = 0; i < out_eles.size(); ++i)
	{
		for (int j = 0; j < 3; ++j)
		{
			out_eles[i].vert[j] = vert_map[out_eles[i].vert[j]]; // 将旧的顶点映射为新的顶点
		}
	}

	// 手动删除新建的顶点
	for (int i = 0; i < extra_nodes.size(); ++i)
	{
		delete extra_nodes[i];
	}

	return;
}

//提取二维三角网络的自网络（一部分）
void gctl::geometry2d::extract_triangular_mesh_2d(const array<triangle2d> &in_eles, const array<bool> &out_list, 
	array<vertex2dc> &out_nodes, array<triangle2d> &out_eles)
{
	if (in_eles.empty() || in_eles.size() != out_list.size())
	{
		throw runtime_error("Invalid input arrays (empty or incompatible sizes). From geometry2d::extract_triangular_mesh_2d(...)");
	}

	int valid_count = 0;
	for (int i = 0; i < in_eles.size(); ++i)
	{
		if (out_list[i]) valid_count++;
	}

	//整理出新的顶点与元素列表
	int count = 0;
	std::vector<vertex2dc*> node_index;
	out_eles.resize(valid_count);
	for (int i = 0; i < in_eles.size(); ++i)
	{
		if (out_list[i])
		{
			out_eles[count] = in_eles[i];
			out_eles[count].id = i;
			for (int j = 0; j < 3; ++j)
			{
				node_index.push_back(in_eles[i].vert[j]);
			}
			count++;
		}
	}

	std::vector<vertex2dc*>::iterator pos;
	std::sort(node_index.begin(), node_index.end()); //对顶点序列由小到大排序
	pos = std::unique(node_index.begin(), node_index.end()); //获取重复序列开始的位置
	node_index.erase(pos, node_index.end()); //删除重复点

	std::map<vertex2dc*, vertex2dc*> vert_map;
	out_nodes.resize(node_index.size());
	for (int i = 0; i < node_index.size(); ++i)
	{
		out_nodes[i] = *node_index[i];
		out_nodes[i].id = i;
		vert_map[node_index[i]] = out_nodes.get(i); // 建立新旧顶点的映射
	}

	for (int i = 0; i < out_eles.size(); ++i)
	{
		for (int j = 0; j < 3; ++j)
		{
			out_eles[i].vert[j] = vert_map[out_eles[i].vert[j]]; // 将旧的顶点映射为新的顶点
		}
	}

	return;
}