From 10288fa170dba0d00d05641b99b2beb193d41201 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Mon, 6 Jan 2025 10:55:37 +0100
Subject: [PATCH] Second order reconstruction
---
src/scheme/HyperplasticSolverO2.cpp | 7 ++++--
src/scheme/PolynomialReconstruction.cpp | 29 ++++++++++++++++++++-----
2 files changed, 28 insertions(+), 8 deletions(-)
diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
index c99d0100c..1d6174dd1 100644
--- a/src/scheme/HyperplasticSolverO2.cpp
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -150,7 +150,10 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
double chi = 0;
const double normS2 = frobeniusNorm(Sd) * frobeniusNorm(Sd);
double limit2 = 2. / 3 * yield * yield;
- chi = std::max(0., std::sqrt(normS2) - std::sqrt(limit2));
+ const double power = 0.5;
+ if ((normS2 - limit2) > 0) {
+ chi = std::pow((normS2 - limit2), power);
+ }
return chi;
}
@@ -189,7 +192,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
symmetry_boundary_descriptor_list);
DiscreteFunctionDPk<Dimension, Rd> DPk_fh_lim = copy(DPk_fh);
const auto xj = mesh_data.xj();
-
+ // return DPk_fh_lim;
const double eps2 = 1E-14;
double cmin = 1. - eps;
double cmax = 1. + eps;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2af32e157..4ae9e11f4 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -25,6 +25,15 @@
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& M)
+{
+ const TinyMatrix<Dimension> Id = identity;
+ const TinyMatrix<Dimension> S = Id - 2. * tensorProduct(normal, normal);
+ return S * M * S;
+}
+
template <size_t Dimension>
PUGS_INLINE auto
symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
@@ -927,12 +936,20 @@ PolynomialReconstruction::_build(
} else if constexpr (is_tiny_matrix_v<DataType>) {
if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
(DataType::NumberOfColumns == MeshType::Dimension)) {
- throw NotImplementedError("TinyMatrix symmetries for reconstruction");
- }
- const DataType& qi_qj = discrete_function[cell_i_id] - qj;
- for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
- for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+ const DataType& qi = discrete_function[cell_i_id];
+ const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ }
+ }
+ } else {
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ }
}
}
}
--
GitLab