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

Remove algorithm legacy files

parent 9aba2979
Branches
No related tags found
No related merge requests found
Showing
with 1 addition and 4278 deletions
# ------------------- Source files -------------------- # ------------------- Source files --------------------
add_library(PugsLanguageAlgorithms add_library(PugsLanguageAlgorithms
ElasticityDiamondAlgorithm.cpp INTERFACE # THIS SHOULD BE REMOVED IF FILES ARE FINALY PROVIDED
UnsteadyElasticity.cpp
Heat5PointsAlgorithm.cpp
HeatDiamondAlgorithm.cpp
HeatDiamondAlgorithm2.cpp
ParabolicHeat.cpp
# INTERFACE # THIS SHOULD BE REMOVED IF FILES ARE FINALY PROVIDED
) )
......
This diff is collapsed.
#ifndef ELASTICITY_DIAMOND_ALGORITHM_HPP
#define ELASTICITY_DIAMOND_ALGORITHM_HPP
#include <algebra/TinyVector.hpp>
#include <language/utils/FunctionSymbolId.hpp>
#include <memory>
#include <mesh/IMesh.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <variant>
#include <vector>
template <size_t Dimension>
class ElasticityDiamondScheme
{
private:
class DirichletBoundaryCondition;
class NormalStrainBoundaryCondition;
class SymmetryBoundaryCondition;
public:
ElasticityDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& lambda_id,
const FunctionSymbolId& mu_id,
const FunctionSymbolId& f_id,
const FunctionSymbolId& U_id);
};
#endif // ELASTICITY_DIAMOND_ALGORITHM2_HPP
#include <language/algorithms/Heat5PointsAlgorithm.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/LeastSquareSolver.hpp>
#include <algebra/LinearSolver.hpp>
#include <algebra/LinearSolverOptions.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <algebra/Vector.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/DualMeshManager.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
template <size_t Dimension>
Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<Dimension>;
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
if constexpr (Dimension == 2) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
CellValue<double> Tj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
NodeValue<double> Tr(mesh->connectivity());
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
CellValuePerNode<double> w_rj{mesh->connectivity()};
for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
SmallVector<double> b{Dimension + 1};
b[0] = 1;
for (size_t i = 1; i < Dimension + 1; i++) {
b[i] = xr[i_node][i - 1];
}
const auto& node_to_cell = node_to_cell_matrix[i_node];
SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
for (size_t j = 0; j < node_to_cell.size(); j++) {
A(0, j) = 1;
}
for (size_t i = 1; i < Dimension + 1; i++) {
for (size_t j = 0; j < node_to_cell.size(); j++) {
const CellId J = node_to_cell[j];
A(i, j) = xj[J][i - 1];
}
}
SmallVector<double> x{node_to_cell.size()};
x = zero;
LeastSquareSolver ls_solver;
ls_solver.solveLocalSystem(A, x, b);
Tr[i_node] = 0;
for (size_t j = 0; j < node_to_cell.size(); j++) {
Tr[i_node] += x[j] * Tj[node_to_cell[j]];
w_rj(i_node, j) = x[j];
}
}
{
std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
std::shared_ptr mapper =
DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
NodeValue<double> Trd{diamond_mesh->connectivity()};
mapper->toDualNode(Tr, Tj, Trd);
CellValue<double> kappaj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
mesh_data.xj());
CellValue<double> dual_kappaj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
diamond_mesh_data
.xj());
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
FaceValue<double> compute_mes_l{mesh->connectivity()};
compute_mes_l.fill(1);
return compute_mes_l;
} else {
return mesh_data.ll();
}
}();
const CellValue<const double> dual_mes_l_j = [=] {
CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
mapper->toDualCell(mes_l, compute_mes_j);
return compute_mes_j;
}();
FaceValue<const double> alpha_l = [&] {
CellValue<double> alpha_j{diamond_mesh->connectivity()};
parallel_for(
diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
alpha_j[diamond_cell_id] =
dual_mes_l_j[diamond_cell_id] * dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
});
FaceValue<double> computed_alpha_l{mesh->connectivity()};
mapper->fromDualCell(alpha_j, computed_alpha_l);
return computed_alpha_l;
}();
const Array<int> non_zeros{mesh->numberOfCells()};
non_zeros.fill(Dimension);
CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
const CellId cell_id1 = primal_face_to_cell[i_cell];
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
const CellId cell_id2 = primal_face_to_cell[j_cell];
if (i_cell == j_cell) {
S(cell_id1, cell_id2) -= beta_l;
} else {
S(cell_id1, cell_id2) += beta_l;
}
}
}
}
CellValue<double> fj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
const CellValue<const double> primal_Vj = mesh_data.Vj();
CRSMatrix A{S.getCRSMatrix()};
Vector<double> b{mesh->numberOfCells()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
Vector<double> T{mesh->numberOfCells()};
T = zero;
LinearSolver solver;
solver.solveLocalSystem(A, T, b);
CellValue<double> Temperature{mesh->connectivity()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_id]; });
Vector<double> error{mesh->numberOfCells()};
parallel_for(
mesh->numberOfCells(),
PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * primal_Vj[cell_id]; });
std::cout << "||Error||_2 = " << std::sqrt(dot(error, error)) << "\n";
}
} else {
throw NotImplementedError("not done in this dimension");
}
}
template Heat5PointsAlgorithm<1>::Heat5PointsAlgorithm(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
template Heat5PointsAlgorithm<2>::Heat5PointsAlgorithm(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
template Heat5PointsAlgorithm<3>::Heat5PointsAlgorithm(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
#ifndef HEAT_5POINTS_ALGORITHM_HPP
#define HEAT_5POINTS_ALGORITHM_HPP
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <memory>
#include <vector>
template <size_t Dimension>
struct Heat5PointsAlgorithm
{
Heat5PointsAlgorithm(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id);
};
#endif // HEAT_5POINTS_ALGORITHM_HPP
This diff is collapsed.
#ifndef HEAT_DIAMOND_ALGORITHM_HPP
#define HEAT_DIAMOND_ALGORITHM_HPP
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <memory>
#include <variant>
#include <vector>
template <size_t Dimension>
class HeatDiamondScheme
{
private:
class DirichletBoundaryCondition;
class FourierBoundaryCondition;
class NeumannBoundaryCondition;
class SymmetryBoundaryCondition;
public:
HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
This diff is collapsed.
#ifndef HEAT_DIAMOND_ALGORITHM2_HPP
#define HEAT_DIAMOND_ALGORITHM2_HPP
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <memory>
#include <variant>
#include <vector>
template <size_t Dimension>
class HeatDiamondScheme2
{
private:
class DirichletBoundaryCondition;
class FourierBoundaryCondition;
class NeumannBoundaryCondition;
class SymmetryBoundaryCondition;
public:
HeatDiamondScheme2(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id);
};
#endif // HEAT_DIAMOND_ALGORITHM2_HPP
This diff is collapsed.
#ifndef PARABOLIC_HEAT_HPP
#define PARABOLIC_HEAT_HPP
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IMesh.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <memory>
#include <variant>
#include <vector>
template <size_t Dimension>
class ParabolicHeatScheme
{
private:
class DirichletBoundaryCondition;
class FourierBoundaryCondition;
class NeumannBoundaryCondition;
class SymmetryBoundaryCondition;
class InterpolationWeightsManager;
public:
ParabolicHeatScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& T_init_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
const double& final_time,
const double& dt);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
This diff is collapsed.
#ifndef UNSTEADY_ELASTICITY_HPP
#define UNSTEADY_ELASTICITY_HPP
#include <algebra/TinyVector.hpp>
#include <language/utils/FunctionSymbolId.hpp>
#include <memory>
#include <mesh/IMesh.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <variant>
#include <vector>
template <size_t Dimension>
class UnsteadyElasticity
{
private:
class DirichletBoundaryCondition;
class NormalStrainBoundaryCondition;
class SymmetryBoundaryCondition;
class InterpolationWeightsManager;
public:
UnsteadyElasticity(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& lambda_id,
const FunctionSymbolId& mu_id,
const FunctionSymbolId& f_id,
const FunctionSymbolId& U_id,
const FunctionSymbolId& U_init_id,
const double& Tf,
const double& dt);
};
#endif // ELASTICITY_DIAMOND_ALGORITHM2_HPP
...@@ -3,12 +3,6 @@ ...@@ -3,12 +3,6 @@
#include <analysis/GaussLegendreQuadratureDescriptor.hpp> #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussLobattoQuadratureDescriptor.hpp> #include <analysis/GaussLobattoQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp> #include <analysis/GaussQuadratureDescriptor.hpp>
#include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
#include <language/algorithms/Heat5PointsAlgorithm.hpp>
#include <language/algorithms/HeatDiamondAlgorithm.hpp>
#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
#include <language/algorithms/ParabolicHeat.hpp>
#include <language/algorithms/UnsteadyElasticity.hpp>
#include <language/modules/BinaryOperatorRegisterForVh.hpp> #include <language/modules/BinaryOperatorRegisterForVh.hpp>
#include <language/modules/MathFunctionRegisterForVh.hpp> #include <language/modules/MathFunctionRegisterForVh.hpp>
#include <language/modules/UnaryOperatorRegisterForVh.hpp> #include <language/modules/UnaryOperatorRegisterForVh.hpp>
...@@ -440,67 +434,6 @@ SchemeModule::SchemeModule() ...@@ -440,67 +434,6 @@ SchemeModule::SchemeModule()
)); ));
this->_addBuiltinFunction("heat", std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 2: {
HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 3: {
HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction("parabolicheat",
std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& T_init_id,
const FunctionSymbolId& kappa_id, const FunctionSymbolId& f_id,
const double& final_time, const double& dt) -> void {
switch (p_mesh->dimension()) {
case 1: {
ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt};
break;
}
case 2: {
ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt};
break;
}
case 3: {
ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction( this->_addBuiltinFunction(
"parabolicheat", "parabolicheat",
std::function( std::function(
...@@ -515,40 +448,6 @@ SchemeModule::SchemeModule() ...@@ -515,40 +448,6 @@ SchemeModule::SchemeModule()
)); ));
this->_addBuiltinFunction("unsteadyelasticity",
std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
const FunctionSymbolId& f_id, const FunctionSymbolId& U_id,
const FunctionSymbolId& U_init_id, const double& final_time,
const double& dt) -> void {
switch (p_mesh->dimension()) {
case 1: {
UnsteadyElasticity<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
U_id, U_init_id, final_time, dt};
break;
}
case 2: {
UnsteadyElasticity<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
U_id, U_init_id, final_time, dt};
break;
}
case 3: {
UnsteadyElasticity<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
U_id, U_init_id, final_time, dt};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction("unsteadyelasticity", this->_addBuiltinFunction("unsteadyelasticity",
std::function( std::function(
...@@ -603,92 +502,6 @@ SchemeModule::SchemeModule() ...@@ -603,92 +502,6 @@ SchemeModule::SchemeModule()
)); ));
this->_addBuiltinFunction("heat2", std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
HeatDiamondScheme2<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 2: {
HeatDiamondScheme2<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 3: {
HeatDiamondScheme2<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction("elasticity",
std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
ElasticityDiamondScheme<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
case 2: {
ElasticityDiamondScheme<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
case 3: {
ElasticityDiamondScheme<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction("heat5points",
std::function(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
Heat5PointsAlgorithm<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 2: {
Heat5PointsAlgorithm<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
case 3: {
Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
default: {
throw UnexpectedError("invalid mesh dimension");
}
}
}
));
this->_addBuiltinFunction("lagrangian", this->_addBuiltinFunction("lagrangian",
std::function( std::function(
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment