diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp index 81fbcea6a803d13bcb85e9b016732a94de7b6895..992f848ca961785a659363bc7f96cd66ff51172a 100644 --- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp +++ b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp @@ -70,22 +70,20 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( case IBoundaryConditionDescriptor::Type::symmetry: { const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); - if (sym_bc_descriptor.name() == "symmetry") { - 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 == sym_bc_descriptor.boundaryDescriptor()) { - if constexpr (Dimension > 1) { - Array<const FaceId> face_list = ref_face_list.list(); - boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list}); - } else { - throw NotImplementedError("Symmetry conditions are not supported in 1d"); - break; - } + 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 == sym_bc_descriptor.boundaryDescriptor()) { + if constexpr (Dimension > 1) { + Array<const FaceId> face_list = ref_face_list.list(); + boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list}); + } else { + throw NotImplementedError("Symmetry conditions are not supported in 1d"); } } } + break; } case IBoundaryConditionDescriptor::Type::dirichlet: { const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = @@ -265,10 +263,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( 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])); + primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]) && + (!primal_face_is_symmetry[face_id])); } - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); - + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr(); const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl(); @@ -278,6 +276,14 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( CellValuePerNode<double> w_rj{mesh->connectivity()}; FaceValuePerNode<double> w_rl{mesh->connectivity()}; + const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr(); + auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> { + TinyVector<Dimension> proj; + const TinyVector<Dimension> nil = primal_nlr(face_id, 0); + proj = x - ((x - xl[face_id]), nil) * nil; + return proj; + }; + for (size_t i = 0; i < w_rl.numberOfValues(); ++i) { w_rl[i] = std::numeric_limits<double>::signaling_NaN(); } @@ -334,7 +340,15 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( 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]; + if (primal_face_is_symmetry[face_id]) { + for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) { + const CellId cell_id = face_to_cell_matrix[face_id][j]; + TinyVector<Dimension> xproj = project_to_face(xj[cell_id], face_id); + A(i, node_to_cell.size() + cpt_face) = xproj[i - 1]; + } + } else { + A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1]; + } cpt_face++; } } @@ -467,12 +481,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( }(); TinyMatrix<Dimension, double> eye = identity; - // for (size_t i = 0; i < Dimension; ++i) { - // eye(i, i) = 1; - // } - const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr(); - const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed(); + const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed(); const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells(); SparseMatrixDescriptor S{number_of_dof * Dimension}; for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { @@ -494,6 +504,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( const CellId j_id = primal_face_to_cell[j_cell]; TinyMatrix<Dimension, double> M = beta_mu_l * eye + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil); + TinyMatrix<Dimension, double> N = tensorProduct(nil, nil); + if (i_cell == j_cell) { for (size_t i = 0; i < Dimension; ++i) { for (size_t j = 0; j < Dimension; ++j) { @@ -501,6 +513,12 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( if (primal_face_is_neumann[face_id]) { S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j); } + if (primal_face_is_symmetry[face_id]) { + S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) += + -((i == j) ? 1 : 0) + N(i, j); + S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) += + (i == j) ? 1 : 0; + } } } } else { @@ -653,26 +671,13 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl()); FaceValue<TinyVector<Dimension>> Speed_face{mesh->connectivity()}; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - for (size_t i = 0; i < Dimension; ++i) { - U[(cell_dof_number[cell_id] * Dimension) + i] = Uj[cell_id][i]; - } - } - for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - for (size_t i = 0; i < Dimension; ++i) { - if (primal_face_is_on_boundary[face_id]) { - U[(face_dof_number[face_id] * Dimension) + i] = Ul[face_id][i]; - } - } - } Vector r = A * U - b; - - std::cout << "initial residu = " << std::sqrt((r, r)) << '\n'; + std::cout << "initial (real) residu = " << std::sqrt((r, r)) << '\n'; BiCGStab{b, A, U, 10000, 1e-15}; r = A * U - b; - std::cout << "real residu = " << std::sqrt((r, r)) << '\n'; + std::cout << "final (real) residu = " << std::sqrt((r, r)) << '\n'; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { for (size_t i = 0; i < Dimension; ++i) { diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp index 2d6c298eb6ea26d6ff626a1c9de3ee4162ae5721..24988916ec799ecaa3e227ebc39d54878bfef8d3 100644 --- a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp +++ b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp @@ -1,11 +1,14 @@ #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 <memory> #include <variant> #include <vector> diff --git a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp index f1b9fd81ea81b8e39efb2be3eaa33b258aacda2b..eb19f89d46be3f5c2d73a6a50097ef82fa815339 100644 --- a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp +++ b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp @@ -12,7 +12,7 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor std::ostream& _write(std::ostream& os) const final { - os << "symmetry(" << *m_boundary_descriptor << ")"; + os << "symmetry(" << m_name << ',' << *m_boundary_descriptor << ")"; return os; }