conf/coupled/gamm-paper/coarse-biventricle.conf

-#WithPrestress = false
-CoupledProblem = BiventricleCoarse
-
-#VolumetricPenalty = 10000
-#PrestressSteps = 20
-#CoupledProblem = BiventricleCoarse
-
 CoupledProblem = DeformedBiventricleCoarse
 ActiveMaterial = Bonet
 PassiveMaterial = Bonet
@@ -22,14 +15,14 @@ GeoPath = ../geo/
 h0=1e-4
 h_min=1e-6
 
-
+#logfile=mitInnerFaces
 activeStress = 0
-Mesh = TestBiventricle
+Mesh = TestBiventricleSimpleExcitation_smooth#TestBiventricle # TestBiventricleSimpleExcitation_smooth
 
 ProblemDimension = 3
 ProblemGeometry = Tet
 
-WithPrestress = true
+WithPrestress = false
 PrestressSteps = 10
 PressureIterations = 20
 PressureDepth = 0
@@ -46,28 +39,30 @@ activeDeformation = Quarteroni2
 
 withPressure = 1
 
-MechDeltaTime  = 0.001 #0.004# 0.001
-ElphyDeltaTime = 0.001 # 0.001 #0.000001
+MechDeltaTime  = 0.001 #0.0004              0.004# 0.001
+ElphyDeltaTime = 0.001 #0.0004              0.001 #0.000001
 
 #StretchFactorSheet = 2.0
 #StretchFactorNormal = 5.0
 
-linearImplicitSolverCell = 1
-withCellValues = 1
-deformConcentrationStretch = 1 #1
+########ANSCHAUEN############
+deformConcentrationStretch = 0 #1
+deformFibersInConductivity=true
 timeStepExtrapolation = 0
 DampExtrapolation = 0.5;
 
 Model = Coupled
 CoupledModel = Segregated
 
-#CoupledProblem = DeformedBiventricleCoarse
-CoupledProblem = BiventricleCoarse
+CoupledProblem = DeformedBiventricleCoarse
+ElphyProblem=DeformedBiventricleCoarse
+#CoupledProblem = BiventricleCoarse
+WithGerachDiffsion=true
 
 ElphyModelClass = MElphyModel
 ElphyModelName = TenTusscher
 TensionModelName = Rossi
-RossiViscosity = 500
+RossiViscosity = 1000#500
 #RossiViscosity = 7000
 RossiForce = 10
 #RossiForce = 7
@@ -82,12 +77,12 @@ ElphyPolynomialDegree = 1
 ######################
 StartTime = 0.0
 
-#EndTime=0.00001
+#EndTime=0.001
 EndTime = 0.8
 #EndTime = 0.8
 #EndTime = 0.02
 
-DeltaTime = 0.001
+DeltaTime = 0.001 #0.0004
 DeltaTimeMin = 0.00001
 DeltaTimeMax = 0.001
 
@@ -101,14 +96,16 @@ PlottingSteps = 10
 #######################
 ClearDistribution = 0
 
-MechDiscretization = Conforming
+MechDiscretization = EG#Conforming
+DGSign = -1
+DGPenalty = 200000#20000000
 MechRuntype = Default
-#MechDynamics = Newmark
+MechDynamics = Newmark #RayleighNewmark (s.u. mit RayleighAlpha = 200 und NewmarkBeta=0.25, NewmarkGamma=0.5)
 
 ## Rayleigh
 
-MechDynamics = RayleighNewmark
-RayleighAlpha = 100
+#MechDynamics = RayleighNewmark
+#RayleighAlpha = 200
 
 ##  KelvinVoigt
 
@@ -120,11 +117,11 @@ RayleighAlpha = 100
 # Newmark Parameters
 #------------------------
 # Damped Newmark
-#NewmarkBeta = 0.3
-#NewmarkGamma = 0.6
+NewmarkBeta = 0.3
+NewmarkGamma = 0.6
 # Undamped Newmark
-NewmarkBeta = 0.25
-NewmarkGamma = 0.5
+#NewmarkBeta = 0.25
+#NewmarkGamma = 0.5
 
 # If Iteration with Pressuresteps fails, recursively restart PressureSolver with Steps*2 until iteration limit is reached
 PressureIterations = 10
@@ -162,10 +159,10 @@ VolumetricPenalty = 1000 # ?
 PermeabilityPenalty = 0 #1000 # ?
 
 
-#ElphyModel = SemiImplicit
-ElphyModel = SemiImplicitOnCells
+ElphyModel = SemiImplicit
+#ElphyModel = SemiImplicitOnCells
 OrderStaggeredScheme = wcV
-ReassembleRHSonce = 0
+ReassembleRHSonce = 1 #hier
 
 CrankNicolsonTheta = 0.5
 ElphySplittingMethod =Godunov  #Strang, Stein,Godunov
@@ -173,7 +170,7 @@ ElphySplittingMethod =Godunov  #Strang, Stein,Godunov
 #Iext splitted in PDE =1, Iext in ODE =0
 IextInPDE=1
 #Glaettung von Anregung  in Zeit: Stepfunction,Arctan,Tanh
-ExternalCurrentSmoothingInTime = Stepfunction#Arctan #Stepfunction
+ExternalCurrentSmoothingInTime = Arctan#Arctan #Stepfunction
 #Glättung Anregung im Raum: discrete, linear
 ExternalCurrentSmoothingInSpace = linear
 
@@ -193,8 +190,8 @@ ElphyIntegrator=ExplicitEuler
 CellModelScheme = ExponentialIntegrator
 IonScheme = Quadratic # ExplicitEuler,RungeKutta2,RungeKutta4
 
-ThresholdVForActivation=-40.0
-ScaleExcitationAmplitude=0.04101
+ThresholdVForActivation=-50.0
+#ScaleExcitationAmplitude=0.04101
 
 EulerReAssemble = 0
 
@@ -202,8 +199,8 @@ SetSolution = 0
 
 # === Linear Solver === #
 LinearSolver = gmres
-LinearEpsilon = 1e-8
-LinearReduction = 1e-8
+LinearEpsilon = 1e-3
+LinearReduction = 1e-4
 LinearReduction = 1e-11
 LinearSteps = 10000
 #LinearPrintSteps = 1000
@@ -221,22 +218,22 @@ GatingReduction = 1e-12
 # === Elphy Solver ================================================
 ElphySolver = gmres
 ElphyPreconditioner=Jacobi
-ElphyEpsilon = 1e-6
-ElphyReduction = 1e-10
+ElphyEpsilon = 1e-3
+ElphyReduction = 1e-4
 
 # === Mechanics Solver ============================================
 MechSolver = gmres
 MechPreconditioner = GaussSeidel
-MechEpsilon = 1e-10
-MechReduction = 1e-10 #1e-12
+MechEpsilon = 1e-6
+MechReduction = 1e-6          #1e-12
 MechSteps = 20000
 
 # === Newton Method === #
-NewtonEpsilon = 1e-5 #1e-8
-NewtonReduction = 1e-12 #1e-5
+NewtonEpsilon = 1e-6#                  1e-8
+NewtonReduction = 1e-6
 NewtonSteps = 100
 NewtonLineSearchSteps = 10
-NewtonLinearizationReduction = 1e-4 #1e-12
+NewtonLinearizationReduction = 1e-6 #                 1e-12
 
 #Smoother = GaussSeidel
 #Smoother = PointBlockGaussSeidel
@@ -275,4 +272,4 @@ DynamicSolverVTK = -10
 ElphySolverVTK = -10
 CoupledSolverVTK = -10
 PlotCalcium=0
-PlotVTK = -1
\ No newline at end of file
+PlotVTK = -1

conf/m++.conf

 #loadconf = electrophysiology/m++.conf;
 #loadconf = elasticity/m++.conf;
-loadconf=coupled/electroMechConPDE.conf
+#loadconf=coupled/electroMechConPDE.conf
 #loadconf = coupled/electroMechCon.conf;
 
