diff --git a/cardmech/conf/elasticity.conf b/cardmech/conf/elasticity.conf index c95371c5b4583398dfe0e5410dfa98765aea5041..f580d58a817a5a365970be74d2950eb5e444d6c3 100644 --- a/cardmech/conf/elasticity.conf +++ b/cardmech/conf/elasticity.conf @@ -46,7 +46,7 @@ PressureIterations = 1 DebugPressure = 1 MechLevel = 0 -MechPolynomialDegree = 2 +MechPolynomialDegree = 1 WithPrestress=false PrestressSteps=5 @@ -59,7 +59,7 @@ PrestressSteps=5 # Meshes # ###################### #Mesh = BenchmarkBeam2DTet -Mesh = FullEllipsoid +#Mesh = FullEllipsoid Mesh = QuarterEllipsoid Mesh = OrientedEllipsoid #Mesh = KovachevaHeart @@ -124,7 +124,7 @@ StretchFactorNormal = 4.0 # Material Parameters# ###################### -#ActiveMaterial = Linear +ActiveMaterial = Linear #ActiveMaterial = NeoHooke #ActiveMaterial = MooneyRivlin #ActiveMaterial = StVenantKirchhoff @@ -147,15 +147,15 @@ Density = 0.00 #------------------------ VolumetricSplit = 0 Incompressible = false -QuasiCompressiblePenalty = None +QuasiCompressiblePenalty = Ciarlet VolumetricPenalty = 20000 -PermeabilityPenalty = 20000 +PermeabilityPenalty = 0 #------------------------ # Linear Material #------------------------ -LinearMat_Lambda = 50 -LinearMat_Mu = 10 +LinearMat_Lambda = 50000 +LinearMat_Mu = 10000 #------------------------ # NeoHooke Material @@ -182,9 +182,9 @@ StVenantMat_Mu = 8000 # Benchmark Parameter GuccioneMat_C=10000 -GuccioneMat_bf=1 -GuccioneMat_bs=1 -GuccioneMat_bfs=1 +GuccioneMat_bf=10 +GuccioneMat_bs=10 +GuccioneMat_bfs=10 # Benchmark Parameter GuccioneMat_C=2000 diff --git a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp b/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp index 8b302a2ede4f7fa617625082e47d6a22d0e03544..cd5d4c8252276bf3613cea4434de2849b660fb52 100644 --- a/cardmech/src/elasticity/assemble/LagrangeElasticity.cpp +++ b/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; } }