Select Git revision
Synchronizer.hpp
HeatDiamondAlgorithm.cpp 32.37 KiB
#include <language/algorithms/HeatDiamondAlgorithm.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/LeastSquareSolver.hpp>
#include <algebra/LinearSolver.hpp>
#include <algebra/LocalRectangularMatrix.hpp>
#include <algebra/SparseMatrixDescriptor.hpp>
#include <algebra/TinyVector.hpp>
#include <algebra/Vector.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItem.hpp>
#include <output/VTKWriter.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/Timer.hpp>
template <size_t Dimension>
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& f_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<Dimension>;
constexpr ItemType FaceType = [] {
if constexpr (Dimension > 1) {
return ItemType::face;
} else {
return ItemType::node;
}
}();
using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition,
SymmetryBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList boundary_condition_list;
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
for (const auto& bc_descriptor : bc_descriptor_list) {
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
throw NotImplementedError("NIY");
break;
}
case IBoundaryConditionDescriptor::Type::dirichlet: {
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if constexpr (Dimension > 1) {
for (size_t i_ref_face_list = 0;
i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
const RefId& ref = ref_face_list.refId();
if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const FaceId> face_list = ref_face_list.list();
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
face_list);
boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
}
}
}
break;
}
case IBoundaryConditionDescriptor::Type::fourier: {
throw NotImplementedError("NIY");
break;
}
case IBoundaryConditionDescriptor::Type::neumann: {
const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
++i_ref_face_list) {
const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
const RefId& ref = ref_face_list.refId();
if (ref == neumann_bc_descriptor.boundaryDescriptor()) {
const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
if constexpr (Dimension > 1) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const FaceId> face_list = ref_face_list.list();
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
face_list);
boundary_condition_list.push_back(NeumannBoundaryCondition{face_list, value_list});
}
}
}
break;
}
default: {
is_valid_boundary_condition = false;
}
}
if (not is_valid_boundary_condition) {
std::ostringstream error_msg;
error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation";
throw NormalError(error_msg.str());
}
}
if constexpr (Dimension > 1) {
const CellValue<const size_t> cell_dof_number = [&] {
CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
return compute_cell_dof_number;
}();
size_t number_of_dof = mesh->numberOfCells();
const FaceValue<const size_t> face_dof_number = [&] {
FaceValue<size_t> compute_face_dof_number{mesh->connectivity()}; compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
for (const auto& boundary_condition : boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
(std::is_same_v<T, FourierBoundaryCondition>) or
(std::is_same_v<T, SymmetryBoundaryCondition>) or
(std::is_same_v<T, DirichletBoundaryCondition>)) {
const auto& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
std::ostringstream os;
os << "The face " << face_id << " is used at least twice for boundary conditions";
throw NormalError(os.str());
} else {
compute_face_dof_number[face_id] = number_of_dof++;
}
}
}
},
boundary_condition);
}
return compute_face_dof_number;
}();
const FaceValue<const bool> primal_face_is_neumann = [&] {
FaceValue<bool> face_is_neumann{mesh->connectivity()};
face_is_neumann.fill(false);
for (const auto& boundary_condition : boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
(std::is_same_v<T, FourierBoundaryCondition>) or
(std::is_same_v<T, SymmetryBoundaryCondition>)) {
const auto& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
face_is_neumann[face_id] = true;
}
}
},
boundary_condition);
}
return face_is_neumann;
}();
const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
NodeValue<bool> primal_node_is_on_boundary(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 < 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;
}
}
}
FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
if (parallel::size() > 1) {
throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
}
primal_face_is_on_boundary.fill(false);
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (face_to_cell_matrix[face_id].size() == 1) {
primal_face_is_on_boundary[face_id] = true;
}
}
FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
if (parallel::size() > 1) {
throw NotImplementedError("Calculation of face_is_neumann is incorrect");
}
primal_face_is_dirichlet.fill(false);
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
}
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
CellValue<double> Tj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
FaceValue<double> Tl =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
NodeValue<double> Tr(mesh->connectivity());
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
CellValuePerNode<double> w_rj{mesh->connectivity()};
FaceValuePerNode<double> w_rl{mesh->connectivity()};
for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
w_rl[i] = std::numeric_limits<double>::signaling_NaN();
}
for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
Vector<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];
if (not primal_node_is_on_boundary[i_node]) {
LocalRectangularMatrix<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];
}
}
Vector<double> x{node_to_cell.size()};
x = 0;
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];
}
} else {
int nb_face_used = 0;
for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
FaceId face_id = node_to_face_matrix[i_node][i_face];
if (primal_face_is_on_boundary[face_id]) {
nb_face_used++;
}
}
LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
for (size_t j = 0; j < node_to_cell.size() + nb_face_used; 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];
}
}
for (size_t i = 1; i < Dimension + 1; i++) {
int cpt_face = 0;
for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
FaceId face_id = node_to_face_matrix[i_node][i_face];
if (primal_face_is_on_boundary[face_id]) {
A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
cpt_face++;
}
}
}
Vector<double> x{node_to_cell.size() + nb_face_used};
x = 0;
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];
}
int cpt_face = node_to_cell.size();
for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
FaceId face_id = node_to_face_matrix[i_node][i_face];
if (primal_face_is_on_boundary[face_id]) {
w_rl(i_node, i_face) = x[cpt_face++];
Tr[i_node] += w_rl(i_node, i_face) * Tl[face_id];
}
}
}
}
{
VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
vtk_writer.write(mesh,
{NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
0, true); // forces last output
}
{
std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh);
MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
std::shared_ptr mapper = DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
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());
VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
vtk_writer.write(diamond_mesh,
{NamedItemValue{"kappaj", dual_kappaj}, NamedItemValue{"coords", diamond_mesh->xr()},
NamedItemValue{"Vj", diamond_mesh_data.Vj()}, NamedItemValue{"xj", diamond_mesh_data.xj()}},
0, true); // forces last output
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 CellId> face_dual_cell_id = [=]() {
FaceValue<CellId> computed_face_dual_cell_id{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{mesh->connectivity()};
cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
NodeValue<NodeId> node_primal_id{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), 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;
}();
CellValue<NodeId> primal_cell_dual_node_id = [=]() { CellValue<NodeId> cell_id{mesh->connectivity()};
NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
parallel_for(
diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
return cell_id;
}();
Timer my_timer;
const auto& dual_Cjr = diamond_mesh_data.Cjr();
FaceValue<TinyVector<Dimension>> dualClj = [&] {
FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
parallel_for(
mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
CellId cell_id = primal_face_to_cell[i];
const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size(); i_dual_cell++) {
const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
if (face_dual_cell_id[face_id] == dual_cell_id) {
for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
i_dual_node++) {
const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
if (final_dual_node_id == dual_node_id) {
computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
}
}
}
}
}
});
return computedClj;
}();
FaceValue<TinyVector<Dimension>> nlj = [&] {
FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
parallel_for(
mesh->numberOfFaces(),
PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
return computedNlj;
}();
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_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
});
FaceValue<double> computed_alpha_l{mesh->connectivity()};
mapper->fromDualCell(alpha_j, computed_alpha_l);
VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
vtk_writer.write(diamond_mesh,
{NamedItemValue{"alphaj", alpha_j}, NamedItemValue{"xj", diamond_mesh_data.xj()},
NamedItemValue{"Vl", diamond_mesh_data.Vj()}, NamedItemValue{"|l|", dual_mes_l_j}},
0,
true); // forces last output
return computed_alpha_l;
}();
SparseMatrixDescriptor S{number_of_dof};
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 = 1. / Dimension * alpha_l[face_id] * mes_l[face_id];
const double beta_l = l2Norm(dualClj[face_id]) * 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_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
}
} else {
S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
}
}
}
}
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const double alpha_face_id = mes_l[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 = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
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];
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
return -nlj[face_id];
} else {
return nlj[face_id];
}
}();
CellId dual_cell_id = face_dual_cell_id[face_id];
for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
if (dual_node_primal_node_id[dual_node_id] == node_id) {
const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
const double a_ir = alpha_face_id * (nil, Clr);
for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
}
}
if (primal_node_is_on_boundary[node_id]) {
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
FaceId l_id = node_to_face_matrix[node_id][l_face];
if (primal_face_is_on_boundary[l_id]) {
S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
}
} }
}
}
}
}
}
}
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (primal_face_is_dirichlet[face_id]) {
S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
}
}
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{number_of_dof};
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
}
// Dirichlet on b^L_D
{
for (const auto& boundary_condition : boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
b[face_dof_number[face_id]] += value_list[i_face]; // sign
}
}
},
boundary_condition);
}
}
// EL b^L
for (const auto& boundary_condition : boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
(std::is_same_v<T, FourierBoundaryCondition>)) {
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
Assert(face_to_cell_matrix[face_id].size() == 1);
b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face]; // sign
}
}
},
boundary_condition);
}
Vector<double> T{number_of_dof};
T = 1;
LinearSolver solver;
solver.solveLocalSystem(A, T, b);
CellValue<double> Temperature{mesh->connectivity()};
FaceValue<double> Temperature_face{mesh->connectivity()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; }); parallel_for(
mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
if (primal_face_is_neumann[face_id]) {
Temperature_face[face_id] = T[face_dof_number[face_id]];
} else {
Temperature_face[face_id] = Tl[face_id];
}
});
Vector<double> error{mesh->numberOfCells()};
CellValue<double> cell_error{mesh->connectivity()};
Vector<double> face_error{mesh->numberOfFaces()};
double error_max = 0.;
size_t cell_max = 0;
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
});
parallel_for(
mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
if (primal_face_is_on_boundary[face_id]) {
face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
} else {
face_error[face_id] = 0;
}
});
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
if (error_max < std::abs(cell_error[cell_id])) {
error_max = std::abs(cell_error[cell_id]);
cell_max = cell_id;
}
}
std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n";
std::cout << "||Error||_2 (face)= " << std::sqrt((face_error, face_error)) << "\n";
std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((face_error, face_error)) << "\n";
{
VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
vtk_writer.write(mesh,
{NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
NamedItemValue{"Texact", Tj}, NamedItemValue{"error", cell_error}},
0,
true); // forces last output
}
}
} else {
throw NotImplementedError("not done in 1d");
}
}
template <size_t Dimension>
class HeatDiamondScheme<Dimension>::DirichletBoundaryCondition
{
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{ return m_value_list;
}
DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}
{
Assert(m_value_list.size() == m_face_list.size());
}
~DirichletBoundaryCondition() = default;
};
template <size_t Dimension>
class HeatDiamondScheme<Dimension>::NeumannBoundaryCondition
{
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}
{
Assert(m_value_list.size() == m_face_list.size());
}
~NeumannBoundaryCondition() = default;
};
template <size_t Dimension>
class HeatDiamondScheme<Dimension>::FourierBoundaryCondition
{
private:
const Array<const double> m_coef_list;
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
const Array<const double>&
coefList() const
{
return m_coef_list;
}
public:
FourierBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& coef_list,
const Array<const double>& value_list)
: m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
{
Assert(m_coef_list.size() == m_face_list.size());
Assert(m_value_list.size() == m_face_list.size());
}
~FourierBoundaryCondition() = default;
};
template <size_t Dimension>
class HeatDiamondScheme<Dimension>::SymmetryBoundaryCondition
{
private:
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
public:
SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
~SymmetryBoundaryCondition() = default;
};
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&);