Skip to content
Snippets Groups Projects
Commit 4267419a authored by PATELA Julie's avatar PATELA Julie
Browse files

Add internal node's matrix to HeatDiamondAlgorithm; Pass second member in parameter

parent d67e7b1b
No related branches found
No related tags found
No related merge requests found
......@@ -21,7 +21,8 @@ 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& kappa_id,
const FunctionSymbolId& f_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
......@@ -176,24 +177,18 @@ Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
}
}
}
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};
Vector<double> b{mesh->numberOfCells()};
double pi = 3.14159265;
for (CellId i_cell = 0; i_cell < mesh->numberOfCells(); ++i_cell) {
if constexpr (Dimension == 2) {
b[i_cell] =
-2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
} else if constexpr (Dimension == 3) {
b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
}
}
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 = 0;
PCG{b, A, A, T, 100, 1e-12};
PCG{b, A, A, T, 1000, 1e-12};
CellValue<double> Temperature{mesh->connectivity()};
......@@ -202,7 +197,8 @@ Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
Vector<double> error{mesh->numberOfCells()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
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((error, error)) << "\n";
......@@ -222,16 +218,19 @@ 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&);
......@@ -14,7 +14,8 @@ 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& kappa_id,
const FunctionSymbolId& f_id);
};
#endif // HEAT_5POINTS_ALGORITHM_HPP
......@@ -21,7 +21,8 @@ HeatDiamondScheme<Dimension>::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& kappa_id,
const FunctionSymbolId& f_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
......@@ -115,26 +116,15 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
0, true); // forces last output
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
const FaceValue<const double> mes_l = [=] {
FaceValue<double> compute_mes_l{m_mesh->connectivity()};
const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
FaceValue<double> compute_mes_l{m_mesh->connectivity()};
compute_mes_l.fill(1);
} else if constexpr (Dimension == 2) {
parallel_for(
m_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& face_to_node = face_to_node_matrix[face_id];
const NodeId node_id1 = face_to_node[0];
const NodeId node_id2 = face_to_node[1];
const TinyVector<Dimension, double> r = xr[node_id1] - xr[node_id2];
compute_mes_l[face_id] = l2Norm(r);
});
return compute_mes_l;
} else {
throw NotImplementedError("Not implemented in 3D");
return mesh_data.ll();
}
return compute_mes_l;
}();
const CellValue<const double> dual_mes_l_j = [=] {
......@@ -167,16 +157,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
return computed_alpha_l;
}();
SparseMatrixDescriptor S{m_mesh->numberOfCells()};
const auto& primal_face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
NodeValue<bool> primal_node_is_on_boundary(m_mesh->connectivity());
if (parallel::size() > 1) {
throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
}
primal_node_is_on_boundary.fill(false);
for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
if (face_to_cell_matrix[face_id].size() == 1) {
for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
primal_node_is_on_boundary[node_id] = true;
}
}
}
SparseMatrixDescriptor S{m_mesh->numberOfCells()};
for (FaceId face_id = 0; face_id < m_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];
const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
// CellValue<double> dual_cell_id{diamond_mesh->connectivity()};
// mapper->toDualCell(face_id, dual_cell_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) {
......@@ -186,24 +190,104 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
} else {
S(cell_id1, cell_id2) += beta_l;
}
// S(cell_id1, cell_id2)+= alpha_l[face_id]*(,Cjr[dual_cell_id,dual_node_j]);
}
}
}
const CellValue<const double> primal_Vj = mesh_data.Vj();
CRSMatrix A{S};
Vector<double> b{m_mesh->numberOfCells()};
double pi = 3.14159265;
for (CellId i_cell = 0; i_cell < m_mesh->numberOfCells(); ++i_cell) {
if constexpr (Dimension == 2) {
b[i_cell] =
-2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
} else if constexpr (Dimension == 3) {
b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
FaceValue<const CellId> face_dual_cell_id = [=]() {
FaceValue<CellId> computed_face_dual_cell_id{m_mesh->connectivity()};
CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
parallel_for(
diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
return computed_face_dual_cell_id;
}();
NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
CellValue<NodeId> cell_ignored_id{m_mesh->connectivity()};
cell_ignored_id.fill(NodeId{std::numeric_limits<NodeId>::max()});
NodeValue<NodeId> node_primal_id{m_mesh->connectivity()};
parallel_for(
m_mesh->numberOfCells(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
return computed_dual_node_primal_node_id;
}();
const auto& dual_Cjr = diamond_mesh_data.Cjr();
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
const auto& primal_face_cell_is_reversed = m_mesh->connectivity().cellFaceIsReversed();
const auto& face_local_numbers_in_their_cells = m_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& primal_node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
if (face_to_cell_matrix[face_id].size() > 1) {
const double alpha_face_id = alpha_l[face_id];
for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
const bool is_face_reversed_for_cell_i =
primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
for (size_t j_face_cell = 0; j_face_cell < face_to_cell_matrix[face_id].size(); ++j_face_cell) {
CellId j_id = face_to_cell_matrix[face_id][j_face_cell];
CellId dual_cell_id = face_dual_cell_id[face_id];
for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
return -primal_nlr(face_id, i_node);
} else {
return primal_nlr(face_id, i_node);
}
}();
NodeId face_node_id = primal_face_to_node_matrix[face_id][i_node];
const auto& primal_node_to_cell = primal_node_to_cell_matrix[face_node_id];
if (not primal_node_is_on_boundary[face_node_id]) {
const double weight_rj = [&] {
for (size_t i_cell = 0; i_cell < primal_node_to_cell.size(); ++i_cell) {
if (j_id == primal_node_to_cell[i_cell]) {
return w_rj(face_node_id, i_cell);
}
}
Assert(false, "could not determine the weight");
return std::numeric_limits<double>::signaling_NaN();
}();
for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
if (dual_node_id == face_node_id) {
const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
S(i_id, j_id) += weight_rj * alpha_face_id * (nil, Clr);
}
}
} else {
std::cout << "*** ignore node " << face_node_id << '\n';
}
}
}
}
}
}
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};
Vector<double> b{m_mesh->numberOfCells()};
parallel_for(
m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
Vector<double> T{m_mesh->numberOfCells()};
T = 0;
......@@ -216,7 +300,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Vector<double> error{m_mesh->numberOfCells()};
parallel_for(
m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
m_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((error, error)) << "\n";
......@@ -228,7 +313,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
} else {
throw NotImplementedError("not done in 3d");
throw NotImplementedError("not done in 1d");
}
}
......@@ -236,16 +321,19 @@ template HeatDiamondScheme<1>::HeatDiamondScheme(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
template HeatDiamondScheme<2>::HeatDiamondScheme(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
template HeatDiamondScheme<3>::HeatDiamondScheme(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&);
......@@ -14,7 +14,8 @@ struct HeatDiamondScheme
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& kappa_id,
const FunctionSymbolId& f_id);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
......@@ -176,23 +176,24 @@ SchemeModule::SchemeModule()
this->_addBuiltinFunction("heat", std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&, const FunctionSymbolId&)>>(
const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
[](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) -> void {
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};
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};
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};
HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
default: {
......@@ -207,23 +208,24 @@ SchemeModule::SchemeModule()
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&, const FunctionSymbolId&)>>(
const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId)>>(
[](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) -> void {
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};
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};
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};
Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
break;
}
default: {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment