"src/language/utils/DataVariant.hpp" did not exist on "a446e8905b30e16152ca4ae8db14c6116ae29ed8"
Select Git revision
MeshSmootherEscobar.cpp 28.93 KiB
#include <mesh/MeshSmootherEscobar.hpp>
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.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;
class AxisBoundaryCondition;
class FixedBoundaryCondition;
class SymmetryBoundaryCondition;
using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList m_boundary_condition_list;
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
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;
}
void
_applyBC(NodeValue<Rd>& shift) 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 Rdxd I = identity;
const Rdxd nxn = tensorProduct(n, n);
const Rdxd P = I - nxn;
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];
shift[node_id] = P * shift[node_id];
});
} else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
if constexpr (Dimension > 1) {
const Rd& t = bc.direction();
const Rdxd txt = tensorProduct(t, t);
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];
shift[node_id] = txt * shift[node_id];
});
} 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];
shift[node_id] = zero;
});
} else {
throw UnexpectedError("invalid boundary condition type");
}
},
boundary_condition);
}
}
NodeValue<Rd>
_getDisplacement() const
{
const ConnectivityType& connectivity = m_given_mesh.connectivity();
NodeValue<const Rd> given_xr = m_given_mesh.xr();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
NodeValue<double> max_delta_xr{connectivity};
parallel_for(
connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
const Rd& x0 = given_xr[node_id];
const auto& node_cell_list = node_to_cell_matrix[node_id];
double min_distance_2 = std::numeric_limits<double>::max();
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
const CellId cell_id = node_cell_list[i_cell];
const auto& cell_node_list = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
if (i_node != i_cell_node) {
const NodeId cell_node_id = cell_node_list[i_node];
const Rd delta = x0 - given_xr[cell_node_id];
min_distance_2 = std::min(min_distance_2, dot(delta, delta));
}
}
}
double max_delta = std::sqrt(min_distance_2);
max_delta_xr[node_id] = max_delta;
});
NodeValue<Rd> shift_r{connectivity};
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
const auto& node_cell_list = node_to_cell_matrix[node_id];
Rd mean_position(zero);
size_t number_of_neighbours = 0;
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
const CellId cell_id = node_cell_list[i_cell];
const auto& cell_node_list = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
if (i_node != i_cell_node) {
const NodeId cell_node_id = cell_node_list[i_node];
mean_position += given_xr[cell_node_id];
number_of_neighbours++;
}
}
}
mean_position = 1. / number_of_neighbours * mean_position;
shift_r[node_id] = mean_position - given_xr[node_id];
});
this->_applyBC(shift_r);
synchronize(shift_r);
return shift_r;
}
public:
std::shared_ptr<const ItemValueVariant>
getQuality() const
{
if constexpr (Dimension == 2) {
const ConnectivityType& connectivity = m_given_mesh.connectivity();
NodeValue<const Rd> xr = m_given_mesh.xr();
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
auto is_boundary_node = connectivity.isBoundaryNode();
NodeValue<double> quality{connectivity};
constexpr double eps = 1E-15;
quality.fill(2);
auto f_inner = [=](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 = 2 * std::acos(-1) / cell_list.size();
const TinyMatrix<Dimension> W{1, std::cos(alpha), //
0, std::sin(alpha)};
const TinyMatrix<Dimension> inv_W = inverse(W);
std::array<TinyMatrix<Dimension>, Dimension> S_gradient =
{TinyMatrix<Dimension>{-1, -1. / std::sin(alpha) + 1. / std::tan(alpha), //
+0, +0}, //
TinyMatrix<Dimension>{+0, +0, //
-1, -1. / std::sin(alpha) + 1. / std::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];
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 = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
const TinyVector b = 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]);
}
const double sigma_min = min(sigma_list);
const 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)); };
// TinyVector<Dimension> f_gradient = zero;
// TinyMatrix<Dimension> f_hessian = zero;
double final_f = 0;
for (size_t i_iter = 0; i_iter < 100; ++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 = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
const TinyVector b = 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 f = 0;
TinyVector<Dimension> f_gradient = zero;
TinyMatrix<Dimension> f_hessian = zero;
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 = sigma * inverse(S);
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;
TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
trace(Sigma * S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> //
Sigma_gradient_old{sigma_gradient[0] * inverse(S) - inverse(S) * S_gradient[0] * Sigma,
sigma_gradient[1] * inverse(S) - inverse(S) * S_gradient[1] * Sigma};
const std::array<TinyMatrix<Dimension>, Dimension> //
Sigma_gradient_new{TinyMatrix<Dimension>{0, 1. / std::sin(alpha - 1. / std::tan(alpha)), //
0, -1},
TinyMatrix<Dimension>{-1. / std::sin(alpha) + 1. / std::tan(alpha), 0, //
1, 0}};
const auto Sigma_gradient = Sigma_gradient_new;
std::cout << "Sigma_gradient_old[0] = " << Sigma_gradient_old[0] << '\n';
std::cout << "Sigma_gradient_new[0] = " << Sigma_gradient_new[0] << '\n';
std::cout << "Sigma_gradient_old[1] = " << Sigma_gradient_old[1] << '\n';
std::cout << "Sigma_gradient_new[1] = " << Sigma_gradient_new[1] << '\n';
// TinyVector<Dimension> h_gradient = h / (h - sigma_list[i_cell]) * sigma_gradient;
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) //
//
- 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] //
+ g[j] * g[i]) *
f_jr;
}
}
f += f_jr;
f_gradient += f_jr_gradient;
f_hessian += f_jr_hessian;
}
std::cout << "f = " << f << '\n';
std::cout << "grad(f) = " << f_gradient << '\n';
std::cout << "hess(f) = " << f_hessian << " | hess(f)^T -hess(f) = " << transpose(f_hessian) - f_hessian
<< '\n';
std::cout << "inv(H) = " << inverse(f_hessian) << '\n';
std::cout << "inv(H)*grad(f) = " << inverse(f_hessian) * f_gradient << '\n';
std::cout << rang::fgB::yellow << "x = " << x << " -> " << x - inverse(f_hessian) * f_gradient
<< rang::fg::reset << '\n';
std::cout << rang::fgB::green << i_iter << ": l2Norm(f_gradient) = " << l2Norm(f_gradient) << rang::fg::reset
<< '\n';
if (l2Norm(f_gradient) < 1E-6) {
break;
}
x -= inverse(f_hessian) * f_gradient;
final_f = f;
}
return final_f;
};
parallel_for(
connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
// auto cell_list = node_to_cell_matrix[node_id];
// auto node_number_in_cell = node_number_in_their_cells[node_id];
if (is_boundary_node[node_id]) {
quality[node_id] = 1;
} else {
TinyVector x = xr[node_id];
quality[node_id] = f_inner(node_id, x);
std::exit(0);
// TinyMatrix<Dimension> B = identity;
}
});
return std::make_shared<ItemValueVariant>(quality);
} else {
throw NotImplementedError("Dimension != 2");
}
}
std::shared_ptr<const IMesh>
getSmoothedMesh() const
{
NodeValue<const Rd> given_xr = m_given_mesh.xr();
NodeValue<Rd> xr = this->_getDisplacement();
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { xr[node_id] += given_xr[node_id]; });
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
std::shared_ptr<const IMesh>
getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
{
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<Rd> xr = this->_getDisplacement();
parallel_for(
m_given_mesh.numberOfNodes(),
PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
std::shared_ptr<const IMesh>
getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const
{
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<Rd> xr = this->_getDisplacement();
parallel_for(
m_given_mesh.numberOfNodes(),
PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
std::shared_ptr<const IMesh>
getSmoothedMesh(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
{
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 < discrete_function_variant_list.size(); ++i_zone) {
auto is_zone_cell = 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 &= (is_zone_cell[cell_id] != 0);
}
if (displace) {
is_displaced[node_id] = true;
}
});
}
synchronize(is_displaced);
NodeValue<Rd> xr = this->_getDisplacement();
parallel_for(
m_given_mesh.numberOfNodes(),
PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
MeshSmootherEscobar(MeshSmootherEscobar&&) = delete;
MeshSmootherEscobar(const MeshType& given_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
: m_given_mesh(given_mesh), m_boundary_condition_list(this->_getBCList(given_mesh, bc_descriptor_list))
{}
~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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), bc_descriptor_list);
return smoother.getSmoothedMesh(discrete_function_variant_list);
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}