-#loadconf = coupled/gamm-paper/coarse-biventricle.conf
\ No newline at end of file
+loadconf = coupled/gamm-paper/coarse-biventricle.conf
\ No newline at end of file

mpp @ 5779d25f

-Subproject commit e448af94f8abb4d860b9ce064dd363c663353690
+Subproject commit 5779d25f3d7623188289550ae238443bf1deef54

src/coupled/solvers/SegregatedSolver.cpp

@@ -36,11 +36,11 @@ std::vector<std::string> metrics = {
     "FrobeniusNorm(gamma_v)",
     "L2Norm(gamma_v)",
     "H1Norm(gamma_v)",
-    "FrobeniusNorm(gamma_u)",
-    "L2Norm(gamma_u)",
-    "H1Norm(gamma_u)",
-    "NormFa(u)",
-    "NormPK(u)"
+    "FrobeniusNorm(gamma_u)"//,
+    //"L2Norm(gamma_u)",
+    //"H1Norm(gamma_u)"//,
+    //"NormFa(u)",
+    //"NormPK(u)"
 };
 
 void
@@ -69,10 +69,10 @@ SegregatedSolver::EvaluateIteration(const IElphyAssemble &elphyAssemble, const I
   evaluationMetrics[metrics[16]].emplace_back(elphyAssemble.L2Norm(v_stretch(elphyValues)));
   evaluationMetrics[metrics[17]].emplace_back(elphyAssemble.H1Norm(v_stretch(elphyValues)));
   evaluationMetrics[metrics[18]].emplace_back(u_stretch(mechValues).ConsistentNorm());
-  evaluationMetrics[metrics[19]].emplace_back(elphyAssemble.L2Norm(u_stretch(mechValues)));
-  evaluationMetrics[metrics[20]].emplace_back(elphyAssemble.H1Norm(u_stretch(mechValues)));
-  evaluationMetrics[metrics[21]].emplace_back(mechAssemble.NormFa(u_displacement(mechValues)));
-  evaluationMetrics[metrics[22]].emplace_back(mechAssemble.NormPK(u_displacement(mechValues)));
+  //evaluationMetrics[metrics[19]].emplace_back(elphyAssemble.L2Norm(u_stretch(mechValues)));
+  //evaluationMetrics[metrics[20]].emplace_back(elphyAssemble.H1Norm(u_stretch(mechValues)));
+  //evaluationMetrics[metrics[21]].emplace_back(mechAssemble.NormFa(u_displacement(mechValues)));
+  //evaluationMetrics[metrics[22]].emplace_back(mechAssemble.NormPK(u_displacement(mechValues)));
 
   std::unordered_map<std::string, double> solverInformation;
   mechSolver->FillInformation(mechAssemble, solverInformation);
@@ -93,7 +93,6 @@ SegregatedSolver::EvaluateIteration(const IElphyAssemble &elphyAssemble, const I
 
 void SegregatedSolver::printStep() {
   if (verbose < 1) return;
-  mout << "Segregated Solver Step Info:" << endl;
   for (const auto &metric: evaluationMetrics) {
     mout << "\t" << metric.first << " = " << metric.second.back() << endl;
   }

src/coupled/transfers/MultiPartTransfer.cpp

@@ -266,7 +266,134 @@ GetCardiacTransfer(const Vector &coarse, const Vector &fine) {
 
   if (coarseType == typeid(LagrangeDiscretization)) {
     return std::make_unique<MultiPartLagrangeTransfer>(coarse, fine);
+  } else if (coarseType == typeid(MixedEGDiscretization)) {
+    return std::make_unique<MultiPartEGTransfer>(coarse, fine);
   }
-
   THROW("Transfer type not implemented")
+}
+
+
+MultiPartEGTransfer::MultiPartEGTransfer(const Vector &coarse, const Vector &fine)
+    : MultiPartDeformationTransfer(), rowList(coarse, fine) {
+  auto &mechDisc = (const MixedEGDiscretization &) coarse.GetDisc();
+  auto &elphyDisc = (const MultiPartDiscretization &) fine.GetDisc();
+
+  mechDim = mechDisc.Size();
+  elphyDim = elphyDisc.Size();
+
+  discs = std::pair<std::shared_ptr<MixedEGDiscretization>, std::shared_ptr<MultiPartDiscretization>>{
+      std::make_shared<MixedEGDiscretization>(mechDisc, mechDisc.Degree(), mechDisc.Size()),
+      std::make_shared<MultiPartDiscretization>(elphyDisc, mechDisc.Size()),
+  };
+
+  constructRowList(coarse, fine);
+}
+
+void MultiPartEGTransfer::constructRowList(
+    const Vector &coarse, const Vector &fine) {
+  //dlevel berechnet die Anzahl der Verfeinerungsstufen auf Basis des Polynomgrades der coarse-Diskretisierung.
+  int dlevel = (int) log2(discs.first->Degree());
+  //lvl bestimmt die Verfeinerungsstufe zwischen coarse und fine durch die checkRefinement-Funktion.
+  double lvl = checkRefinement({discs.first->Degree(), discs.second->Degree()},
+                               {coarse.SpaceLevel(), fine.SpaceLevel()});
+
+  //Bestimmung gemeinsamer Subdomains
+  std::vector<int> coarseSD = coarse.GetMesh().IncludedSubdomains();
+  std::vector<int> fineSD = fine.GetMesh().IncludedSubdomains();
+  std::vector<int> intersectSD;
+  std::set_intersection(coarseSD.begin(), coarseSD.end(), fineSD.begin(), fineSD.end(),
+                        std::back_inserter(intersectSD));
+
+  std::vector<Point> childrenVerteces{};
+  std::vector<Point> childVerteces{};
+  std::vector<Point> coarseNodalPoints{};
+  for (cell c = coarse.cells(); c != coarse.cells_end(); ++c) {
+    //Überspringt Zellen, die nicht in den gemeinsamen Subdomains enthalten sind.
+    if (std::find(intersectSD.begin(), intersectSD.end(), c.Subdomain()) ==
+        intersectSD.end()) {
+      continue;
+    }
+
+    //Jedes parent wird rekursiv verfeinert, bis die gewünschte dlevel-Tiefe erreicht ist.
+    //Die neuen feineren Zellen (childVerteces) werden erstellt und zu recursiveCells hinzugefügt.
+    std::vector<std::vector<std::unique_ptr<Cell>>> recursiveCells(dlevel + 1);
+    recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
+    for (int l = 1; l <= dlevel; ++l) {
+      for (const auto &parent: recursiveCells[l - 1]) {
+        const std::vector<Rule> &R = parent->Refine(childrenVerteces);
+        for (int i = 0; i < R.size(); ++i) {
+          R[i](childrenVerteces, childVerteces);
+          recursiveCells[l].emplace_back(
+              CreateCell(R[i].type(), c.Subdomain(), childVerteces, false));
+        }
+      }
+
+    }
+    //Sobald die Zellen auf der feinsten Stufe (dlevel) erstellt sind, werden ihre Lagrange-Nodalpunkte berechnet.
+    //child ist eine verfeinerte (refinierte) Zelle enstanden aus einer coarse-Zelle
+    for (const auto &child: recursiveCells[dlevel]) {
+      LagrangeNodalPoints(*child, coarseNodalPoints, 1);
+      std::vector<int> coarseIds = mpp::nodalIds(coarse, coarseNodalPoints);
+      //linearTransfer überträgt die Werte durch Interpolation auf den fine-Vektor und speichert sie in rowList.
+      mpp::linearTransfer(child->ReferenceType(), coarseNodalPoints, coarseIds, lvl, fine, rowList);
+    }
+  }
+}
+
+void MultiPartEGTransfer::Project(Vector &coarse, const Vector &fine) const {
+  if (coarse.GetDisc().DiscName() != discs.first->DiscName() ||
+      fine.GetDisc().DiscName() != discs.second->DiscName()) {
+    THROW("Discretizations of transfer and vectors do not match")
+  }
+
+  coarse = 0.0;
+  for (row r = coarse.rows(); r != coarse.rows_end(); ++r) {
+    row fineR = fine.find_row(r());
+    if (fineR == fine.rows_end()) continue;
+    if (r.NumberOfDofs()==1) continue;
+
+    double scale = ((double) elphyDim) / ((double) fineR.NumberOfDofs());
+    for (int j = 0; j < fineR.NumberOfDofs(); ++j) {
+      coarse(r, (j * elphyDim) / fineR.NumberOfDofs()) += fine(fineR, j) * scale;
+    }
+  }
+  //coarse.ClearDirichletValues();
+}
+
+// Überträgt eine Lösung von einem groben Gitter auf ein feineres Gitter
+void MultiPartEGTransfer::Prolongate(const Vector &coarse, Vector &fine) const {
+  coarse.GetMesh().PrintInfo();
+
+  if (coarse.GetDisc().DiscName() != discs.first->DiscName() ||
+      fine.GetDisc().DiscName() != discs.second->DiscName()) {
+    THROW("Discretizations of transfer and vectors do not match")
+  }
+
+  fine = 0.0;
+  for (int i = 0; i < rowList.IdSize(); ++i) {
+    int n = fine.Dof(i);
+    int scale = n / elphyDim;
+    for (int j = 0; j < n; ++j) {
+      fine(i, j) = rowList.CalculateFineEntry(i, coarse, j / scale);
+    }
+  }
+  fine.ClearDirichletValues();
+}
+void MultiPartEGTransfer::ProlongateCell(const Vector &coarse, Vector &fine) const {
+  THROW("Not implemented yet!")
+}
+void MultiPartEGTransfer::ProlongateTransposed(Vector &coarse, const Vector &fine) const {
+  THROW("Not implemented yet!")
+}
+void MultiPartEGTransfer::Restrict(Vector &coarse, const Vector &fine) const {
+  THROW("Not implemented yet!")
+}
+void MultiPartEGTransfer::ProjectCells(Vector &coarse, const Vector &fine) const {
+  THROW("Not implemented yet!")
+}
+RMatrix MultiPartEGTransfer::AsMatrix() const {
+  THROW("Not implemented yet!")
+}
+std::unique_ptr<Vector> MultiPartEGTransfer::CreateFineDeformation() const {
+  THROW("Not implemented yet!")
 }
\ No newline at end of file

src/coupled/transfers/MultiPartTransfer.hpp

@@ -8,6 +8,7 @@
 
 #include "MultiPartDiscretization.hpp"
 #include "LagrangeDiscretization.hpp"
+#include "MixedEGDiscretization.hpp"
 
 double checkRefinement(std::pair<int, int> degrees, std::pair<int, int> levels);
 
@@ -71,6 +72,36 @@ public:
 
 };
 
+class MultiPartEGTransfer : public MultiPartDeformationTransfer {
+  int mechDim{1};
+  int elphyDim{1};
+
+  std::pair<std::shared_ptr<MixedEGDiscretization>, std::shared_ptr<MultiPartDiscretization>> discs{};
+  RowList rowList;
+
+
+  void constructRowList(const Vector &coarse, const Vector &fine);
+public:
+  MultiPartEGTransfer(const Vector &coarse, const Vector &fine);
+
+  void Prolongate(const Vector &coarse, Vector &fine) const override;
+  void ProlongateCell(const Vector &coarse, Vector &fine) const override;
+
+  void Project(Vector &coarse, const Vector &fine) const override;
+  void ProjectCells(Vector &coarse, const Vector &fine) const override;
+
+
+  [[nodiscard]] std::unique_ptr<Vector> CreateFineDeformation() const override;
+
+  void ProlongateTransposed(Vector &coarse, const Vector &fine) const override;
+
+  void Restrict(Vector &coarse, const Vector &fine) const override;
+
+  [[nodiscard]] RMatrix AsMatrix() const override;
+
+};
+
+
 std::unique_ptr<MultiPartDeformationTransfer>
 GetCardiacTransfer(const Vector &coarse, const Vector &fine);
 

src/elasticity/assemble/EGElasticity.hpp

@@ -8,12 +8,13 @@
 #include "IElasticity.hpp"
 #include "MixedEGDiscretization.hpp"
 #include "MixedEGVectorFieldElement.hpp"
-
+//TODO Reihenfolge so wie in CG machen
 class EGElasticity : public IElasticity {
   double sign{-1};
   double penalty{1};
   int size{3};
 
+  std::string name;
   std::shared_ptr<const MixedEGDiscretization> disc;
 
   void fixationBnd(const cell &c, const Vector &U, RowBndValues &u_c, int face) const;
@@ -34,32 +35,36 @@ public:
 
   std::shared_ptr<const IDiscretization> GetSharedDisc() const override {return disc; }
 
-  //void Initialize(Vector &u) const override {};
-
   void Initialize(const cell &c, Vector &U) const override;
 
-  void Initialize(Vectors &u) const override;
+  void Initialize(Vectors &vecs) const override;
 
   double Energy(const Vector &u) const override;
 
   void Energy(const cell &c, const Vector &u, double &energy) const override;
 
-
-
   void Residual(const cell &c, const Vector &u, Vector &r) const override;
 
+  void ViscoResidual(const cell &c, const Vector &u, const Vector &v, Vector &r) const override;
+
   void PressureBoundary(const cell &c, int face, int bc, const Vector &U, Vector &r) const;
 
   using IAssemble::Jacobi;
 
   void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
 
-  void TractionStiffness(Matrix &matrix) const override;
+  void ViscoJacobi(const cell &c, const Vector &u, const Vector &v, Matrix &A) const override;
 
   void MassMatrix(Matrix &massMatrix) const override;
 
   void SystemMatrix(Matrix &systemMatrix) const override;
 
+  void TractionStiffness(Matrix &matrix) const override;
+
+  void TractionViscosity(Matrix &matrix) const override;
+
+  void EvaluateTractionStiffness(Vector &u) const override;
+
   void GetInvariant(const Vector &displacement, Vector &iota4) const override;
 
   double TractionEnergy(const cell &c, int face, int bc, const Vector &U) const;
@@ -116,13 +121,16 @@ public:
 
   std::pair<double, double> detF(const Vector &u) const override;
 
+  void UpdateStretch(const Vector &gamma_f) const override;
+
+  void UpdateStretch(const Vector &gamma_f, const Vector &gamma_f_c) const override;
+
   void SetInitialValue(Vector &u) override;
 
   void Scale (double s) const {
     Material &cellMat = eProblem.GetMaterial();
-    //cellMat.ScalePenalty(s*s);
+    cellMat.ScalePenalty(s*s);
   }
-  
 };
 
 

src/elasticity/assemble/IElasticity.hpp

@@ -81,8 +81,8 @@ public:
   using IAssemble::Residual;
 
   double Residual(const Vector &u, Vector &r) const override {
-    r=0;
     r = *Prestress;
+    r = 0;
     for (cell ce = u.cells(); ce != u.cells_end(); ++ce) {
       Residual(ce, u, r);
     }

src/elasticity/assemble/LagrangeElasticity.cpp

@@ -219,6 +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;?
 
 
@@ -257,11 +258,19 @@ void LagrangeElasticity::PressureBoundary(const cell &c, int face, int bc, const
                                           Vector &r) const {
   VectorFieldFaceElement FE(U, *c, face);
   RowValues r_c(r, FE);
+
+  Tensor Q = OrientationMatrix(c.GetData());
   const Material &cellMat = eProblem.GetMaterial(c);
 
   for (int q = 0; q < FE.nQ(); ++q) {
     double w = FE.QWeight(q);
+    VectorField stretch =  FE.VectorValue(q, *gamma);
+    if (withCellValues) stretch = VectorField((*gamma_c)(c(), 0), (*gamma_c)(c(), 1), (*gamma_c)(c(), 2));
+
     Tensor F = DeformationGradient(FE, q, U);
+    Tensor inverseFA = Contracts(c) ? mat::inverse_active_deformation(Q,stretch)
+                                    : One;
+    F = F * inverseFA;
     VectorField N = mat::cofactor(F) * FE.QNormal(q);
     if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) {
        N = FE.QNormal(q);
@@ -351,14 +360,20 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
       for (int q = 0; q < FE.nQ(); ++q) {
         double w = FE.QWeight(q);
         VectorField N = FE.QNormal(q);
+        VectorField stretch =  FE.VectorValue(q, *gamma);
+        if (withCellValues) stretch = VectorField((*gamma_c)(c(), 0), (*gamma_c)(c(), 1), (*gamma_c)(c(), 2));
+
         Tensor F = DeformationGradient(FE, q, U);
+        Tensor inverseFA = Contracts(c) ? mat::inverse_active_deformation(Q,stretch)
+                                        : One;
+        F = F * inverseFA;
         double localPressure = -eProblem.Pressure(Time(), FE.QPoint(q), E.Bnd(face));
         for (int i = 0; i < FE.NodalPoints(); ++i) {
           for (int k = 0; k < c.dim(); ++k) {
             auto u_ik = FE.VectorComponentValue(q, i, k);
             for (int j = 0; j < FE.NodalPoints(); ++j) {
               for (int l = 0; l < c.dim(); ++l) {
-                auto F_jl = FE.VectorRowGradient(q, j, l);
+                auto F_jl = FE.VectorRowGradient(q, j, l);// * inverseFA;
                 A_c(i, j, k, l) -= w * (localPressure) * (Cross(F, F_jl) * N) * u_ik;
               }
             }
@@ -429,14 +444,22 @@ void LagrangeElasticity::ViscoJacobi(const cell &c, const Vector &u, const Vecto
       for (int q = 0; q < FE.nQ(); ++q) {
         double w = FE.QWeight(q);
         VectorField N = FE.QNormal(q);
-        Tensor F = DeformationGradient(FE, q, u);
+
+        VectorField stretch =  FE.VectorValue(q, *gamma);
+        if (withCellValues) stretch = VectorField((*gamma_c)(c(), 0), (*gamma_c)(c(), 1), (*gamma_c)(c(), 2));
+
+        Tensor inverseFA = Contracts(c) ?
+                           mat::inverse_active_deformation(Q, stretch) : One;
+        Tensor F = DeformationGradient(FE, q, u) * inverseFA;
+        auto Fiso = eProblem.IsochoricPart(F);
+
         double localPressure = -eProblem.Pressure(Time(), FE.QPoint(q), E.Bnd(face));
         for (int i = 0; i < FE.NodalPoints(); ++i) {
           for (int k = 0; k < c.dim(); ++k) {
             auto u_ik = FE.VectorComponentValue(q, i, k);
             for (int j = 0; j < FE.NodalPoints(); ++j) {
               for (int l = 0; l < c.dim(); ++l) {
-                auto F_jl = FE.VectorRowGradient(q, j, l);
+                auto F_jl = FE.VectorRowGradient(q, j, l);// * inverseFA;
                 A_c(i, j, k, l) -= w * (localPressure) * (Cross(F, F_jl) * N) * u_ik;
               }
             }

src/elasticity/materials/Material.hpp

@@ -290,6 +290,16 @@ public:
     return viscoMat->SecondDerivative(F, L, H, G);
   }
 
+  double ViscoSecondDerivative(const Tensor &F, const Tensor &L, const Tensor &H,
+                               const TensorRow &G) const {
+    return viscoMat->SecondDerivative(F, L, H, Tensor(G));
+  }
+
+  double ViscoSecondDerivative(const Tensor &F, const Tensor &L, const TensorRow &H,
+                               const Tensor &G) const {
+    return viscoMat->SecondDerivative(F, L, Tensor(H), G);
+  }
+
   const std::string &Name() const { return hyperelMat->Name(); }
 
   double GetParameter(int i) const {