Merge branch 'debug' into '105-calculate-material-volume'

# Conflicts:
#   cardmech/conf/elasticity.conf

cardmech/conf/elasticity.conf

@@ -45,7 +45,7 @@ PressureIterations = 1
 
 DebugPressure = 1
 
-MechLevel = 1
+MechLevel = 0
 MechPolynomialDegree = 1
 
 WithPrestress=false
@@ -148,8 +148,8 @@ Density = 0.00
 VolumetricSplit = 0
 Incompressible = false
 QuasiCompressiblePenalty = Ciarlet
-VolumetricPenalty = 2000000
-PermeabilityPenalty = 2000000
+VolumetricPenalty = 20000
+PermeabilityPenalty = 0
 
 #------------------------
 # Linear Material

cardmech/src/elasticity/assemble/LagrangeElasticity.cpp

@@ -71,15 +71,13 @@ double LagrangeElasticity::PressureEnergy(const cell &c, int face, int bc, const
   double W_Pressure = 0.0;
   VectorFieldFaceElement E(U, c, face);
 
-  Tensor Q = RotationMatrix(c.GetData());
-  Tensor QI = Invert(Q);
-
   for (int q = 0; q < E.nQ(); ++q) {
     double w = E.QWeight(q);
-    VectorField FN = E.QNormal(q);
+    //VectorField N = mat::cofactor(DeformationGradient(E, q, U)) * N;
+    VectorField N = E.QNormal(q);
 
     double localPressure = -pProblem.Pressure(currentTime, E.QPoint(q), bc);
-    W_Pressure += w * localPressure * (FN * E.VectorValue(q, U));
+    W_Pressure += w * localPressure * (N * E.VectorValue(q, U));
   }
 
   return W_Pressure;
@@ -97,7 +95,6 @@ void LagrangeElasticity::Residual(const cell &c, const Vector &U, Vector &r) con
   VectorFieldElement E(U, c);
   RowBndValues r_c(r, c);
 
-  // TODO: Should get the actual material used in this cell.
   Tensor Q = OrientationMatrix(c.GetData());
   const Material &cellMat = eProblem.GetMaterial(c);
 
@@ -166,23 +163,17 @@ void LagrangeElasticity::PressureBoundary(const cell &c, int face, int bc, const
   VectorFieldFaceElement E(U, c, face);
   RowValues r_c(r, E);
 
-  Tensor Q = RotationMatrix(c.GetData());
-  Tensor QT = transpose(Q);
-
   for (int q = 0; q < E.nQ(); ++q) {
     double w = E.QWeight(q);
+    //Tensor F = DeformationGradient(E, q, U);
+    //VectorField N = mat::cofactor(F) * E.QNormal(q);
     VectorField N = E.QNormal(q);
 
     double localPressure = -pProblem.Pressure(currentTime, E.QPoint(q), bc);
 
-    Tensor F_0 = One + Q * E.VectorGradient(q, *u_pre) * QT;
-
     for (int i = 0; i < E.size(); ++i) {
       for (int k = 0; k < c.dim(); ++k) {
         VectorField u_ik = E.VectorValue(q, i, k);
-        Tensor F_ik = Product(u_ik, N);
-
-
         r_c(i, k) -= w * (localPressure) * N * u_ik;
       }
     }