Couette-Meshless-Multigrid/general_computation_functions.cpp at main · comp-physics/Couette-Meshless-Multigrid · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include "general_computation_functions.h"
using std::cout;
using std::endl;
double distance(Point refPoint, Point queryPoint) {
	return std::sqrt(std::pow(std::get<0>(refPoint) - std::get<0>(queryPoint), 2) + std::pow(std::get<1>(refPoint) - std::get<1>(queryPoint), 2));
}

std::pair<double, double> minMaxCoord(vector<Point> points_, char coord){
	double min;
	double max;
	double tempCoord;
	if (coord == 'x') {
		min = std::get<0>(points_[0]);
		max = min;
	}
	else if (coord == 'y') {
		min = std::get<1>(points_[0]);
		max = min;
	}
	for (size_t i = 0; i < points_.size(); i++) {
		if (coord == 'x') {
			tempCoord = std::get<0>(points_[i]);
		}
		else if (coord == 'y') {
			tempCoord = std::get<1>(points_[i]);
		}
		if (tempCoord > max) {
			max = tempCoord;
		}
		else if (tempCoord < min) {
			min = tempCoord;
		}
	}
	return std::pair<double, double>(min, max);
}
Point vec_from_pts(Point point1, Point point2) {
	double xdiff = std::get<0>(point2) - std::get<0>(point1);
	double ydiff = std::get<1>(point2) - std::get<1>(point1);
	double zdiff = std::get<2>(point2) - std::get<2>(point1);

	//cout << "xdiff: " << xdiff << " ydiff: " << ydiff << endl;
	return std::make_tuple(xdiff, ydiff, zdiff);
}
Point midpoint(Point point1, Point point2) {
	double xdiff = std::get<0>(point2) + std::get<0>(point1);
	double ydiff = std::get<1>(point2) + std::get<1>(point1);
	double zdiff = std::get<2>(point2) + std::get<2>(point1);

	//cout << "xdiff: " << xdiff << " ydiff: " << ydiff << endl;
	return std::make_tuple(xdiff/2, ydiff/2, zdiff/2);
}

Point centroid(Point point1, Point point2, Point point3) {
	double xdiff = std::get<0>(point3) + std::get<0>(point2) + std::get<0>(point1);
	double ydiff = std::get<0>(point3) + std::get<1>(point2) + std::get<1>(point1);

	return std::make_tuple(xdiff / 3, ydiff / 3, 0);
}

Point unit_normal_vec(Point vec, bool cw) {
	double x = std::get<0>(vec);
	double y = std::get<1>(vec);

	double norm = std::sqrt(x*x + y * y);
	x /= norm;
	y /= norm;
	return cw? std::make_tuple(-y, x, 0) : std::make_tuple(y, -x, 0);
}
Point avg_unit_norm_vec(Point norm1, Point norm2) {
	Point vec = norm1;
	vec = std::make_tuple(std::get<0>(norm1) + std::get<0>(norm2),
		std::get<1>(norm1) + std::get<1>(norm2), 0);

	double x = std::get<0>(vec);
	double y = std::get<1>(vec);

	double norm = std::sqrt(x*x + y * y);
	x /= norm;
	y /= norm;
	return std::make_tuple(-y, x, 0);
}
std::vector<Point> shifting_scaling(vector<Point> points_, Point evalPoint) {
	std::pair<double, double> minMaxX, minMaxY;
	minMaxX = minMaxCoord(points_, 'x');
	minMaxY = minMaxCoord(points_, 'y');
	double x, y, x_ss, y_ss, scale, minX, minY, maxX, maxY;
	minX = minMaxX.first;
	maxX = minMaxX.second;
	minY = minMaxY.first;
	maxY = minMaxY.second;
	scale = std::max(maxX - minX, maxY - minY);
	std::vector<Point> scaledPoints;
	for (size_t i = 0; i < points_.size(); i++) {
		x = std::get<0>(points_[i]);
		y = std::get<1>(points_[i]);
		x_ss = (x - minX) / scale;
		y_ss = (y - minY) / scale;
		scaledPoints.push_back(Point(x_ss, y_ss, 0));
	}
	scaledPoints.push_back(Point(scale, scale, scale));
	x = std::get<0>(evalPoint);
	y = std::get<1>(evalPoint);
	x_ss = (x - minX) / scale;
	y_ss = (y - minY) / scale;
	scaledPoints.push_back(Point(x_ss, y_ss, 0));
	return scaledPoints;
}
void cuthill_mckee_ordering(vector<vector<int>> &adjacency, vector<int> &order, vector <int> &bcFlags) {

	vector<bool> visited(order.size(), false);
	std::queue<int> queue;
	vector<int> new_order;
	int startPoint = 0;
	while(bcFlags[startPoint] != 0){
		startPoint++;
	}
	visited[startPoint] = true;
	queue.push(startPoint);
	int currPoint, adjPoint;
	while (!queue.empty()) {
		currPoint = queue.front();
		new_order.push_back(currPoint);
		queue.pop();
		for (size_t i = 0; i < adjacency[currPoint].size(); i++) {
			adjPoint = adjacency[currPoint].at(i);
			if (visited[adjPoint] == false) {
				visited[adjPoint] = true;
				queue.push(adjPoint);
			}
		}
	}

	for (int i = 0; i < order.size(); i++){
		if(bcFlags[i] == 2){
			new_order.push_back(i);
		}
	}
	order = new_order;
}
void reverse_cuthill_mckee_ordering(vector<vector<int>> &adjacency, vector<int> &order , vector<int> &bcFlags) {
	cuthill_mckee_ordering(adjacency, order, bcFlags);
	std::reverse(order.begin(), order.end());
}