diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index 6036a30f69fce9a6849c46d0f3382fa8b5d8c54e..54b61b70a7d31444afe7347adffd57919a66bbe3 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -63,14 +63,17 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
{
private:
using Rdxd = TinyMatrix<Dimension>;
+ using R3x3 = TinyMatrix<3>;
using Rd = TinyVector<Dimension>;
+ using R3 = TinyVector<3>;
using MeshType = Mesh<Connectivity<Dimension>>;
using MeshDataType = MeshData<Dimension>;
- using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
- using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
- using DiscreteTensorFunction = DiscreteFunctionP0<Dimension, const Rdxd>;
+ using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
+ using DiscreteTensorFunction = DiscreteFunctionP0<Dimension, const Rdxd>;
+ using DiscreteTensorFunction3d = DiscreteFunctionP0<Dimension, const R3x3>;
class FixedBoundaryCondition;
class PressureBoundaryCondition;
@@ -86,6 +89,86 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
using BoundaryConditionList = std::vector<BoundaryCondition>;
+ R3x3
+ to3D(const Rdxd tensor) const
+ {
+ R3x3 tensor3d = zero;
+ for (size_t i = 0; i < Dimension; i++) {
+ for (size_t j = 0; j < Dimension; j++) {
+ tensor3d(i, j) = tensor(i, j);
+ }
+ }
+ return tensor3d;
+ }
+
+ R3
+ to3D(const Rd vector) const
+ {
+ R3 vector3d = zero;
+ for (size_t i = 0; i < Dimension; i++) {
+ vector3d[i] = vector[i];
+ }
+ return vector3d;
+ }
+
+ CellValue<const R3x3>
+ to3D(const MeshType& mesh, const CellValue<Rdxd>& tensor) const
+ {
+ CellValue<R3x3> tensor3d{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor3d[j] = to3D(tensor[j]); });
+ return tensor3d;
+ }
+
+ CellValue<const R3>
+ to3D(const MeshType& mesh, const CellValue<Rd> vector) const
+ {
+ CellValue<R3> vector3d{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { vector3d[j] = to3D(vector[j]); });
+ return vector3d;
+ }
+
+ Rdxd
+ toDimension(const R3x3 tensor3d) const
+ {
+ Rdxd tensor = zero;
+ for (size_t i = 0; i < Dimension; i++) {
+ for (size_t j = 0; j < Dimension; j++) {
+ tensor(i, j) = tensor3d(i, j);
+ }
+ }
+ return tensor;
+ }
+
+ Rd
+ toDimension(const R3 vector3d) const
+ {
+ Rd vector = zero;
+ for (size_t i = 0; i < Dimension; i++) {
+ vector[i] = vector3d[i];
+ }
+ return vector;
+ }
+
+ CellValue<const Rdxd>
+ toDimension(const MeshType& mesh, const CellValue<R3x3>& tensor3d) const
+ {
+ CellValue<Rdxd> tensor{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor[j] = toDimension(tensor3d[j]); });
+ return tensor;
+ }
+
+ CellValue<const Rd>
+ toDimension(const MeshType& mesh, const CellValue<R3>& vector3d) const
+ {
+ CellValue<Rd> vector{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { vector[j] = to3D(vector3d[j]); });
+ return vector;
+ }
+
NodeValuePerCell<const Rdxd>
_computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
{
@@ -412,15 +495,22 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
throw NormalError("hypoelastic solver expects P0 functions");
}
- const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
- const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
- const DiscreteTensorFunction& sigmad = sigmad_v->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
- const Rdxd I = identity;
- const DiscreteTensorFunction& sigma = -p * I + sigmad;
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction3d& sigmad = sigmad_v->get<DiscreteTensorFunction3d>();
+ const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
+ const R3x3 I = identity;
+
+ std::shared_ptr<const MeshType> p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ CellValue<R3x3> tensor_values3d{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { tensor_values3d[j] = -p[j] * I + sigmad[j]; });
+ CellValue<const Rdxd> tensor_values = toDimension(mesh, tensor_values3d);
+ DiscreteTensorFunction sigma(p_mesh, tensor_values);
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
@@ -450,7 +540,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
- const DiscreteTensorFunction& sigmad,
+ const DiscreteTensorFunction3d& sigmad,
const DiscreteScalarFunction& mu,
const DiscreteScalarFunction& yieldy,
const NodeValue<const Rd>& ur,
@@ -475,7 +565,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
- CellValue<Rdxd> new_sigmad = copy(sigmad.cellValues());
+ CellValue<R3x3> new_sigmad = copy(sigmad.cellValues());
auto chi = [&](const double a, const double b) -> double {
if ((a >= 0) and (b > 0))
return 1;
@@ -496,24 +586,33 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
+ R3x3 gradv3d = to3D(gradv);
const double eps = 1.e-16;
- const Rdxd I = identity;
+ const R3x3 I = identity;
const double norms = l2Norm(sigmad[j]);
const double yieldf = norms - std::sqrt(2. / 3.) * yieldy[j];
- const Rdxd N = 1. / (norms + eps) * sigmad[j];
- const Rdxd D = 0.5 * (gradv + transpose(gradv));
- const Rdxd Dd = D - 1. / Dimension * trace(D) * I;
+ const R3x3 N = 1. / (norms + eps) * sigmad[j];
+ const R3x3 D = 0.5 * (gradv3d + transpose(gradv3d));
+ const R3x3 Dd = D - 1. / 3. * trace(D) * I;
const double scalarND = scalarProduct(N, D);
- const Rdxd Dp = chi(yieldf, scalarND) * scalarND * N;
- const Rdxd gradva = 0.5 * (gradv - transpose(gradv));
- const Rdxd jaumann = sigmad[j] * gradva - gradva * sigmad[j];
+ const R3x3 Dp = chi(yieldf, scalarND) * scalarND * N;
+ const R3x3 gradva = 0.5 * (gradv3d - transpose(gradv3d));
+ const R3x3 jaumann = sigmad[j] * gradva - gradva * sigmad[j];
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
+ // const R3x3 sigmad_tr = sigmad[j] + dt_over_Vj * (2 * mu[j] * Dd - jaumann);
+ // const double normsd = l2Norm(sigmad_tr);
+ // if (chi(yieldf, scalarND)) {
+ // new_sigmad[j] = std::sqrt(2. / 3.) * yieldy[j] / (normsd + eps) * sigmad_tr;
+ // } else {
+ // new_sigmad[j] = sigmad_tr;
+ // }
new_sigmad[j] += dt_over_Vj * (2 * mu[j] * (Dd - Dp) - jaumann);
if (trace(new_sigmad[j]) >= 1.e-6) {
- std::cout << j << " trace " << trace(new_sigmad[j]) << "\n";
+ std::cout << j << " trace " << trace(new_sigmad[j]) << " x-strain " << new_sigmad[j](0, 0) << " y-strain "
+ << new_sigmad[j](1, 1) << " z-strain " << new_sigmad[j](2, 2) << "\n";
}
});
@@ -525,7 +624,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
- std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_sigmad))};
+ std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction3d(new_mesh, new_sigmad))};
}
std::tuple<std::shared_ptr<const IMesh>,
@@ -552,15 +651,15 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
throw NormalError("hypoelastic solver expects P0 functions");
}
- return this->apply_fluxes(dt, //
- dynamic_cast<const MeshType&>(*i_mesh), //
- rho->get<DiscreteScalarFunction>(), //
- u->get<DiscreteVectorFunction>(), //
- E->get<DiscreteScalarFunction>(), //
- sigmad->get<DiscreteTensorFunction>(), //
- mu->get<DiscreteScalarFunction>(), //
- yieldy->get<DiscreteScalarFunction>(), //
- ur->get<NodeValue<const Rd>>(), //
+ return this->apply_fluxes(dt, //
+ dynamic_cast<const MeshType&>(*i_mesh), //
+ rho->get<DiscreteScalarFunction>(), //
+ u->get<DiscreteVectorFunction>(), //
+ E->get<DiscreteScalarFunction>(), //
+ sigmad->get<DiscreteTensorFunction3d>(), //
+ mu->get<DiscreteScalarFunction>(), //
+ yieldy->get<DiscreteScalarFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}