Select Git revision
MeshSmootherEscobar.cpp
MeshSmootherEscobar.cpp 54.23 KiB
#include <mesh/MeshSmootherEscobar.hpp>
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <geometry/TriangleTransformation.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/ItemValueVariant.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshCellZone.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshLineNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/RandomEngine.hpp>
#include <variant>
template <size_t Dimension>
class MeshSmootherEscobarHandler::MeshSmootherEscobar
{
private:
using Rd = TinyVector<Dimension>;
using Rdxd = TinyMatrix<Dimension>;
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
const MeshType& m_given_mesh;
std::shared_ptr<const IMesh> m_p_given_mesh;
class AxisBoundaryCondition;
class FixedBoundaryCondition;
class SymmetryBoundaryCondition;
using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList m_boundary_condition_list;
NodeValue<const bool> m_is_fixed_node;
NodeValue<const bool> m_is_smoothed_node;
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
BoundaryConditionList bc_list;
for (const auto& bc_descriptor : bc_descriptor_list) {
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::axis: {
if constexpr (Dimension == 1) {
bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
} else {
bc_list.emplace_back(
AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
}
break;
}
case IBoundaryConditionDescriptor::Type::symmetry: {
bc_list.emplace_back(
SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
break;
}
case IBoundaryConditionDescriptor::Type::fixed: {
bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
break;
}
default: {
std::ostringstream error_msg;
error_msg << *bc_descriptor << " is an invalid boundary condition for mesh smoother";
throw NormalError(error_msg.str());
}
}
}
return bc_list;
}
NodeValue<const bool>
_getFixedNodes() const
{
NodeValue<bool> is_fixed{m_given_mesh.connectivity()};
is_fixed.fill(false);
for (auto&& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr ((std::is_same_v<BCType, FixedBoundaryCondition>) or
((Dimension == 2) and std::is_same_v<BCType, AxisBoundaryCondition>)) {
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) { is_fixed[node_list[i_node]] = true; });
}
},
boundary_condition);
}
if constexpr (Dimension == 3) {
throw NotImplementedError("treat sliding axis node kind");
}
synchronize(is_fixed);
NodeValue<int> bc_number{m_given_mesh.connectivity()};
bc_number.fill(-1);
{
int bc_id = 0;
for (auto&& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
if (not is_fixed[node_id]) {
if (bc_number[node_id] < 0) {
bc_number[node_id] = bc_id;
} else {
bc_number[node_id] = -2; // error tag
}
}
});
}
},
boundary_condition);
bc_id++;
}
if (min(bc_number) < -1) {
throw NormalError("Smoothing boundary conditions overlap. One must precise fixed points and smoothing axis");
}
}
bool missing_bc = false;
auto is_boundary = m_given_mesh.connectivity().isBoundaryNode();
parallel_for(m_given_mesh.numberOfNodes(), [&](NodeId node_id) {
if (is_boundary[node_id]) {
if (not is_fixed[node_id] and (bc_number[node_id] == -1)) {
missing_bc = true;
}
}
});
if (parallel::allReduceOr(missing_bc)) {
throw NormalError("Some boundary conditions are missing (boundary nodes without BC)");
}
return is_fixed;
}
void
_applyBC(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
{
for (auto&& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
const Rd& n = bc.outgoingNormal();
const auto node_list = bc.nodeList();
const ConnectivityType& connectivity = m_given_mesh.connectivity();
auto is_boundary_node = connectivity.isBoundaryNode();
auto node_is_owned = connectivity.nodeIsOwned();
if constexpr (Dimension == 2) {
const TinyVector<Dimension> t{-n[1], n[0]};
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
constexpr double eps = 1E-14;
auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
auto cell_list = node_to_cell_matrix[node_id];
auto node_number_in_cell = node_number_in_their_cells[node_id];
const double alpha = std::acos(-1) / cell_list.size();
const double cos_alpha = std::cos(alpha);
const double sin_alpha = std::sin(alpha);
const double inv_sin_alpha = 1. / sin_alpha;
const double inv_tan_alpha = 1. / std::tan(alpha);
const TinyMatrix<Dimension> W{1, cos_alpha, //
0, sin_alpha};
const TinyMatrix<Dimension> inv_W = inverse(W);
TinyMatrix<Dimension, Dimension> dtheta_S =
TinyMatrix<Dimension>{t[0], t[0] * (inv_sin_alpha - inv_tan_alpha), //
t[1], t[1] * (inv_sin_alpha - inv_tan_alpha)};
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const size_t i_cell_node = node_number_in_cell[i_cell];
const auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]];
const size_t cell_nb_nodes = cell_node_list.size();
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 + cell_nb_nodes - 1) % cell_nb_nodes]];
const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
S_list[i_cell] = A * inv_W;
}
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]);
}
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 final_f = 0;
for (size_t i_iter = 0; i_iter < 200; ++i_iter) {
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
for (size_t i_cell = 0; i_cell < cell_list.size(); ++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 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 + cell_nb_nodes - 1) % cell_nb_nodes]];
const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], //
a[1] - x[1], b[1] - x[1]};
S_list[i_cell] = A * inv_W;
}
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]);
}
{
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;
}
}
double f = 0;
double dtheta_f = 0;
double d2theta_f = 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];
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;
}();
const double S_norm = frobenius(S);
const double Sigma_norm = frobenius(Sigma);
const double S_norm2 = S_norm * S_norm;
const double Sigma_norm2 = Sigma_norm * Sigma_norm;
const double h = sigma + std::sqrt(sigma * sigma + 4 * delta * delta);
const double f_jr = S_norm * Sigma_norm / h;
const double dtheta_sigma = trace(Sigma * dtheta_S);
TinyMatrix<Dimension, Dimension> dtheta_Sigma =
TinyMatrix<Dimension>{t[1] * (inv_sin_alpha - inv_tan_alpha),
-t[0] * (inv_sin_alpha - inv_tan_alpha), //
-t[1], t[0]};
const double g = trace(transpose(S) * dtheta_S) / S_norm2 //
+ trace(transpose(Sigma) * dtheta_Sigma) / Sigma_norm2 //
- dtheta_sigma / (h - sigma);
const double dtheta_f_jr = f_jr * g;
const double dtheta2_f_jr =
(trace(transpose(dtheta_S) * dtheta_S) / S_norm2 //
- 2 * trace(transpose(S) * dtheta_S) * trace(transpose(S) * dtheta_S) / (S_norm2 * S_norm2) //
//
+ trace(transpose(dtheta_Sigma) * dtheta_Sigma) / Sigma_norm2 // + 0
- 2 * trace(transpose(Sigma) * dtheta_Sigma) * trace(transpose(Sigma) * dtheta_Sigma) /
(Sigma_norm2 * Sigma_norm2) //
//
- 2 * trace(dtheta_Sigma * dtheta_S) / (h - sigma) //
+ 2 * sigma / (std::pow(h - sigma, 3)) * dtheta_sigma * dtheta_sigma //
+ g * g) *
f_jr;
f += f_jr;
dtheta_f += dtheta_f_jr;
d2theta_f += dtheta2_f_jr;
}
if (std::abs(dtheta_f) < 1E-12) {
break;
}
const double max_edge_lenght = [=] {
auto node_edge_list = node_to_edge_matrix[node_id];
double max_lenght = 0;
for (size_t i_edge = 0; i_edge < node_edge_list.size(); ++i_edge) {
const EdgeId edge_id = node_edge_list[i_edge];
auto edge_node_list = edge_to_node_matrix[edge_id];
const TinyVector<Dimension>& x0 = old_xr[edge_node_list[0]];
const TinyVector<Dimension>& x1 = old_xr[edge_node_list[1]];
const double lenght = l2Norm(x0 - x1);
if (lenght > max_lenght) {
max_lenght = lenght;
}
}
return max_lenght;
}();
const double delta_x_candidate = dtheta_f / std::abs(d2theta_f);
const double delta_x_norm = std::abs(delta_x_candidate);
if (delta_x_norm > 0.3 * max_edge_lenght) {
x += (0.3 * max_edge_lenght / delta_x_norm) * delta_x_candidate * t;
} else {
x += delta_x_candidate * t;
}
if (delta_x_norm < 1E-2 * max_edge_lenght) {
break;
}
final_f = f;
}
return final_f;
};
parallel_for(
node_list.size(), PUGS_LAMBDA(const size_t i_node) {
const NodeId node_id = node_list[i_node];
if (m_is_smoothed_node[node_id] and node_is_owned[node_id]) {
smooth(node_id, new_xr[node_id]);
}
});
} else {
throw NotImplementedError("Dimension != 2");
}
} else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
if constexpr (Dimension > 1) {
throw NotImplementedError("Escobar: axis boundary conditions");
} else {
throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
}
} else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(const size_t i_node) {
const NodeId node_id = node_list[i_node];
new_xr[node_id] = old_xr[node_id];
});
} else {
throw UnexpectedError("invalid boundary condition type");
}
},
boundary_condition);
}
}
CellValue<const int>
_getNonConvexCellTag(const NodeValue<const Rd>& xr) const
{
CellValue<int> non_convex_cell_tag{m_given_mesh.connectivity()};
auto cell_to_node_matrix = m_given_mesh.connectivity().cellToNodeMatrix();
if constexpr (Dimension == 2) {
parallel_for(
m_given_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto cell_node_list = cell_to_node_matrix[cell_id];
bool is_convex = true;
for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
const TriangleTransformation<Dimension> T{xr[cell_node_list[i_node]],
xr[cell_node_list[(i_node + 1) % cell_node_list.size()]],
xr[cell_node_list[(i_node + 2) % cell_node_list.size()]]};
if (T.jacobianDeterminant() <= 0) {
is_convex = false;
break;
}
}
non_convex_cell_tag[cell_id] = (is_convex == false);
});
} else {
throw NotImplementedError("non convex cell calculation in dimension != 2");
}
return non_convex_cell_tag;
}
NodeValue<const bool>
_getSmoothedNode(CellValue<const int>& non_convex_cell_tag, const NodeValue<const bool>& m_is_fixed_node) const
{
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
NodeValue<bool> is_smoothed{m_given_mesh.connectivity()};
is_smoothed.fill(false);
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
if (not m_is_fixed_node[node_id]) {
auto node_cell_list = node_to_cell_matrix[node_id];
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
const CellId cell_id = node_cell_list[i_cell];
if (non_convex_cell_tag[cell_id]) {
is_smoothed[node_id] = true;
break;
}
}
}
});
return is_smoothed;
}
public:
std::shared_ptr<const ItemValueVariant>
getQuality() const
{
return std::make_shared<ItemValueVariant>(m_given_mesh.xr());
}
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();
auto node_is_owned = connectivity.nodeIsOwned();
if constexpr (Dimension == 2) {
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
NodeValue<double> quality{connectivity};
constexpr double eps = 1E-14;
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];
const double alpha = 2 * std::acos(-1) / cell_list.size();
const double sin_alpha = std::sin(alpha);
const double cos_alpha = std::cos(alpha);
const double inv_sin_alpha = 1. / sin_alpha;
const double inv_tan_alpha = 1. / std::tan(alpha);
const TinyMatrix<Dimension> W{1, cos_alpha, //
0, sin_alpha};
const TinyMatrix<Dimension> inv_W = inverse(W);
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) {
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 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 + cell_nb_nodes - 1) % cell_nb_nodes]];
const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], //
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());
for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
sigma_list[i_cell] = det(S_list[i_cell]);
}
{
double current_sigma_min = min(sigma_list);
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);
}
if (nb_medium_range_simplices > 0) {
SmallArray<double> medium_range_sigma_list(nb_medium_range_simplices);
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> 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 2 ===\n";
std::cout << "MC a = " << a << '\n';
std::cout << "MC b = " << b << '\n';
std::cout << "MC x = " << x << '\n';
}
}
}
}
}
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;
}
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) {
for (size_t j = 0; j < Dimension; ++j) {
transposed_comS(i, j) = cofactor(S, j, i);
}
}
return transposed_comS;
}();
const double S_norm = frobenius(S);
const double Sigma_norm = frobenius(Sigma);
const double S_norm2 = S_norm * S_norm;
const double Sigma_norm2 = Sigma_norm * Sigma_norm;
double h = sigma + std::sqrt(sigma * sigma + 4 * delta * delta);
const double f_jr = S_norm * Sigma_norm / h;
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),
//
trace(transpose(S) * S_gradient[1]) / S_norm2 //
+ trace(transpose(Sigma) * Sigma_gradient[1]) / Sigma_norm2 //
- 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) {
f_jr_hessian(i, j) = //
(trace(transpose(S_gradient[j]) * S_gradient[i]) / S_norm2 //
- 2 * trace(transpose(S) * S_gradient[j]) * trace(transpose(S) * S_gradient[i]) /
(S_norm2 * S_norm2) //
//
+ trace(transpose(Sigma_gradient[j]) * Sigma_gradient[i]) / Sigma_norm2 // + 0
- 2 * trace(transpose(Sigma) * Sigma_gradient[j]) * trace(transpose(Sigma) * Sigma_gradient[i]) /
(Sigma_norm2 * Sigma_norm2) //
- 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';
}
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 = [=] {
auto node_edge_list = node_to_edge_matrix[node_id];
double max_lenght = 0;
for (size_t i_edge = 0; i_edge < node_edge_list.size(); ++i_edge) {
const EdgeId edge_id = node_edge_list[i_edge];
auto edge_node_list = edge_to_node_matrix[edge_id];
const TinyVector<Dimension>& x0 = old_xr[edge_node_list[0]];
const TinyVector<Dimension>& x1 = old_xr[edge_node_list[1]];
const double lenght = l2Norm(x0 - x1);
if (lenght > max_lenght) {
max_lenght = lenght;
}
}
return max_lenght;
}();
const TinyVector delta_x_candidate = -inverse(f_hessian) * f_gradient;
const double delta_x_norm = l2Norm(delta_x_candidate);
double factor = 1;
if (delta_x_norm > 0.3 * max_edge_lenght) {
factor = (0.3 * max_edge_lenght / delta_x_norm);
}
x += factor * delta_x_candidate;
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;
}
final_f = f;
}
return final_f;
};
// 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");
}
this->_applyBC(old_xr, new_xr);
}
std::shared_ptr<const IMesh>
_smoothedMesh()
{
NodeValue<const Rd> xr = m_given_mesh.xr();
size_t nb_non_convex_cells = 0;
size_t i_convexify = 1;
do {
std::cout << "- convexfication step " << i_convexify++ << '\n';
NodeValue<Rd> new_xr = copy(xr);
this->_getNewPositions(xr, new_xr);
xr = new_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';
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);
}
std::shared_ptr<const IMesh>
getSmoothedMesh()
{
NodeValue<bool> is_smoothed_node(m_given_mesh.connectivity());
parallel_for(
m_given_mesh.numberOfNodes(),
PUGS_LAMBDA(const NodeId node_id) { is_smoothed_node[node_id] = not m_is_fixed_node[node_id]; });
m_is_smoothed_node = is_smoothed_node;
return _smoothedMesh();
}
std::shared_ptr<const IMesh>
getSmoothedMesh(const FunctionSymbolId& function_symbol_id)
{
NodeValue<const Rd> given_xr = m_given_mesh.xr();
NodeValue<const bool> is_displaced =
InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
NodeValue<bool> is_fixed_node(m_given_mesh.connectivity());
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
});
synchronize(is_fixed_node);
m_is_fixed_node = is_fixed_node;
return _smoothedMesh();
}
std::shared_ptr<const IMesh>
getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list)
{
NodeValue<const Rd> given_xr = m_given_mesh.xr();
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
is_displaced.fill(false);
for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
const auto cell_list = cell_zone.cellList();
CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
is_zone_cell.fill(false);
parallel_for(
cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
auto node_cell_list = node_to_cell_matrix[node_id];
bool displace = true;
for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
const CellId cell_id = node_cell_list[i_node_cell];
displace &= is_zone_cell[cell_id];
}
if (displace) {
is_displaced[node_id] = true;
}
});
}
synchronize(is_displaced);
NodeValue<bool> is_fixed_node(m_given_mesh.connectivity());
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
});
synchronize(is_fixed_node);
m_is_fixed_node = is_fixed_node;
return _smoothedMesh();
}
std::shared_ptr<const IMesh>
getSmoothedMesh(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
{
NodeValue<const Rd> given_xr = m_given_mesh.xr();
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
is_zone_cell.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>>();
parallel_for(
m_given_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
if (characteristic_function[cell_id] != 0) {
is_zone_cell[cell_id] = true;
}
});
}
CellValue<const int> non_convex_cell_tag = this->_getNonConvexCellTag(given_xr);
const bool has_non_convex_cell = max(non_convex_cell_tag) > 0;
if (not has_non_convex_cell) {
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>>();
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
auto node_cell_list = node_to_cell_matrix[node_id];
bool displace = true;
for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
const CellId cell_id = node_cell_list[i_node_cell];
displace &= (characteristic_function[cell_id] != 0);
}
if (displace) {
is_displaced[node_id] = true;
}
});
}
synchronize(is_displaced);
NodeValue<bool> is_fixed_node(m_given_mesh.connectivity());
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
});
synchronize(is_fixed_node);
m_is_fixed_node = is_fixed_node;
m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
return _smoothedMesh();
}
MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
MeshSmootherEscobar(MeshSmootherEscobar&&) = delete;
MeshSmootherEscobar(const MeshType& given_mesh,
const std::shared_ptr<const IMesh>& p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
: m_given_mesh{given_mesh},
m_p_given_mesh{p_mesh},
m_boundary_condition_list{this->_getBCList(given_mesh, bc_descriptor_list)},
m_is_fixed_node{this->_getFixedNodes()}
{}
~MeshSmootherEscobar() = default;
};
template <size_t Dimension>
class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::AxisBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
public:
const Rd&
direction() const
{
return m_mesh_line_node_boundary.direction();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_line_node_boundary.nodeList();
}
AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
: m_mesh_line_node_boundary(mesh_line_node_boundary)
{
;
}
~AxisBoundaryCondition() = default;
};
template <>
class MeshSmootherEscobarHandler::MeshSmootherEscobar<1>::AxisBoundaryCondition
{
public:
AxisBoundaryCondition() = default;
~AxisBoundaryCondition() = default;
};
template <size_t Dimension>
class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::FixedBoundaryCondition
{
private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
~FixedBoundaryCondition() = default;
};
template <size_t Dimension>
class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
std::shared_ptr<const ItemValueVariant>
MeshSmootherEscobarHandler::getQuality(
const std::shared_ptr<const IMesh>& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getQuality();
}
case 2: {
constexpr size_t Dimension = 2;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getQuality();
}
case 3: {
constexpr size_t Dimension = 3;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getQuality();
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
std::shared_ptr<const IMesh>
MeshSmootherEscobarHandler::getSmoothedMesh(
const std::shared_ptr<const IMesh>& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh();
}
case 2: {
constexpr size_t Dimension = 2;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh();
}
case 3: {
constexpr size_t Dimension = 3;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh();
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
std::shared_ptr<const IMesh>
MeshSmootherEscobarHandler::getSmoothedMesh(
const std::shared_ptr<const IMesh>& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& function_symbol_id) const
{
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(function_symbol_id);
}
case 2: {
constexpr size_t Dimension = 2;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(function_symbol_id);
}
case 3: {
constexpr size_t Dimension = 3;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(function_symbol_id);
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
std::shared_ptr<const IMesh>
MeshSmootherEscobarHandler::getSmoothedMesh(
const std::shared_ptr<const IMesh>& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
{
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(smoothing_zone_list);
}
case 2: {
constexpr size_t Dimension = 2;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(smoothing_zone_list);
}
case 3: {
constexpr size_t Dimension = 3;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(smoothing_zone_list);
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
std::shared_ptr<const IMesh>
MeshSmootherEscobarHandler::getSmoothedMesh(
const std::shared_ptr<const IMesh>& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
if (not hasSameMesh(discrete_function_variant_list)) {
throw NormalError("discrete functions are not defined on the same mesh");
}
std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
if (common_mesh != mesh) {
throw NormalError("discrete functions are not defined on the smoothed mesh");
}
switch (mesh->dimension()) {
case 1: {
constexpr size_t Dimension = 1;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(discrete_function_variant_list);
}
case 2: {
constexpr size_t Dimension = 2;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(discrete_function_variant_list);
}
case 3: {
constexpr size_t Dimension = 3;
using MeshType = Mesh<Connectivity<Dimension>>;
MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
return smoother.getSmoothedMesh(discrete_function_variant_list);
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}