Commit 7c9303a0 authored by uheqb's avatar uheqb
Browse files

moved some methods from FiniteElasticity to LagrangeElasticity

parent 81e45e4f
Pipeline #218936 passed with stages
in 35 minutes and 48 seconds
#include <VectorFieldFaceElement.hpp>
#include <LagrangeDiscretization.hpp>
//#include <LagrangeDiscretization.hpp>
#include "Plotting.hpp"
#include "FiniteElasticity.hpp"
#include "VectorFieldElement.hpp"
......@@ -32,41 +32,6 @@ FiniteElasticity::FiniteElasticity(const ElasticityProblem &eP, const PressurePr
}
void
FiniteElasticity::dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const {
ScalarElement elem(U, c);
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
VectorField D(eProblem.Deformation(timeStep, elem.NodalPoint(k)));
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = D[d];
u_c.D(k, d) = true;
}
}
}
void FiniteElasticity::planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face, int d) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
void
FiniteElasticity::fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
}
void FiniteElasticity::Initialize(Vector &u, double t, double dt) {
currentTime = t;
timeStep = dt;
......@@ -123,8 +88,6 @@ void FiniteElasticity::Initialize(const cell &c, Vector &U) const {
}
}
void FiniteElasticity::Finalize(Vector &u, double t, double dt) {
UpdateAcceleration(u, dt);
}
......@@ -151,15 +114,6 @@ void FiniteElasticity::Initialize(Vectors &vecs, double t, double dt) {
UpdateStretch(vecs[1]);
}
void FiniteElasticity::UpdateStretch(const Vector &gamma_f) {
for (row r = gamma->rows(); r != gamma->rows_end(); ++r) {
double g = gamma_f(r, 0);
(*gamma)(r, 0) = g;
(*gamma)(r, 2) = factor_gn * g;
// To ensure det(FA)=1 [See Bonet et al.(2019)]
(*gamma)(r, 1) = (1.0 / (1.0 + g)) * (1.0 / (1.0 + factor_gn * g)) - 1.0;
}
}
void FiniteElasticity::Finalize(Vectors &vecs, double t, double dt) {
Finalize(vecs[0], t, dt);
......
......@@ -8,11 +8,17 @@
class FiniteElasticity : public CardiacAssemble {
void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
virtual void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const{
Warning("dirichletBnd is not defined")
};
void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const;
virtual void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const{
Warning("dirichletBnd is not defined")
};
void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
virtual void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const{
Warning("dirichletBnd is not defined")
};
protected:
const ElasticityProblem &eProblem;
......@@ -124,7 +130,9 @@ public:
return 0.0;
}
void UpdateStretch(const Vector &gamma_f);
virtual void UpdateStretch(const Vector &gamma_f) const {
Warning("UpdateStretch is not defined")
};
virtual void SetDisplacement(Vector &u) const {
Warning("SetDisplacement is not defined")
......
......@@ -465,6 +465,39 @@ double LagrangeElasticity::EnergyError(const Vector &u, const Vector &reference)
return EnergyNorm(u - projection);
}
void LagrangeElasticity::planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face, int d) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
void
LagrangeElasticity::fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c,
int face) const {
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = 0.0;
u_c.D(k, d) = true;
}
}
}
void LagrangeElasticity::dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const {
ScalarElement elem(U, c);
for (int j = 0; j < U.NodalPointsOnFace(c, face); ++j) {
int k = U.NodalPointOnFace(c, face, j);
VectorField D(eProblem.Deformation(timeStep, elem.NodalPoint(k)));
for (int d = 0; d < c.dim(); ++d) {
u_c(k, d) = D[d];
u_c.D(k, d) = true;
}
}
}
void LagrangeElasticity::InitializePrestress(const Vector &U) {
Prestress = std::make_unique<Vector>(0.0, U);
Prestress->Clear();
......@@ -824,6 +857,16 @@ void LagrangeElasticity::SetDisplacement(Vector &u) const {
}
}
void LagrangeElasticity::UpdateStretch(const Vector &gamma_f) const {
for (row r = gamma->rows(); r != gamma->rows_end(); ++r) {
double g = gamma_f(r, 0);
(*gamma)(r, 0) = g;
(*gamma)(r, 2) = factor_gn * g;
// To ensure det(FA)=1 [See Bonet et al.(2019)]
(*gamma)(r, 1) = (1.0 / (1.0 + g)) * (1.0 / (1.0 + factor_gn * g)) - 1.0;
}
}
void LagrangeElasticity::InitializeTimeStepping(double t, double dt, const Vector &U) {
auto u = new Vector(U.GetDisc(), U.SpaceLevel(), U.TimeLevel());
u->Clear();
......
......@@ -7,6 +7,12 @@ class LagrangeElasticity : public FiniteElasticity {
public:
LagrangeElasticity(const ElasticityProblem &eP, const PressureProblem &pP, bool isStatic = false);
void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const override;
void planeFixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face, int d) const override;
void dirichletBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const override;
void InitializePrestress(const Vector &u) override;
void InitializeTimeStepping(double t, double dt, const Vector &U) override;
......@@ -96,6 +102,8 @@ public:
void SetDisplacement(Vector &u) const override;
void UpdateStretch(const Vector &gamma_f) const override;
VectorField CalculateStretch(double t, const Point &xi) const;
};
......
Subproject commit c57a5a8d9610674fc2920fa8396852a9206c8fc8
Subproject commit e2c67db1a9a3d16accdd964f96b351c3a5fb672f
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment