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

Remove alpha and construction of the matrix for Neumann BC

parent 54913220
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/language/modules/SchemeModule.cpp
+ 15
− 16
View file @ 6a6be095
......@@ -603,21 +603,20 @@ SchemeModule::SchemeModule()
}
));
this->_addBuiltinFunction(
"nodalheat",
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
this->_addBuiltinFunction("nodalheat",
std::make_shared<BuiltinFunctionEmbedder<
std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(
const std::vector<std::shared_ptr<
const IBoundaryConditionDescriptor>>&)>>(
[](const std::shared_ptr<const IDiscreteFunction>& alpha,
const std::shared_ptr<const IDiscreteFunction>& mub_dual,
const std::shared_ptr<const IDiscreteFunction>& mu_dual, const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-> const std::shared_ptr<const IDiscreteFunction> {
return ScalarNodalSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
[](const std::shared_ptr<const IDiscreteFunction>& mub_dual,
const std::shared_ptr<const IDiscreteFunction>& mu_dual,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
return ScalarNodalSchemeHandler{mub_dual, mu_dual, f, bc_descriptor_list}.solution();
}
));
......
src/scheme/ScalarNodalScheme.cpp
+ 101
− 108
View file @ 6a6be095
......@@ -27,8 +27,15 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
const Array<const NodeId> m_node_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_node_list;
}
const Array<const FaceId>&
faceList() const
{
......@@ -41,8 +48,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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}
DirichletBoundaryCondition(const Array<const FaceId>& face_list,
const Array<const NodeId>& node_list,
const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
{
Assert(m_value_list.size() == m_face_list.size());
}
......@@ -134,7 +143,6 @@ 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>>>& cell_k_bound,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
......@@ -157,6 +165,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
if constexpr (Dimension > 1) {
MeshFaceBoundary<Dimension> mesh_face_boundary =
getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
MeshNodeBoundary<Dimension> mesh_node_boundary =
getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
......@@ -167,7 +177,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
mesh_face_boundary
.faceList());
boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
boundary_condition_list.push_back(
DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
} else {
throw NotImplementedError("Dirichlet BC in 1d");
......@@ -451,86 +462,101 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return kappa_invBeta;
}();
const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
NodeValue<TinyVector<Dimension>> compute_sum_Cjr;
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
compute_sum_Cjr[node_id] = zero;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node = node_local_number_in_its_cell[i_cell];
compute_sum_Cjr[node_id] += Cjr(cell_id, i_node);
}
});
return compute_sum_Cjr;
}();
const NodeValue<const TinyVector<Dimension>> v = [&] {
NodeValue<TinyVector<Dimension>> compute_v;
parallel_for(
mesh->numberOfNodes(),
PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kappar[node_id] * exterior_normal[node_id]; });
return compute_v;
}();
const NodeValuePerCell<const double> theta = [&] {
NodeValuePerCell<double> compute_theta;
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_nodes = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
const NodeId node_id = cell_nodes[i_node];
compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]);
}
});
return compute_theta;
}();
const NodeValue<const double> sum_theta = [&] {
NodeValue<double> compute_sum_theta;
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
compute_sum_theta[node_id] = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node = node_local_number_in_its_cell[i_cell];
compute_sum_theta[node_id] += theta(cell_id, i_node);
}
});
return compute_sum_theta;
}();
const Array<int> non_zeros{mesh->numberOfCells()};
non_zeros.fill(Dimension); // Modif pour optimiser
CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
if (not is_boundary_node[node_id]) {
for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
const CellId J = node_to_cell[cell_id_j];
const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
const size_t node_id_j = node_local_number_in_j_cell[node_id];
for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
const CellId cell_id_j = node_to_cell[i_cell_j];
const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
const CellId K = node_to_cell[cell_id_k];
const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
const size_t node_id_k = node_local_number_in_k_cell[node_id];
for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
S(cell_id_j, cell_id_k) +=
dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
}
}
} else if (not node_is_corner[node_id]) {
}
}
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>) {
const auto& node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
const CellId J = node_to_cell[cell_id_j];
const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
const size_t node_id_j = node_local_number_in_j_cell[node_id];
for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
const CellId K = node_to_cell[cell_id_k];
const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
const size_t node_id_k = node_local_number_in_k_cell[node_id];
const TinyMatrix<Dimension> I = identity;
const TinyVector<Dimension> n = exterior_normal[node_id];
const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
const TinyMatrix<Dimension> P = I - nxn;
S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), P * Cjr(J, node_id_j));
}
}
}
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& node_list = bc.faceList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
const CellId cell_id_j = node_to_cell[i_cell_j];
const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
const CellId J = node_to_cell[cell_id_j];
const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
const size_t node_id_j = node_local_number_in_j_cell[node_id];
for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
const CellId K = node_to_cell[cell_id_k];
const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
const size_t node_id_k = node_local_number_in_k_cell[node_id];
S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
S(cell_id_j, cell_id_k) +=
dot(kappar_invBetar[node_id] *
(Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
P * Cjr(cell_id_k, i_node_j));
}
}
}
}
},
boundary_condition);
}
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
S(cell_id, cell_id) += Vj[cell_id];
......@@ -545,38 +571,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
b[cell_id] = fj[cell_id] * Vj[cell_id];
}
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>) { // To do
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_id] += value_list[i_face];
}
} else if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
const auto& node_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
const TinyVector<Dimension> n = exterior_normal[node_id];
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
CellId J = node_to_cell[i_cell];
const auto& node_local_number_in_cell = node_local_numbers_in_their_cells.itemArray(J);
const size_t node_id_j = node_local_number_in_cell[node_id];
b[J] -= dot(Cjr(J, node_id_j), n) * value_list[i_node];
}
}
}
},
boundary_condition);
}
Vector<double> T{mesh->numberOfCells()};
T = zero;
......@@ -598,18 +592,17 @@ ScalarNodalSchemeHandler::solution() const
}
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
const std::shared_ptr<const IDiscreteFunction>& cell_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)
{
const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
if (not i_mesh) {
throw NormalError("primal discrete functions are not defined on the same mesh");
}
checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
checkDiscretizationType({cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
......@@ -624,7 +617,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
std::make_unique<ScalarNodalScheme<2>>(mesh,
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
......@@ -639,7 +632,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
std::make_unique<ScalarNodalScheme<3>>(mesh,
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
......
src/scheme/ScalarNodalScheme.hpp
+ 4
− 7
View file @ 6a6be095
......@@ -40,15 +40,12 @@ class ScalarNodalSchemeHandler
std::shared_ptr<const IDiscreteFunction> solution() const;
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,
ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_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);
~ScalarNodalSchemeHandler();
};
#endif // SCALAR_NODAL_SCHEME_HPP
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment