Skip to content
Snippets Groups Projects
Commit cd11ade9 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Begin EMC algorithms [ci-skip] 

Debugging still in progress
parent f5979379
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/mesh/MeshSmootherEscobar.cpp
+ 333
− 62
View file @ cd11ade9
......@@ -460,6 +460,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
void
_getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
{
constexpr bool multi_connection = true;
const ConnectivityType& connectivity = m_given_mesh.connectivity();
auto is_boundary_node = connectivity.isBoundaryNode();
......@@ -478,6 +480,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
quality.fill(2);
auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
const bool print = (node_id == 4);
auto cell_list = node_to_cell_matrix[node_id];
auto node_number_in_cell = node_number_in_their_cells[node_id];
......@@ -492,11 +496,33 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const TinyMatrix<Dimension> inv_W = inverse(W);
std::array<TinyMatrix<Dimension>, Dimension> S_gradient =
{TinyMatrix<Dimension>{-1, -inv_sin_alpha + inv_tan_alpha, //
+0, +0}, //
TinyMatrix<Dimension>{+0, +0, //
-1, -inv_sin_alpha + inv_tan_alpha}};
std::array<TinyMatrix<Dimension>, Dimension> //
S_gradient = {TinyMatrix<Dimension>{-1, -1, //
+0, +0},
TinyMatrix<Dimension>{+inv_tan_alpha, +inv_tan_alpha, //
-inv_sin_alpha, -inv_sin_alpha}};
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
double sigma_min = std::numeric_limits<double>::max();
double delta = 0;
auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
double final_f = 0;
std::cout.precision(15);
for (size_t i_iter = 0; i_iter < 1; ++i_iter) {
if (print) {
std::cout << "------------------------------\n";
std::cout << "iter = " << i_iter << '\n';
std::cout << "x = " << x << '\n';
}
double f = 0;
TinyVector<Dimension> f_gradient = zero;
TinyMatrix<Dimension> f_hessian = zero;
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
......@@ -511,6 +537,15 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
a[1] - x[1], b[1] - x[1]};
S_list[i_cell] = A * inv_W;
if (print) {
std::cout << "--- S --- cell = " << i_cell << " (" << node_to_cell_matrix[node_id][i_cell] << ")\n";
std::cout << "a = " << a << '\n';
std::cout << "b = " << b << '\n';
std::cout << "x = " << x << '\n';
std::cout << "A = " << A << '\n';
std::cout << "S = " << S_list[i_cell] << '\n';
}
}
SmallArray<double> sigma_list(S_list.size());
......@@ -518,54 +553,104 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
sigma_list[i_cell] = det(S_list[i_cell]);
}
double sigma_min = min(sigma_list);
double delta =
(sigma_min < eps) ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min)) : 0;
auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
{
double current_sigma_min = min(sigma_list);
double final_f = 0;
if constexpr (multi_connection) {
size_t nb_medium_range_simplices = 0;
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
nb_medium_range_simplices += 2 * (cell_to_node_matrix[cell_id].size() - 3);
}
for (size_t i_iter = 0; i_iter < 200; ++i_iter) {
double f = 0;
TinyVector<Dimension> f_gradient = zero;
TinyMatrix<Dimension> f_hessian = zero;
if (nb_medium_range_simplices > 0) {
SmallArray<double> medium_range_sigma_list(nb_medium_range_simplices);
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
const size_t i_cell_node = node_number_in_cell[i_cell];
auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]];
const size_t cell_nb_nodes = cell_node_list.size();
const double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
const double mr_sin_alpha = std::sin(mr_alpha);
const double mr_cos_alpha = std::cos(mr_alpha);
const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha, //
0, mr_sin_alpha};
const TinyMatrix mr_inv_W = inverse(mr_W);
if (print) {
std::cout << "MC === i_cell = " << i_cell << '\n';
}
{ // first simplex
const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
const TinyVector b = old_xr[cell_node_list[(i_cell_node + 2) % cell_nb_nodes]];
const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
const TinyMatrix mr_S = mr_A * mr_inv_W;
const double mr_sigma = det(mr_S);
if (mr_sigma < current_sigma_min) {
current_sigma_min = mr_sigma;
}
if (print) {
std::cout << "MC triangle 1 ===\n";
std::cout << "MC a = " << a << '\n';
std::cout << "MC b = " << b << '\n';
std::cout << "MC x = " << x << '\n';
}
}
{ // second simplex
const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], //
const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
S_list[i_cell] = A * inv_W;
const TinyMatrix mr_S = mr_A * mr_inv_W;
const double mr_sigma = det(mr_S);
if (mr_sigma < current_sigma_min) {
current_sigma_min = mr_sigma;
}
SmallArray<double> sigma_list(S_list.size());
for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
sigma_list[i_cell] = det(S_list[i_cell]);
if (print) {
std::cout << "MC triangle 2 ===\n";
std::cout << "MC a = " << a << '\n';
std::cout << "MC b = " << b << '\n';
std::cout << "MC x = " << x << '\n';
}
}
}
}
}
{
double current_sigma_min = min(sigma_list);
if (current_sigma_min < sigma_min) {
sigma_min = current_sigma_min;
delta = (sigma_min < eps)
? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
: 0;
}
}
for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
const double sigma = sigma_list[i_cell];
const TinyMatrix<Dimension> S = S_list[i_cell];
if (print) {
std::cout << " ==== \n";
std::cout << "delta = " << delta << '\n';
}
}
auto escobarSSigma = [&frobenius, delta, //
//
&print, &node_to_cell_matrix,
&node_id](const double sigma, const TinyMatrix<Dimension>& S,
const TinyVector<Dimension>& sigma_gradient,
const std::array<TinyMatrix<Dimension>, Dimension>& S_gradient,
const std::array<TinyMatrix<Dimension>, Dimension>& Sigma_gradient, double& f,
TinyVector<Dimension>& f_gradient, TinyMatrix<Dimension>& f_hessian) {
const TinyMatrix<Dimension> Sigma = [&S] {
TinyMatrix<Dimension> transposed_comS;
for (size_t i = 0; i < Dimension; ++i) {
......@@ -585,15 +670,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double f_jr = S_norm * Sigma_norm / h;
TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
trace(Sigma * S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> //
Sigma_gradient{TinyMatrix<Dimension>{0, inv_sin_alpha - inv_tan_alpha, //
0, -1},
TinyMatrix<Dimension>{-inv_sin_alpha + inv_tan_alpha, 0, //
1, 0}};
TinyVector<Dimension> g{trace(transpose(S) * S_gradient[0]) / S_norm2 //
+ trace(transpose(Sigma) * Sigma_gradient[0]) / Sigma_norm2 //
- trace(Sigma * S_gradient[0]) / (h - sigma),
......@@ -603,6 +679,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
- trace(Sigma * S_gradient[1]) / (h - sigma)};
const TinyVector<Dimension> f_jr_gradient = f_jr * g;
TinyMatrix<Dimension> f_jr_hessian = zero;
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
......@@ -615,20 +692,189 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
- 2 * trace(transpose(Sigma) * Sigma_gradient[j]) * trace(transpose(Sigma) * Sigma_gradient[i]) /
(Sigma_norm2 * Sigma_norm2) //
- 2 * trace(Sigma_gradient[j] * S_gradient[i]) / (h - sigma) //
+ 2 * trace(Sigma * S_gradient[i]) * sigma / (std::pow(h - sigma, 3)) * sigma_gradient[j] //
- trace(Sigma_gradient[j] * S_gradient[i]) / (h - sigma) //
+ sigma / (std::pow(h - sigma, 3)) * sigma_gradient[i] * sigma_gradient[j] //
+ g[j] * g[i]) *
f_jr;
}
}
if (print) {
std::cout << "S = " << S << '\n';
std::cout << "Sigma = " << Sigma << '\n';
std::cout << "normS = " << S_norm << '\n';
std::cout << "normSigma = " << Sigma_norm << '\n';
std::cout << "sigma = " << sigma << '\n';
std::cout << "f_jr = " << f_jr << '\n';
std::cout << "g = " << g << '\n';
std::cout << "grad_f_jr = " << f_jr_gradient << '\n';
std::cout << "hess_f_jr = " << f_jr_hessian << '\n';
}
f += f_jr;
f_gradient += f_jr_gradient;
f_hessian += f_jr_hessian;
};
auto escobar = escobarSSigma;
for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
if (print) {
std::cout << "------ cell = " << i_cell << '\n';
}
if (l2Norm(f_gradient) < 1E-12) {
break;
const double sigma = sigma_list[i_cell];
const TinyMatrix<Dimension> S = S_list[i_cell];
const TinyMatrix<Dimension> Sigma = [&S] {
TinyMatrix<Dimension> transposed_comS;
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
transposed_comS(i, j) = cofactor(S, j, i);
}
}
return transposed_comS;
}();
TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
trace(Sigma * S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> //
Sigma_gradient{TinyMatrix<Dimension>{+0, +1, //
+0, -1},
TinyMatrix<Dimension>{-inv_sin_alpha, -inv_tan_alpha, //
+inv_sin_alpha, +inv_tan_alpha}};
escobar(sigma, S, sigma_gradient, S_gradient, Sigma_gradient, f, f_gradient, f_hessian);
}
if constexpr (multi_connection) {
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]];
const size_t cell_nb_nodes = cell_node_list.size();
if (print) {
std::cout << "=== MC cell " << i_cell << " ===\n";
}
if (cell_nb_nodes <= 3) {
continue;
}
const size_t i_cell_node = node_number_in_cell[i_cell];
const double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
if (print) {
std::cout << "=== MC mr_alpha " << mr_alpha << " ===\n";
}
const double mr_sin_alpha = std::sin(mr_alpha);
const double mr_cos_alpha = std::cos(mr_alpha);
const double mr_inv_sin_alpha = 1. / mr_sin_alpha;
const double mr_inv_tan_alpha = 1. / std::tan(mr_alpha);
const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha, //
0, mr_sin_alpha};
const TinyMatrix mr_inv_W = inverse(mr_W);
std::array<TinyMatrix<Dimension>, Dimension> //
mr_S_gradient = {TinyMatrix<Dimension>{-1, (mr_cos_alpha - 1) / mr_sin_alpha, //
+0, +0},
TinyMatrix<Dimension>{+0, +0, //
-1, (mr_cos_alpha - 1) / mr_sin_alpha}};
{ // first simplex
const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
const TinyVector b = old_xr[cell_node_list[(i_cell_node + 2) % cell_nb_nodes]];
const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
const TinyMatrix mr_S = mr_A * mr_inv_W;
const double mr_sigma = det(mr_S);
const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
TinyMatrix<Dimension> transposed_comS;
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
transposed_comS(i, j) = cofactor(mr_S, j, i);
}
}
return transposed_comS;
}();
TinyVector<Dimension> mr_sigma_gradient{trace(mr_Sigma * mr_S_gradient[0]), //
trace(mr_Sigma * mr_S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> //
mr_Sigma_gradient{TinyMatrix<Dimension>{+0, (1 - mr_cos_alpha) / mr_sin_alpha, //
+0, -1},
TinyMatrix<Dimension>{(mr_cos_alpha - 1) / mr_sin_alpha, +0, //
+1, +0}};
if (print) {
std::cout << "=== MC triangle 1 ===\n";
std::cout << " A = " << mr_A << '\n';
std::cout << " W = " << mr_W << '\n';
std::cout << " inv(W) = " << mr_inv_W << '\n';
std::cout << " A * inv(W) = " << mr_A * mr_inv_W << '\n';
std::cout << " mr_S = " << mr_S << '\n';
std::cout << " mr_Sigma = " << mr_Sigma << '\n';
std::cout << " mr_S_gradient = " << mr_S_gradient[0] << " || " << mr_S_gradient[1] << '\n';
std::cout << " mr_Sigma_gradient = " << mr_Sigma_gradient[0] << " || " << mr_Sigma_gradient[1]
<< '\n';
}
escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
}
{ // second simplex
const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
const TinyMatrix mr_S = mr_A * mr_inv_W;
const double mr_sigma = det(mr_S);
const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
TinyMatrix<Dimension> transposed_comS;
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
transposed_comS(i, j) = cofactor(mr_S, j, i);
}
}
return transposed_comS;
}();
TinyVector<Dimension> mr_sigma_gradient{trace(mr_Sigma * mr_S_gradient[0]), //
trace(mr_Sigma * mr_S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> //
mr_Sigma_gradient{TinyMatrix<Dimension>{0, +1, //
0, -1},
TinyMatrix<Dimension>{-mr_inv_sin_alpha, -mr_inv_tan_alpha, //
+mr_inv_sin_alpha, +mr_inv_tan_alpha}};
if (print) {
std::cout << "=== MC triangle 2 ===\n";
std::cout << " A = " << mr_A << '\n';
std::cout << " W = " << mr_W << '\n';
std::cout << " inv(W) = " << mr_inv_W << '\n';
std::cout << " mr_S = " << mr_S << '\n';
std::cout << " mr_Sigma = " << mr_Sigma << '\n';
}
escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
}
}
}
const double max_edge_lenght = [=] {
......@@ -651,13 +897,25 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double delta_x_norm = l2Norm(delta_x_candidate);
double factor = 1;
if (delta_x_norm > 0.3 * max_edge_lenght) {
x += (0.3 * max_edge_lenght / delta_x_norm) * delta_x_candidate;
} else {
x += delta_x_candidate;
factor = (0.3 * max_edge_lenght / delta_x_norm);
}
x += factor * delta_x_candidate;
if (delta_x_norm < 1E-2 * max_edge_lenght) {
if (print) {
std::cout << "==============\n";
std::cout << "f_gradient = " << f_gradient << '\n';
std::cout << "f_hessian = " << f_hessian << '\n';
std::cout << "inv(f_hes) = " << inverse(f_hessian) << '\n';
std::cout << "factor = " << factor << '\n';
std::cout << "delta_x= " << delta_x_candidate << '\n';
}
if (l2Norm(f_gradient) < 1E-12) {
break;
}
if (delta_x_norm < 1E-8 * max_edge_lenght) {
break;
}
......@@ -667,12 +925,21 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return final_f;
};
parallel_for(
connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
// parallel_for(
// connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
// if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
// quality[node_id] = smooth(node_id, new_xr[node_id]);
// }
// });
size_t count = 0;
for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
quality[node_id] = smooth(node_id, new_xr[node_id]);
count++;
}
});
}
std::cout << "nb smoothed = " << count << '\n';
} else {
throw NotImplementedError("Dimension != 2");
}
......@@ -691,9 +958,13 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
NodeValue<Rd> new_xr = copy(xr);
this->_getNewPositions(xr, new_xr);
xr = new_xr;
nb_non_convex_cells = sum(this->_getNonConvexCellTag(xr));
auto non_convex_cell_tag = this->_getNonConvexCellTag(xr);
nb_non_convex_cells = sum(non_convex_cell_tag);
std::cout << " remaining non convex cells: " << nb_non_convex_cells << '\n';
} while ((nb_non_convex_cells > 0) and (i_convexify < 120));
if (nb_non_convex_cells > 0) {
m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
}
} while ((nb_non_convex_cells > 0) and (i_convexify < 10));
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
......@@ -784,9 +1055,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
is_displaced.fill(false);
CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
is_zone_cell.fill(false);
......@@ -809,6 +1077,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return m_p_given_mesh;
}
NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
is_displaced.fill(false);
for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
auto characteristic_function =
discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment