From a78ec7ee82a90ee6fd3d5bf56cd1e151f88004b6 Mon Sep 17 00:00:00 2001
From: Alexandre GANGLOFF <alexandre.gangloff@cea.fr>
Date: Thu, 16 Jan 2025 14:37:31 +0100
Subject: [PATCH] Add bugged order 2 reconstruction
---
src/language/modules/SchemeModule.cpp | 56 +++-
src/scheme/Order2HyperelasticSolver.cpp | 296 ++++++++++++------
.../Order2LocalDtHyperelasticSolver.cpp | 190 ++++++-----
.../Order2LocalDtHyperelasticSolver.hpp | 11 +-
src/scheme/PolynomialReconstruction.cpp | 31 +-
5 files changed, 385 insertions(+), 199 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 33c1f7b50..751278276 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -999,7 +999,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver_order2",
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -1020,7 +1020,9 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2,
- const double& mu, const double& lambda, const double& dt1)
+ const double& lambda1, const double& mu1,
+ const double& lambda2, const double& mu2,
+ const double& gamma1, const double& gamma2 , const double& dt1)
-> std::tuple<
std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1035,9 +1037,55 @@ SchemeModule::SchemeModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+ .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
- bc_descriptor_list2, mu, lambda);
+ bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list2,
+ const double& lambda1, const double& mu1,
+ const double& lambda2, const double& mu2,
+ const double& gamma1, const double& p_inf1,
+ const double& gamma2, const double& p_inf2, const double& dt1)
+ -> std::tuple<
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+ {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+ getCommonMesh(
+ {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+ .solver()
+ .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, p_inf1, gamma2, p_inf2);
}
));
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index c249ab60f..83b802f51 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -185,7 +185,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
_computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
- NodeValuePerCell<const Rdxd> DPk_sigma) const
+ DiscreteFunctionDPk<Dimension, const Rdxd> DPk_sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -206,7 +206,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
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_cells[j];
- br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma(J, R) * Cjr(J, R);
+ br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma[j](xr[r]) * Cjr(J, R);
}
b[r] = br;
@@ -338,7 +338,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
const NodeValuePerCell<const Rdxd>& Ajr,
const NodeValue<const Rd>& ur,
DiscreteFunctionDPk<Dimension,const Rd> DPk_uh,
- NodeValuePerCell<const Rdxd> DPk_sigma) const
+ DiscreteFunctionDPk<Dimension,const Rdxd> DPk_sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -357,13 +357,97 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
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_cells[j];
- F(J,R) = -Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma(J,R) * Cjr(J,R);
+ F(J,R) = -Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J,R);
}
});
return F;
}
+ void
+ _tensor_limiter(const MeshType& mesh,
+ const DiscreteFunctionP0<const Rdxd>& f,
+ DiscreteFunctionDPk<Dimension, Rdxd>& DPk_fh) const
+ {
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+ StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+
+ const auto xj = mesh_data.xj();
+
+ parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+ const Rdxd fj = f[cell_id];
+
+ Rdxd f_min = fj;
+ Rdxd f_max = fj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ for(size_t row = 0; row < Dimension; ++row){
+ for(size_t col = 0; col < Dimension; ++col){
+ f_min(row,col) = std::min(f_min(row,col), f[cell_stencil[i_cell]](row,col));
+ f_max(row,col) = std::max(f_max(row,col), f[cell_stencil[i_cell]](row,col));
+ }
+ }
+ }
+
+ Rdxd f_bar_min = fj;
+ Rdxd f_bar_max = fj;
+
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ const CellId cell_k_id = cell_stencil[i_cell];
+ const Rdxd f_xk = DPk_fh[cell_id](xj[cell_k_id]);
+
+ for(size_t row = 0; row < Dimension; ++row){
+ for(size_t col = 0; col < Dimension; ++col){
+ f_bar_min(row,col) = std::min(f_bar_min(row,col), f_xk(row,col));
+ f_bar_max(row,col) = std::max(f_bar_max(row,col), f_xk(row,col));
+ }
+ }
+ }
+
+ const double eps = 1E-14;
+ Rdxd coef1;
+ for(size_t row = 0; row < Dimension; ++row){
+ for(size_t col = 0; col < Dimension; ++col){
+ coef1(row,col) = 1;
+ if(std::abs(f_bar_min(row,col) - fj(row,col)) > eps){
+ coef1(row,col) = (f_max(row,col) - fj(row,col)) / (f_bar_max(row,col) - fj(row,col));
+ }
+ }
+ }
+
+ Rdxd coef2;
+ for(size_t row = 0; row < Dimension; ++row){
+ for(size_t col = 0; col < Dimension; ++col){
+ coef2(row,col) = 1;
+ if (std::abs(f_bar_min(row,col) - fj(row,col)) > eps) {
+ coef2(row,col) = (f_min(row,col) - fj(row,col)) / ((f_bar_min(row,col) - fj(row,col)));
+ }
+ }
+ }
+
+ double min_coef1 = coef1(0,0);
+ double min_coef2 = coef2(0,0);
+ for(size_t row = 0; row < Dimension; ++row){
+ for(size_t col = 0; col < Dimension; ++col){
+ min_coef1 = std::min(min_coef1,coef1(row,col));
+ min_coef2 = std::min(min_coef2,coef2(row,col));
+ }
+ }
+
+ const double lambda = std::max(0., std::min(1.,std::min(min_coef1, min_coef2)));
+
+ auto coefficients = DPk_fh.coefficients(cell_id);
+
+ coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ });
+ }
+
void
_vector_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rd>& f,
@@ -530,7 +614,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
const TinyMatrix<Dimension> I = identity;
const DiscreteScalarFunction& eps = E - 0.5 * dot(u,u);
- const DiscreteScalarFunction& p_d = (1-chi_solid)*(gamma_f-1) * rho * eps - p_inf_f*gamma_f + chi_solid * (gamma_s - 1) * rho * eps - gamma_s * p_inf_s;
+ const DiscreteScalarFunction& p_d = (1-chi_solid)*((gamma_f-1) * rho * eps)+ chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
const DiscreteTensorFunction& sigma_d = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma_f * p_d + chi_solid * (gamma_s*(p_d+p_inf_s))) / rho);
const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
@@ -567,7 +651,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
if(law == 1){
return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,mu,gamma_f,p_inf_f);
} else {
- return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,gamma_f,p_inf_f,gamma_s,p_inf_s);
+ return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid_d,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
}
}
@@ -584,7 +668,8 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
- }
+ } // - p* identity; // - chi_solid*P_zero_air_repo*identity;
+
if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
@@ -608,100 +693,109 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
symmetry_boundary_descriptor_list);
- CellValue<double> limiter_j{mesh.connectivity()};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
- limiter_j[j] = 1;
- }
- );
-
- const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
- std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
- std::vector<size_t> row;
- std::vector<size_t> col;
- for(size_t i = 0; i < Dimension; ++i){
- for(size_t k = i; k < Dimension; ++k){
-
- CellValue<double> coef{mesh.connectivity()};
- for(CellId j = 0; j <mesh.numberOfCells(); ++j){
- coef[j] = sigma_j[j](i,k);
- }
-
- const auto& coef_scalar_function = DiscreteScalarFunction(mesh_v, coef);
- auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
- auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
-
- const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef);
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
- });
-
- sigma_coef.push_back(copy(DPk_coef));
- row.push_back(i);
- col.push_back(k);
- }
- }
-
- NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
- const auto& cell_nodes = cell_to_node_matrix[j];
- for(size_t i = 0; i < Dimension; ++i){
- for(size_t k = 0; k < Dimension; ++k){
- for(size_t R = 0; R < cell_nodes.size(); ++R){
- sigma_lim(j,R)(i,k) = 0;
- }
- }
- }
- }
- );
-
- const NodeValue<const Rd>& xr = copy(mesh.xr());
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
-
- for(size_t R = 0; R < cell_nodes.size(); ++R){
- const NodeId r = cell_nodes[R];
-
- for(size_t l = 0; l < sigma_coef.size(); ++l){
- const size_t& i = row[l];
- const size_t& k = col[l];
-
- auto coefficients = sigma_coef[l].coefficients(j);
-
- coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
- for (size_t indice = 1; indice < coefficients.size(); ++indice) {
- coefficients[indice] *= limiter_j[j];
- }
-
- sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
- if(i != k){
- sigma_lim(j,R)(i,k) *= 2;
- }
- }
-
- sigma_lim(j,R) += transpose(sigma_lim(j,R));
- sigma_lim(j,R) *= 0.5;
- }
-
- }
- );
-
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
- auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-
- DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
+ //CellValue<double> limiter_j{mesh.connectivity()};
+ //parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+ // limiter_j[j] = 1;
+ // }
+ //);
+ //
+ //const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
+ //std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
+ //std::vector<size_t> row;
+ //std::vector<size_t> col;
+ //for(size_t i = 0; i < Dimension; ++i){
+ // for(size_t k = i; k < Dimension; ++k){
+ //
+ // CellValue<double> coef{mesh.connectivity()};
+ // for(CellId j = 0; j <mesh.numberOfCells(); ++j){
+ // coef[j] = sigma_j[j](i,k);
+ // }
+ //
+ // const auto& coef_scalar_function = DiscreteScalarFunction(mesh_v, coef);
+ // auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
+ // auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
+ // auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
+ //
+ // const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef);
+ // parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
+ // });
+ //
+ // sigma_coef.push_back(copy(DPk_coef));
+ // row.push_back(i);
+ // col.push_back(k);
+ // }
+ //}
+ //
+ //NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
+ //const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ //
+ //parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+ // const auto& cell_nodes = cell_to_node_matrix[j];
+ // for(size_t i = 0; i < Dimension; ++i){
+ // for(size_t k = 0; k < Dimension; ++k){
+ // for(size_t R = 0; R < cell_nodes.size(); ++R){
+ // sigma_lim(j,R)(i,k) = 0;
+ // }
+ // }
+ // }
+ // }
+ //);
+ //
+ //const NodeValue<const Rd>& xr = copy(mesh.xr());
+ //parallel_for(
+ // mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // const auto& cell_nodes = cell_to_node_matrix[j];
+ //
+ // for(size_t R = 0; R < cell_nodes.size(); ++R){
+ // const NodeId r = cell_nodes[R];
+ //
+ // for(size_t l = 0; l < sigma_coef.size(); ++l){
+ // const size_t& i = row[l];
+ // const size_t& k = col[l];
+ //
+ // auto coefficients = sigma_coef[l].coefficients(j);
+ //
+ // coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
+ // for (size_t indice = 1; indice < coefficients.size(); ++indice) {
+ // coefficients[indice] *= limiter_j[j];
+ // }
+ //
+ // sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
+ // if(i != k){
+ // sigma_lim(j,R)(i,k) *= 2;
+ // }
+ // }
+ //
+ // sigma_lim(j,R) += transpose(sigma_lim(j,R));
+ // sigma_lim(j,R) *= 0.5;
+ // }
+ //
+ // }
+ //);
+
+ std::cout << "Reconstruction" << "\n";
+ auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v,sigma_v);
+ auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto DPk_sigmah = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+ std::cout << "Reconstruction Ok ! " << "\n";
+
+ std::cout << "Limitation" << "\n";
+ DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
+ DiscreteFunctionDPk<Dimension, Rdxd> sigma_lim = copy(DPk_sigmah);
this->_vector_limiter(mesh,u,u_lim);
+ this->_tensor_limiter(mesh,sigma,sigma_lim);
+ std::cout << "Limitation OK !" << "\n";
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ std::cout << "br" << "\n";
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+ std::cout << "br Ok !" << "\n";
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -710,7 +804,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+ std::cout << "Fjr" << "\n";
NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+ std::cout << "Fjr Ok ! " << "\n";
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -956,14 +1052,14 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
std::shared_ptr<const DiscreteFunctionVariant> CG_app = CG;
auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
- std::tie(m_app,rho_app,u_app,E_app,CG_app) = apply_fluxes(dt/2, rho, u, E, CG, ur, Fjr);
+ //auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, Fjr);
- std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
+ //std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
- auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
+ //auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
- auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
- return order2_apply_fluxes(dt, rho, u, E, CG, rho_app, CG_app, ur_app, Fjr_app);
+ //auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
+ return apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
}
Order2HyperelasticSolver() = default;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index d74835e0b..a80ae83fe 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -672,6 +672,77 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return sigma_lim;
}
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_NeoHook(const DiscreteScalarFunction rho,
+ const DiscreteVectorFunction u,
+ const DiscreteScalarFunction E,
+ const DiscreteTensorFunction CG,
+ const size_t& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
+ {
+ const TinyMatrix<Dimension> I = identity;
+
+ const DiscreteScalarFunction& eps = E - 0.5 * dot(u,u);
+ const DiscreteScalarFunction& p_d = (1-chi_solid)*(gamma-1) * rho * eps - p_inf*gamma;
+ const DiscreteTensorFunction& sigma_d = chi_solid * (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
+ const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
+
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
+ std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ }
+
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_NeoHook2(const DiscreteScalarFunction rho,
+ const DiscreteVectorFunction u,
+ const DiscreteScalarFunction E,
+ const DiscreteTensorFunction CG,
+ const size_t& chi_solid,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
+ {
+ const TinyMatrix<Dimension> I = identity;
+
+ const DiscreteScalarFunction& eps = E - 0.5 * dot(u,u);
+ const DiscreteScalarFunction& p_d = (gamma-1) * rho * eps - p_inf*gamma;
+ const DiscreteTensorFunction& sigma_d = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
+ const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma * p_d + chi_solid * (gamma*(p_d+p_inf))) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
+
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
+ std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ }
+
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_state_law(const size_t law,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
+ const size_t& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
+ {
+ const DiscreteScalarFunction& rho_d = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u_d = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E_d = E_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& CG_d = CG_v->get<DiscreteTensorFunction>();
+
+ if(law == 1){
+ return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid,lambda,mu,gamma,p_inf);
+ } else {
+ return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid,lambda,gamma,p_inf);
+ }
+ }
+
public:
std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes(const SolverType& solver_type,
@@ -1307,6 +1378,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
_compute_classic_midpoint_method(const SolverType& solver_type,
+ const size_t& law,
const double& dt1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
@@ -1319,9 +1391,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
NodeValue<Rd> CR_ur,
NodeValuePerCell<Rd> CR_Fjr,
std::vector<NodeId>& map2,
- size_t fluid,
+ size_t chi_solid,
const double& lambda,
- const double& mu) const
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
{
std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma});
@@ -1339,15 +1413,12 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr<const DiscreteFunctionVariant> new_aL = aL;
std::shared_ptr<const DiscreteFunctionVariant> new_aT = aT;
- DiscreteScalarFunction new_aL_d = aL->get<DiscreteScalarFunction>();
- DiscreteScalarFunction new_aT_d = aT->get<DiscreteScalarFunction>();
-
- const TinyMatrix<Dimension> I = identity;
- const double& gamma = 1.4;
double sum_dt = 0;
while(sum_dt < dt1){
+ DiscreteScalarFunction new_aL_d = new_aL->get<DiscreteScalarFunction>();
+ DiscreteScalarFunction new_aT_d = new_aT->get<DiscreteScalarFunction>();
const std::shared_ptr<const DiscreteFunctionVariant>& new_c = std::make_shared<const DiscreteFunctionVariant>(new_aL_d + new_aT_d);
double dt2 = 0.4 * hyperelastic_dt(new_c);
@@ -1358,39 +1429,12 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
auto [ur, Fjr] = compute_fluxes(solver_type,new_rho,new_aL,new_aT,new_u,new_sigma,bc_descriptor_list,CR_ur,CR_Fjr,map2);
std::tie(m_app, rho_app, u_app, E_app, CG_app) = apply_fluxes(dt2/2,new_rho,new_u,new_E,new_CG,ur,Fjr);
- const DiscreteScalarFunction& rho_d = rho_app->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u_d = u_app->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& E_d = E_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteTensorFunction& CG_d = CG_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p_d = fluid*(gamma-1) * rho_d * eps;
- const DiscreteTensorFunction& sigma_d = (mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I) - p_d * I;
- const DiscreteScalarFunction& aL_d_app = sqrt((lambda + 2 * mu) / rho_d + gamma * p_d / rho_d);
- const DiscreteScalarFunction& aT_d_app = sqrt(mu / rho_d);
-
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma_app =
- std::make_shared<const DiscreteFunctionVariant>(sigma_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& aL_app =
- std::make_shared<const DiscreteFunctionVariant>(aL_d_app);
- const std::shared_ptr<const DiscreteFunctionVariant>& aT_app =
- std::make_shared<const DiscreteFunctionVariant>(aT_d_app);
+ const auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid,lambda,mu,gamma,p_inf);
auto [ur_app, Fjr_app] = compute_fluxes(solver_type,rho_app,aL_app,aT_app,u_app,sigma_app,bc_descriptor_list,CR_ur,CR_Fjr,map2);
std::tie(new_m,new_rho,new_u,new_E,new_CG) = order2_apply_fluxes(dt2,new_rho,new_u,new_E,new_CG,rho_app,CG_app,ur_app,Fjr_app);
- const DiscreteScalarFunction& new_rho_d = new_rho->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& new_u_d = new_u->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& new_E_d = new_E->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& new_eps = new_E_d - 0.5 * dot(new_u_d,new_u_d);
- const DiscreteTensorFunction& new_CG_d = CG_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& new_p_d = fluid*(gamma-1) * new_rho_d * new_eps;
- const DiscreteTensorFunction& new_sigma_d = (mu / sqrt(det(new_CG_d)) * (new_CG_d - I) + lambda / sqrt(det(new_CG_d)) * log(sqrt(det(new_CG_d))) * I) - new_p_d * I;
- new_aL_d = sqrt((lambda + 2 * mu) / new_rho_d + gamma * new_p_d / new_rho_d);
- new_aT_d = sqrt(mu / new_rho_d);
-
- new_sigma = std::make_shared<const DiscreteFunctionVariant>(new_sigma_d);
- new_aL = std::make_shared<const DiscreteFunctionVariant>(new_aL_d);
- new_aT = std::make_shared<const DiscreteFunctionVariant>(new_aT_d);
+ std::tie(new_sigma,new_aL,new_aT) = _apply_state_law(law,new_rho,new_u,new_E,new_CG,chi_solid,lambda,mu,gamma,p_inf);
sum_dt += dt2;
@@ -1507,6 +1551,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -1524,8 +1569,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const double& mu,
- const double& lambda) const
+ const double& lambda_f,
+ const double& mu_f,
+ const double& lambda_s,
+ const double& mu_s,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const
{
std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
@@ -1564,35 +1615,25 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::tie(m2_app,rho2_app,u2_app,E2_app,CG2_app) = apply_fluxes(dt2, rho2, u2, E2, CG2, ur2, Fjr2); //Solid t^n + dt2
- const double gamma = 1.4; //à modif
- const TinyMatrix<Dimension> I = identity;
- double sum_dt = dt2;
+ double sum_dt = dt2;
- size_t fluid = 0;
- if((mu==0) and (lambda==0)){
- fluid = 1;
+ size_t chi_solid = 0;
+ if((mu_s!=0) and (lambda_s!=0)){
+ chi_solid = 1;
}
+ const size_t chi_fluid = 1-chi_solid;
while(sum_dt < dt1/2){ //Partie a revoir avec loi de comportement diff
- const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u2_d = u2_app->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& E2_d = E2_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d, u2_d);
- const DiscreteTensorFunction& CG2_d = CG2_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
- const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
- const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
- const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
-
- sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
- aL2_app = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
- aT2_app = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+ std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+
+ DiscreteScalarFunction aL2_d_app = aL2_app->get<DiscreteScalarFunction>();
+ DiscreteScalarFunction aT2_d_app = aT2_app->get<DiscreteScalarFunction>();
const std::shared_ptr<const DiscreteFunctionVariant>& c_app =
std::make_shared<const DiscreteFunctionVariant>(aL2_d_app + aT2_d_app);
dt2 = 0.4 * hyperelastic_dt(c_app);
if(sum_dt + dt2 > dt1/2){
- dt2 = dt1/2 - sum_dt;
+ dt2 = dt1/2 - sum_dt;
}
auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, rho2_app, aL2_app, aT2_app, u2_app, sigma2_app,
@@ -1604,40 +1645,15 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
sum_dt += dt2;
}
- const DiscreteScalarFunction& rho1_d = rho1_app->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u1_d = u1_app->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& E1_d = E1_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps1 = E1_d - 0.5 * dot(u1_d,u1_d);
- const DiscreteTensorFunction& CG1_d = CG1_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p1_d = (1-fluid)*(gamma-1) * rho1_d * eps1;
- const DiscreteTensorFunction& sigma1_d = fluid * (mu / sqrt(det(CG1_d)) * (CG1_d - I) + lambda / sqrt(det(CG1_d)) * log(sqrt(det(CG1_d))) * I) - p1_d * I;
- const DiscreteScalarFunction& aL1_d = sqrt(fluid * (lambda + 2 * mu) / rho1_d + gamma * p1_d / rho1_d);
- const DiscreteScalarFunction& aT1_d = sqrt(fluid * mu / rho1_d);
-
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_app = std::make_shared<const DiscreteFunctionVariant>(sigma1_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& aL1_app = std::make_shared<const DiscreteFunctionVariant>(aL1_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& aT1_app = std::make_shared<const DiscreteFunctionVariant>(aT1_d);
-
- const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u2_d = u2_app->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& E2_d = E2_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d, u2_d);
- const DiscreteTensorFunction& CG2_d = CG2_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
- const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
- const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
- const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
-
- sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
- aL2_app = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
- aT2_app = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+ const auto& [sigma1_app, aL1_app, aT1_app] = _apply_state_law(law,rho1_app,u1_app,E1_app,CG1_app,chi_fluid,lambda_f,mu_f,gamma_f,p_inf_f);
+ std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes2(solver_type,rho1_app, aL1_app, aT1_app, u1_app, sigma1_app, bc_descriptor_list1,
rho2_app, aL2_app, aT2_app, u2_app, sigma2_app, bc_descriptor_list2, map1, map2); //Fluxes t^n + dt*0.5
const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, Fjr1_app);
- auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
- CR_ur_app,CR_Fjr_app,map2,fluid,lambda,mu);
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
+ CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index 58fe0b852..da09ae613 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -93,6 +93,7 @@ class Order2LocalDtHyperelasticSolverHandler
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -110,8 +111,14 @@ class Order2LocalDtHyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const double& mu,
- const double& lambda) const = 0;
+ const double& lambda_f,
+ const double& mu_f,
+ const double& lambda_s,
+ const double& mu_s,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const = 0;
IOrder2LocalDtHyperelasticSolver() = default;
IOrder2LocalDtHyperelasticSolver(IOrder2LocalDtHyperelasticSolver&&) = default;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 2af32e157..9c8a7bf18 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -41,6 +41,16 @@ symmetrize_coordinates(const TinyVector<Dimension>& origin,
return u - 2 * dot(u - origin, normal) * normal;
}
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_tensor(const TinyVector<Dimension>& normal,
+ const TinyMatrix<Dimension>& A)
+{
+ const TinyMatrix<Dimension>& I = identity;
+ const TinyMatrix<Dimension>& S = (I - 2*tensorProduct(normal,normal));
+ return S*A*S;
+}
+
class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
{
public:
@@ -927,12 +937,21 @@ 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_tensor(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