diff --git a/conf/elasticity/benchmarks/LandBenchmarkBiVentricle.conf b/conf/elasticity/benchmarks/LandBenchmarkBiVentricle.conf
index c6a978d6f5bb32362fa00a339602ce6f1e586f6d..5e4a8e97e670859dd71bdfb6b4abd6daf25996d1 100644
--- a/conf/elasticity/benchmarks/LandBenchmarkBiVentricle.conf
+++ b/conf/elasticity/benchmarks/LandBenchmarkBiVentricle.conf
@@ -6,13 +6,13 @@ GeoPath = ../geo/
# ========================================================
h0=1e-4
h_min=1e-6
-#logfile = Unloaded_l1
+#logfile = AS_T4_1_deg1_gamma07
-#activeDeformation = Quarteroni2
-#g = 0.8 #gamma for updateStretch
-#scale_gamma = 0.8 #0.9, 1.0
+activeDeformation = Quarteroni2
+g = 0.8 #gamma for updateStretch
+scale_gamma = 0.8 #0.9, 1.0
#StretchFactorNormal =1
-activeStress = 0 #35
+#activeStress = 0 #35
Model = ActiveStrainElasticity
MechRuntype = Default
MechDynamics = Static
@@ -22,7 +22,7 @@ CalculateReference = 0
ReferenceLevel = 3
InterpolateStartVector = true
-Mesh = TestBiventricleInstantActivation
+Mesh = TestBiventricle
MechProblem = LandBiventricle
MechMeshPart = vessel_base
@@ -30,7 +30,7 @@ ProblemDimension = 3
ProblemGeometry = Tet
# Amount of pressure steps per iteration
-PressureSteps = 10
+PressureSteps = 4
PressureIterations = 20
PressureDepth = 0
diff --git a/src/elasticity/assemble/LagrangeElasticity.cpp b/src/elasticity/assemble/LagrangeElasticity.cpp
index e6d5c58d85d23fdc8449d1a6159515ee3ee9ab97..81e96654aa04e925d0c1209508d3849fbfc458e8 100644
--- a/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/src/elasticity/assemble/LagrangeElasticity.cpp
@@ -219,7 +219,7 @@ void LagrangeElasticity::ViscoResidual(const cell &c, const Vector &u, const Vec
mat::inverse_active_deformation(Q, stretch) : One;
Tensor F = DeformationGradient(E, q, u) * inverseFA;
auto Fiso = eProblem.IsochoricPart(F);
- Tensor L = E.VectorGradient(q, v) * inverseFA;
+ Tensor L = E.VectorGradient(q, v);// * inverseFA;?
VectorField f = eProblem.Load(Time(), E.QPoint(q), *c);
@@ -313,6 +313,9 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
VectorField stretch = E.VectorValue(q, *gamma);
if (cell_values) stretch = VectorField((*gamma_c)(c(), 0),(*gamma_c)(c(), 1),(*gamma_c)(c(), 2));
+// For Debugging
+// if ((*gamma_c)(c(), 0) != 0)
+// mout << "J " << c() << " gamma_f " << stretch[0] << endl;
Tensor F = DeformationGradient(E, q, U);
Tensor inverseFA = Contracts(c) ? mat::inverse_active_deformation(Q,stretch)
: One;
@@ -381,10 +384,13 @@ void LagrangeElasticity::ViscoJacobi(const cell &c, const Vector &u, const Vecto
VectorField stretch = E.VectorValue(q, *gamma);
if (cell_values) stretch = VectorField((*gamma_c)(c(), 0),(*gamma_c)(c(), 1),(*gamma_c)(c(), 2));
+// For Debugging
+// if ((*gamma_c)(c(), 0) != 0)
+// mout << "VJ " << c() << " gamma_f " << stretch[0] << endl;
Tensor inverseFA = Contracts(c) ?
mat::inverse_active_deformation(Q,stretch) : One;
Tensor F = DeformationGradient(E, q, u) * inverseFA;
- Tensor L = E.VectorGradient(q, v) * inverseFA;
+ Tensor L = E.VectorGradient(q, v);// * inverseFA;?
auto Fiso = eProblem.IsochoricPart(F);
@@ -400,6 +406,15 @@ void LagrangeElasticity::ViscoJacobi(const cell &c, const Vector &u, const Vecto
// Viscoelasic response
A_c(i, j, k, l) += w * cellMat.ViscoSecondDerivative(F, L, F_ik, F_jl);
+// For Debugging
+// double h = 0.0000001;
+// if (stretch[0] < -0.01)
+// mout << "c " << c()
+// << "D^2 " << cellMat.ViscoSecondDerivative(F, L, F_ik, F_jl)
+// << "D^2_h " << (cellMat.ViscoDerivative(F + h * F_jl, L, F_ik)
+// - cellMat.ViscoDerivative(F, L, F_ik)) / h
+// << endl;
+
// Passive Material response
A_c(i, j, k, l) += w * cellMat.TotalSecondDerivative(Fiso, F, Q, F_ik, F_jl);
}
diff --git a/src/elasticity/materials/ViscoelasticDamping.hpp b/src/elasticity/materials/ViscoelasticDamping.hpp
index 6f4307f2c548380f9efbf43fb45a1a002655732b..03883b16a3fd73ab1ab190d621f2b6f197967a85 100644
--- a/src/elasticity/materials/ViscoelasticDamping.hpp
+++ b/src/elasticity/materials/ViscoelasticDamping.hpp
@@ -37,25 +37,25 @@ public:
double Functional(const Tensor &F, const Tensor &L) const override {
auto traceEDot = trace(sym(transpose(F) * L));
- return eta * traceEDot * traceEDot;
- };
+ return 0.5 * eta * Frobenius(F, L) * Frobenius(F, L);
+ }
double Derivative(const Tensor &F, const Tensor &L, const Tensor &H) const override {
- return 2.0 * eta * Frobenius(F, L) * Frobenius(H, L);
+ return eta * Frobenius(F, L) * Frobenius(H, L);
}
double SecondDerivative(const Tensor &F, const Tensor &L, const Tensor &H,
const Tensor &G) const override {
- return 2.0 * eta * Frobenius(H, L) * Frobenius(G, L);
+ return eta * Frobenius(H, L) * Frobenius(G, L);
}
double Derivative(const Tensor &F, const Tensor &L, const TensorRow &H) const override {
- return 2.0 * eta * Frobenius(F, L) * Frobenius(H, L);
+ return eta * Frobenius(F, L) * Frobenius(H, L);
}
double SecondDerivative(const Tensor &F, const Tensor &L, const TensorRow &H,
const TensorRow &G) const override {
- return 2.0 * eta * Frobenius(H, L) * Frobenius(G, L);
+ return eta * Frobenius(H, L) * Frobenius(G, L);
}
};
diff --git a/tools/evaluation/EvaluationGAMMTests.py b/tools/evaluation/EvaluationGAMMTests.py
index 395a580a6b1577b4b46ca73e8f9e5e0f2078427c..f14495b5be6ba1069a1b6128a168c4f44bfb3b0e 100644
--- a/tools/evaluation/EvaluationGAMMTests.py
+++ b/tools/evaluation/EvaluationGAMMTests.py
@@ -68,11 +68,11 @@ def plotVolumeCurves(volume, volumeLabel, timeOrStep, timeOrStepLabel,lvRvLabelc
plt.xlabel(str(timeOrStepLabel))
plt.ylabel(str(volumeLabel))
plt.xlim(0, 800)
- plt.ylim(110, 160)
+ plt.ylim(60, 180)
for i in range(len(volume)):
plt.plot(timeOrStep, volume[i], label = lvRvLabel[i] + str(coupledProblem[i]) + " cV" + str(cV) + "dCV" + str(dCS))
- plt.legend(bbox_to_anchor=(1.8,0), loc='lower right', frameon=True)
- #plt.show()
+ plt.legend( loc='lower right', frameon=True)#bbox_to_anchor=(1.8,0),
+ plt.show()
#plt.savefig(foldername + 'VolumePlot.png') #Todo anpassen an data/LogFileName
@@ -80,8 +80,8 @@ def plotVolumeCurves(volume, volumeLabel, timeOrStep, timeOrStepLabel,lvRvLabelc
CoupledProblem = ["BiventricleCoarse","DeformedBiventricleCoarse"]
CoupledProblemDeformed = ["DeformedBiventricleCoarse"]
-cell_values = [0, 1]
-deformConcentrationStretch = [0, 1]
+cell_values = [0] #, 1
+deformConcentrationStretch = [0]#, 1
def checkIfFolderExists(path):
folderPath =pathToMakeFolders+path
@@ -97,7 +97,8 @@ if __name__=="__main__":
for cP in CoupledProblem:
for cV in cell_values:
for dCS in deformConcentrationStretch:
- foldername = path + pathExperiment + 'log' + str(cP) + '_cellV'+ str(cV) +'_deformCS' + str(dCS) + '.log'
+ foldername = 'DeformedRayleigh_l1.log'
+ #foldername = path + pathExperiment + 'log' + str(cP) + '_cellV'+ str(cV) +'_deformCS' + str(dCS) + '.log'
#foldernameDeformed = path + pathExperiment + 'log' + str(CoupledProblemDeformed[0]) + '_cellV' + str(cV) + '_deformCS' + str(dCS) + '.log'
vol, dataTime, dataStep = readDataFromLog(pathExperiment, foldername)