Skip to content
Snippets Groups Projects
Commit c34c9549 authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Add construction of node_kappar and node_betar

parent 24e24145
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/scheme/ScalarNodalScheme.cpp
+ 76
− 147
View file @ c34c9549
......@@ -128,8 +128,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k_bound,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
......@@ -246,7 +246,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_face_dof_number;
}();
const FaceValue<const bool> primal_face_is_neumann = [&] {
const FaceValue<const bool> face_is_neumann = [&] {
FaceValue<bool> face_is_neumann{mesh->connectivity()};
face_is_neumann.fill(false);
for (const auto& boundary_condition : boundary_condition_list) {
......@@ -269,64 +269,27 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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);
const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
const NodeValue<const TinyVector<Dimension>>& xr = mesh_data.xr();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
{
std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
{
std::shared_ptr mapper =
DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
CellValue<const TinyMatrix<Dimension>> nod_kappar = nod_k->cellValues();
CellValue<const TinyMatrix<Dimension>> nod_kapparb = nod_k_bound->cellValues();
CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
const CellValue<const double> Vj = mesh_data.Vj();
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
......@@ -338,92 +301,59 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}();
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()};
const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
NodeValue<const TinyMatrix<Dimension>> kappa;
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;
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
kappa[node_id] = zero;
const auto& node_to_cell = node_to_cell_matrix[node_id];
double weight = 0;
for (size_t j = 0; j < node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
kappa[node_id] += cell_kappaj[J] * local_weight;
weight += local_weight;
}
kappa[node_id] /= weight;
}
);
return kappa;
}();
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()};
const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
NodeValue<const TinyMatrix<Dimension>> kappa;
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;
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
kappa[node_id] = zero;
const auto& node_to_cell = node_to_cell_matrix[node_id];
double weight = 0;
for (size_t j = 0; j < node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
kappa[node_id] += cell_kappajb[J] * local_weight;
weight += local_weight;
}
kappa[node_id] /= weight;
});
return kappa;
}();
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()};
const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
NodeValue<const TinyMatrix<Dimension>> beta;
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;
}();
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);
}
}
}
}
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
beta[node_id] = zero;
for (size_t j = 0; j < node_to_cell.size(); ++j) {
const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cell[j];
beta[node_id] += tensorProduct(Cjr(J,R),xr[node_id]-xj[J]);
}
});
return computedClj;
return beta;
}();
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 = [&] {
FaceValue<double> alpha_j{mesh->connectivity()};
......@@ -447,23 +377,22 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return alpha_lb;
}();
const CellValue<const double> primal_Vj = mesh_data.Vj();
const Array<int> non_zeros{number_of_dof};
non_zeros.fill(Dimension);
non_zeros.fill(Dimension); //Modif pour optimiser
CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
const auto& face_to_cell = face_to_cell_matrix[face_id];
const double beta_l = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
const double betab_l = l2Norm(dualClj[face_id]) * alphab_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];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id1 = face_to_cell[i_cell];
for (size_t j_cell = 0; j_cell < face_to_cell.size(); ++j_cell) {
const CellId cell_id2 = 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]) {
if (face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
}
} else {
......@@ -475,12 +404,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const size_t j = cell_dof_number[cell_id];
S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
S(j, j) += (*alpha)[cell_id] * Vj[cell_id];
}
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (primal_face_is_dirichlet[face_id]) {
if (face_is_dirichlet[face_id]) {
S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
}
}
......@@ -491,7 +420,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
Vector<double> b{number_of_dof};
b = zero;
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
b[cell_dof_number[cell_id]] = fj[cell_id] * Vj[cell_id];
}
// Dirichlet on b^L_D
{
......@@ -552,8 +481,8 @@ ScalarNodalSchemeHandler::solution() const
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
const std::shared_ptr<const IDiscreteFunction>& nod_k_bound,
const std::shared_ptr<const IDiscreteFunction>& nod_k,
const std::shared_ptr<const IDiscreteFunction>& k_bound,
const std::shared_ptr<const IDiscreteFunction>& cell_k,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
......@@ -561,11 +490,11 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
if (not i_mesh) {
throw NormalError("primal discrete functions are not defined on the same mesh");
}
const std::shared_ptr i_dual_mesh = getCommonMesh({nod_k_bound, nod_k});
const std::shared_ptr i_dual_mesh = getCommonMesh({cell_k_bound, cell_k});
if (not i_dual_mesh) {
throw NormalError("dual discrete functions are not defined on the same mesh");
}
checkDiscretizationType({alpha, nod_k_bound, nod_k, f}, DiscreteFunctionType::P0);
checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
......@@ -581,8 +510,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;
......@@ -600,8 +529,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;
......@@ -619,8 +548,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment