diff --git a/Kinematics.py b/Kinematics.py
index d520ca427952f88a2ca6b5a064eb6f87bbb90da7..da44204184b65b74633b68be352dbe653281a164 100644
--- a/Kinematics.py
+++ b/Kinematics.py
@@ -14,7 +14,6 @@ import dolfin
import numpy
import dolfin_cm as dcm
-from .Problem import Problem
################################################################################
diff --git a/Kinematics_Spherical.py b/Kinematics_Spherical.py
new file mode 100644
index 0000000000000000000000000000000000000000..ff675357e3a055d875832b21814908b4b3673747
--- /dev/null
+++ b/Kinematics_Spherical.py
@@ -0,0 +1,69 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2018-2020 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+# from builtins import *
+
+import dolfin
+import numpy
+
+import dolfin_cm as dcm
+from .Kinematics import Kinematics
+
+################################################################################
+
+class SphericalKinematics(Kinematics):
+
+
+
+ def __init__(self,
+ R,
+ Rho,
+ Rho_old,
+ dim=3):
+
+ self.I = dolfin.Identity(dim)
+
+ # self.U = dolfin.as_vector([R*Rho, 0., 0.])
+
+ self.Ft = dolfin.as_matrix([
+ [1+Rho+R*Rho.dx(0), 0 , 0 ],
+ [ 0 , 1+Rho, 0 ],
+ [ 0 , 0 , 1+Rho]])
+ self.Ft_old = dolfin.as_matrix([
+ [1+Rho_old+R*Rho_old.dx(0), 0 , 0 ],
+ [ 0 , 1+Rho_old, 0 ],
+ [ 0 , 0 , 1+Rho_old]])
+
+ self.Jt = dolfin.det(self.Ft )
+ self.Jt_old = dolfin.det(self.Ft_old)
+
+ self.Ct = self.Ft.T * self.Ft
+ self.Ct_old = self.Ft_old.T * self.Ft_old
+
+ self.Et = (self.Ct - self.I)/2
+ self.Et_old = (self.Ct_old - self.I)/2
+
+ self.Fe = self.Ft
+ self.Fe_old = self.Ft_old
+
+ self.Je = dolfin.det(self.Fe )
+ self.Je_old = dolfin.det(self.Fe_old)
+
+ self.Ce = self.Fe.T * self.Fe
+ self.Ce_old = self.Fe_old.T * self.Fe_old
+ self.Ce_inv = dolfin.inv(self.Ce)
+ self.ICe = dolfin.tr(self.Ce)
+
+ self.Ee = (self.Ce - self.I)/2
+ self.Ee_old = (self.Ce_old - self.I)/2
+
+ self.Fe_mid = (self.Fe_old + self.Fe)/2
+ self.Ce_mid = (self.Ce_old + self.Ce)/2
+ self.Ee_mid = (self.Ee_old + self.Ee)/2
diff --git a/Material_Elastic_Bulk_CiarletGeymonat.py b/Material_Elastic_Bulk_CiarletGeymonat.py
index c084522bc7f6f7d0b24306e8e4b8888504ff7e65..b780d12cd266572298d881c2e12513c665ef8362 100644
--- a/Material_Elastic_Bulk_CiarletGeymonat.py
+++ b/Material_Elastic_Bulk_CiarletGeymonat.py
@@ -27,9 +27,9 @@ class CiarletGeymonatBulkElasticMaterial(BulkElasticMaterial):
if ("lambda" in parameters):
self.lmbda = dolfin.Constant(parameters["lambda"])
elif ("E" in parameters) and ("nu" in parameters):
- self.E = dolfin.Constant(parameters["E"])
- self.nu = dolfin.Constant(parameters["nu"])
- self.lmbda = self.E*self.nu/(1+self.nu)/(1-2*self.nu) # MG20180516: in 2d, plane strain
+ E = dolfin.Constant(parameters["E"])
+ nu = dolfin.Constant(parameters["nu"])
+ self.lmbda = E*nu/(1+nu)/(1-2*nu) # MG20180516: in 2d, plane strain
@@ -37,13 +37,16 @@ class CiarletGeymonatBulkElasticMaterial(BulkElasticMaterial):
U=None,
C=None):
- if (C is None):
+ assert (U is not None) or (C is not None), "Must provide U or C. Aborting."
+ if (U is not None):
dim = U.ufl_shape[0]
I = dolfin.Identity(dim)
F = I + dolfin.grad(U)
+ JF = dolfin.det(F)
C = F.T * F
+ elif (C is not None):
+ JF = dolfin.sqrt(dolfin.det(C)) # MG20200207: Watch out! This is well defined for inverted elements!
- JF = dolfin.sqrt(dolfin.det(C))
IC = dolfin.tr(C)
C_inv = dolfin.inv(C)
diff --git a/Material_Elastic_Bulk_Pneumo.py b/Material_Elastic_Bulk_Pneumo.py
index 80125871037961d4fd186b4e575298e9d6ef4348..b42b692992c9098f5a73c2d3703e8dc70e9c9040 100644
--- a/Material_Elastic_Bulk_Pneumo.py
+++ b/Material_Elastic_Bulk_Pneumo.py
@@ -35,9 +35,19 @@ class PneumoBulkElasticMaterial(BulkElasticMaterial):
def get_free_energy(self,
- C):
+ U=None,
+ C=None):
+
+ assert (U is not None) or (C is not None), "Must provide U or C. Aborting."
+ if (U is not None):
+ dim = U.ufl_shape[0]
+ I = dolfin.Identity(dim)
+ F = I + dolfin.grad(U)
+ JF = dolfin.det(F)
+ C = F.T * F
+ elif (C is not None):
+ JF = dolfin.sqrt(dolfin.det(C)) # MG20200207: Watch out! This is well defined for inverted elements!
- JF = dolfin.sqrt(dolfin.det(C))
IC = dolfin.tr(C)
C_inv = dolfin.inv(C)
diff --git a/Material_Elastic_CiarletGeymonatNeoHookean.py b/Material_Elastic_CiarletGeymonatNeoHookean.py
index 525e2a285cb4e51a2e837f858dba96a8d184fb78..e2122b027ee98f6367302e572a50818ee814c358 100644
--- a/Material_Elastic_CiarletGeymonatNeoHookean.py
+++ b/Material_Elastic_CiarletGeymonatNeoHookean.py
@@ -25,7 +25,7 @@ class CiarletGeymonatNeoHookeanElasticMaterial(ElasticMaterial):
parameters):
self.bulk = dcm.CiarletGeymonatBulkElasticMaterial(parameters)
- self.dev = dcm.NeoHookeanDevElasticMaterial(parameters)
+ self.dev = dcm.NeoHookeanDevElasticMaterial(parameters)
diff --git a/Material_Elastic_CiarletGeymonatNeoHookeanMooneyRivlin.py b/Material_Elastic_CiarletGeymonatNeoHookeanMooneyRivlin.py
new file mode 100644
index 0000000000000000000000000000000000000000..9d7ee388d424c51da0a4e793da65102e880b5c5c
--- /dev/null
+++ b/Material_Elastic_CiarletGeymonatNeoHookeanMooneyRivlin.py
@@ -0,0 +1,46 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2018-2020 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+# from builtins import *
+
+import dolfin
+
+import dolfin_cm as dcm
+from .Material_Elastic import ElasticMaterial
+
+################################################################################
+
+class CiarletGeymonatNeoHookeanMooneyRivlinElasticMaterial(ElasticMaterial):
+
+
+
+ def __init__(self,
+ parameters):
+
+ self.bulk = dcm.CiarletGeymonatBulkElasticMaterial(parameters)
+ self.dev = dcm.NeoHookeanMooneyRivlinDevElasticMaterial(parameters)
+
+
+
+ def get_free_energy(self,
+ *args,
+ **kwargs):
+
+ Psi_bulk, Sigma_bulk = self.bulk.get_free_energy(
+ *args,
+ **kwargs)
+ Psi_dev, Sigma_dev = self.dev.get_free_energy(
+ *args,
+ **kwargs)
+
+ Psi = Psi_bulk + Psi_dev
+ Sigma = Sigma_bulk + Sigma_dev
+
+ return Psi, Sigma
diff --git a/Material_Elastic_Dev_MooneyRivlin.py b/Material_Elastic_Dev_MooneyRivlin.py
index 3fa58f5b97fce258f2896d2f69d996f23fc7f8f9..44fb24b12a26967342658db725fbf0131354c4b4 100644
--- a/Material_Elastic_Dev_MooneyRivlin.py
+++ b/Material_Elastic_Dev_MooneyRivlin.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2018-2019 ###
+### Created by Martin Genet, 2018-2020 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Material_Elastic_Dev_NeoHookean.py b/Material_Elastic_Dev_NeoHookean.py
index 1bde7aa7521eee6870254f0055750ffb02f7e98b..34d4dd40f286e5a140a318c24b82ce3ae4601be4 100644
--- a/Material_Elastic_Dev_NeoHookean.py
+++ b/Material_Elastic_Dev_NeoHookean.py
@@ -24,12 +24,16 @@ class NeoHookeanDevElasticMaterial(DevElasticMaterial):
def __init__(self,
parameters):
- if ("mu" in parameters):
- self.mu = dolfin.Constant(parameters["mu"])
+ if ("C1" in parameters):
+ self.C1 = dolfin.Constant(parameters["C1"])
+ elif ("mu" in parameters):
+ mu = dolfin.Constant(parameters["mu"])
+ self.C1 = mu/2
elif ("E" in parameters) and ("nu" in parameters):
- self.E = dolfin.Constant(parameters["E"])
- self.nu = dolfin.Constant(parameters["nu"])
- self.mu = self.E/2/(1+self.nu)
+ E = dolfin.Constant(parameters["E"])
+ nu = dolfin.Constant(parameters["nu"])
+ mu = E/2/(1+nu)
+ self.C1 = mu/2
@@ -37,21 +41,27 @@ class NeoHookeanDevElasticMaterial(DevElasticMaterial):
U=None,
C=None):
- if (C is None):
+ assert (U is not None) or (C is not None), "Must provide U or C. Aborting."
+ if (U is not None):
dim = U.ufl_shape[0]
I = dolfin.Identity(dim)
F = I + dolfin.grad(U)
+ JF = dolfin.det(F)
C = F.T * F
- else:
+ elif (C is not None):
assert (C.ufl_shape[0] == C.ufl_shape[1])
dim = C.ufl_shape[0]
I = dolfin.Identity(dim)
+ JF = dolfin.sqrt(dolfin.det(C)) # MG20200207: Watch out! This is well defined for inverted elements!
- JF = dolfin.sqrt(dolfin.det(C))
IC = dolfin.tr(C)
C_inv = dolfin.inv(C)
- Psi = (self.mu/2) * (IC - dim - 2*dolfin.ln(JF)) # MG20180516: in 2d, plane strain
- Sigma = self.mu * (I - C_inv)
+ if (dim == 2):
+ Psi = self.C1 * (IC - 2 - 2*dolfin.ln(JF)) # MG20200206: plane strain
+ Sigma = 2*self.C1 * (I - C_inv) # MG20200206: plane strain
+ elif (dim == 3):
+ Psi = self.C1 * (IC - 3 - 2*dolfin.ln(JF))
+ Sigma = 2*self.C1 * (I - C_inv)
return Psi, Sigma
diff --git a/Material_Elastic_Dev_NeoHookeanMooneyRivlin.py b/Material_Elastic_Dev_NeoHookeanMooneyRivlin.py
new file mode 100644
index 0000000000000000000000000000000000000000..200368b6eeb2b8e15812d36da65bf1febfdd5978
--- /dev/null
+++ b/Material_Elastic_Dev_NeoHookeanMooneyRivlin.py
@@ -0,0 +1,57 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2018-2020 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+# from builtins import *
+
+import dolfin
+
+import dolfin_cm as dcm
+from .Material_Elastic_Dev import DevElasticMaterial
+
+################################################################################
+
+class NeoHookeanMooneyRivlinDevElasticMaterial(DevElasticMaterial):
+
+
+
+ def __init__(self,
+ parameters):
+
+ if ("mu" in parameters):
+ mu = parameters["mu"]
+ parameters["C1"] = mu/4
+ parameters["C2"] = mu/4
+ elif ("E" in parameters) and ("nu" in parameters):
+ E = parameters["E"]
+ nu = parameters["nu"]
+ mu = E/2/(1+nu)
+ parameters["C1"] = mu/4
+ parameters["C2"] = mu/4
+
+ self.nh = dcm.NeoHookeanDevElasticMaterial(parameters)
+ self.mr = dcm.MooneyRivlinDevElasticMaterial(parameters)
+
+
+
+ def get_free_energy(self,
+ *args,
+ **kwargs):
+
+ Psi_nh, Sigma_nh = self.nh.get_free_energy(
+ *args,
+ **kwargs)
+ Psi_mr, Sigma_mr = self.mr.get_free_energy(
+ *args,
+ **kwargs)
+
+ Psi = Psi_nh + Psi_mr
+ Sigma = Sigma_nh + Sigma_mr
+
+ return Psi, Sigma
diff --git a/Material_Elastic_Jerome.py b/Material_Elastic_Jerome.py
index c26afb0595da1b0e0bc8fc135a490b5d2b0840f3..437d99b55377ee8ae6f6647e946f63802cbc7bc8 100644
--- a/Material_Elastic_Jerome.py
+++ b/Material_Elastic_Jerome.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2018-2019 ###
+### Created by Martin Genet, 2018-2020 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Material_Inelastic_Damage.py b/Material_Inelastic_Damage.py
index 08a855abe45f898b7fbeb371293600074da7f2ae..61aa483ae03f8b8199f8fd57ebfc39bd419b3b4a 100644
--- a/Material_Inelastic_Damage.py
+++ b/Material_Inelastic_Damage.py
@@ -50,7 +50,7 @@ class DamageInelasticMaterial(InelasticMaterial):
def set_internal_variables_mixed(self):
- self.d = self.problem.subsols["d"].func
+ self.d = self.problem.subsols["d"].func # MG20200117: why not subfunc?
self.d_old = self.problem.subsols["d"].func_old
self.problem.psi = (1 - self.d) * self.problem.psi_eff
@@ -119,7 +119,7 @@ class DamageInelasticMaterial(InelasticMaterial):
jac_form = dolfin.inner(
dolfin.diff(
self.problem.sigma,
- self.d),
+ ),
dolfin.derivative(
self.problem.kinematics.epsilon,
self.problem.subsols["u"].subfunc,
diff --git a/NonlinearSolver.py b/NonlinearSolver.py
index fdd86dcda5874aa4461e04bc3831c655de7f7c03..e5c16106403cad87f8674aad1e5dbeb36816337b 100644
--- a/NonlinearSolver.py
+++ b/NonlinearSolver.py
@@ -219,7 +219,7 @@ class NonlinearSolver():
self.res_old_norm = self.res_norm
# linear system: Assembly
- if (len(self.problem.directional_penalties)): #MG20190513: Cannot use point integral within assemble_system
+ if (len(self.problem.directional_penalties)): # MG20190513: Cannot use point integral within assemble_system
self.printer.print_str("Assembly…",newline=False)
timer = time.time()
dolfin.assemble_system(
diff --git a/Problem.py b/Problem.py
index a74764d989b9f9cf7e3d51326ad6b948be670271..83e5dad8df7b4b5554e12f37819a2692db553607 100644
--- a/Problem.py
+++ b/Problem.py
@@ -76,10 +76,10 @@ class Problem():
degree=1)
self.eR_expr = dolfin.Expression(
- cppcode=("+x[0]/sqrt(pow(x[0],2)+pow(x[1],2))", "+x[1]/sqrt(pow(x[0],2)+pow(x[1],2))"),
+ ("+x[0]/sqrt(pow(x[0],2)+pow(x[1],2))", "+x[1]/sqrt(pow(x[0],2)+pow(x[1],2))"),
element=self.local_basis_fe)
self.eT_expr = dolfin.Expression(
- cppcode=("-x[1]/sqrt(pow(x[0],2)+pow(x[1],2))", "+x[0]/sqrt(pow(x[0],2)+pow(x[1],2))"),
+ ("-x[1]/sqrt(pow(x[0],2)+pow(x[1],2))", "+x[0]/sqrt(pow(x[0],2)+pow(x[1],2))"),
element=self.local_basis_fe)
self.Q_expr = dolfin.as_matrix([[self.eR_expr[0], self.eR_expr[1]],
diff --git a/Problem_Hyperelasticity.py b/Problem_Hyperelasticity.py
index 1cf9728cd7d861298f612207f99b897c6cab3d8a..adf9b020467c982c94bb185877f8e296c14dc1b3 100644
--- a/Problem_Hyperelasticity.py
+++ b/Problem_Hyperelasticity.py
@@ -71,7 +71,7 @@ class HyperelasticityProblem(Problem):
def set_solution_degree(self,
- U_degree=1): #MG20190513: Should have different name, right?
+ U_degree=1): # MG20190513: Should have different name, right?
self.set_subsols(
U_degree=U_degree)
@@ -180,7 +180,7 @@ class HyperelasticityProblem(Problem):
self.subsols["U"].subfunc,
self.mesh_normals)**2 * loading.measure
- # for loading in directional_penalties: #MG20190513: Cannot use point integral within assemble_system
+ # for loading in directional_penalties: # MG20190513: Cannot use point integral within assemble_system
# self.Pi += (loading.val/2) * dolfin.inner(
# self.subsols["U"].subfunc,
# loading.N)**2 * loading.measure
diff --git a/Problem_RelaxedGrowth.py b/Problem_RelaxedGrowth.py
index b7eee410a698bd0d9c2c69c6c996371210e97c04..9f918d02351a0bcbcafa708e335a8439952010d8 100644
--- a/Problem_RelaxedGrowth.py
+++ b/Problem_RelaxedGrowth.py
@@ -29,7 +29,7 @@ class RelaxedGrowthProblem(HyperelasticityProblem):
w_relaxation=False,
w_unloaded_configuration=False):
- Problem.__init__(self)
+ HyperelasticityProblem.__init__(self)
self.w_incompressibility = w_incompressibility
self.w_growth = w_growth
@@ -71,11 +71,11 @@ class RelaxedGrowthProblem(HyperelasticityProblem):
name="thetag",
family="DG",
degree=degree)
- #self.add_tensor_subsol(
- #name="Fg",
- #family="DG",
- #degree=degree,
- #init_val=numpy.eye(self.dim))
+ # self.add_tensor_subsol(
+ # name="Fg",
+ # family="DG",
+ # degree=degree,
+ # init_val=numpy.eye(self.dim))
@@ -327,13 +327,13 @@ class RelaxedGrowthProblem(HyperelasticityProblem):
self.subsols["Up"].subfunc,
self.mesh_normals)**2 * loading.measure
- # for loading in directional_penalties: #MG20190513: Cannot use point integral within assemble_system
+ # for loading in directional_penalties: # MG20190513: Cannot use point integral within assemble_system
# self.Pi += (loading.val/2) * dolfin.inner(
# self.subsols["U"].subfunc,
# loading.N)**2 * loading.measure
#
# if (self.w_unloaded_configuration):
- # for loading in directional_penalties: #MG20190513: Cannot use point integral within assemble_system
+ # for loading in directional_penalties: # MG20190513: Cannot use point integral within assemble_system
# self.Pi += (loading.val/2) * dolfin.inner(
# self.subsols["Up"].subfunc,
# loading.N)**2 * loading.measure
diff --git a/Problem_Spherical.py b/Problem_Spherical.py
new file mode 100644
index 0000000000000000000000000000000000000000..cbf4bccc1a25582a674cdb87fadd80575ab3b2bc
--- /dev/null
+++ b/Problem_Spherical.py
@@ -0,0 +1,329 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2018-2020 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import math
+import numpy
+
+import dolfin_cm as dcm
+from .Problem import Problem
+
+################################################################################
+
+class SphericalProblem(Problem):
+
+
+
+ def __init__(self,
+ w_incompressibility=False):
+
+ Problem.__init__(self)
+
+ self.dim = 3
+
+ self.w_incompressibility = w_incompressibility
+
+
+
+ def set_mesh(self,
+ Ri,
+ Re,
+ N):
+
+ self.Ri = Ri
+ self.Re = Re
+ self.N = N
+
+ self.mesh = dolfin.IntervalMesh(self.N, self.Ri, self.Re)
+
+ self.dR = dolfin.Measure(
+ "dx",
+ domain=self.mesh)
+
+ # self.R = dolfin.MeshCoordinates(self.mesh)
+ self.R_fe = dolfin.FiniteElement(
+ family="CG",
+ cell=self.mesh.ufl_cell(),
+ degree=1)
+ self.R = dolfin.Expression(
+ "x[0]",
+ element=self.R_fe)
+
+ self.Vi0 = 4/3 * math.pi * self.Ri**3
+ self.Ve0 = 4/3 * math.pi * self.Re**3
+ self.Vm0 = dolfin.assemble(4 * math.pi * self.R**2 * self.dR)
+
+ self.Ri_sd = dolfin.CompiledSubDomain("near(x[0], x0) && on_boundary", x0=self.Ri)
+ self.Re_sd = dolfin.CompiledSubDomain("near(x[0], x0) && on_boundary", x0=self.Re)
+
+ self.Ri_id = 1
+ self.Re_id = 2
+
+ self.boundaries_mf = dolfin.MeshFunction("size_t", self.mesh, self.mesh.topology().dim()-1)
+ self.boundaries_mf.set_all(0)
+ self.Ri_sd.mark(self.boundaries_mf, self.Ri_id)
+ self.Re_sd.mark(self.boundaries_mf, self.Re_id)
+
+ self.dS = dolfin.Measure(
+ "ds",
+ domain=self.mesh,
+ subdomain_data=self.boundaries_mf)
+
+
+
+ def add_rho_subsol(self,
+ degree):
+
+ self.add_scalar_subsol(
+ name="Rho",
+ family="CG",
+ degree=degree)
+
+
+
+ def get_rho_function_space(self):
+
+ if (len(self.subsols) == 1):
+ return self.sol_fs
+ else:
+ return self.get_subsol_function_space(name="Rho")
+
+
+
+ def add_pressure_subsol(self,
+ degree):
+
+ if (degree == 0):
+ self.add_scalar_subsol(
+ name="P",
+ family="DG",
+ degree=0)
+ else:
+ self.add_scalar_subsol(
+ name="P",
+ family="CG",
+ degree=degree)
+
+
+
+ def get_pressure_function_space(self):
+
+ assert (len(self.subsols) > 1)
+ return self.get_subsol_function_space(name="P")
+
+
+
+ def set_subsols(self,
+ Rho_degree=1):
+
+ self.add_rho_subsol(
+ degree=Rho_degree)
+
+ if (self.w_incompressibility):
+ self.add_pressure_subsol(
+ degree=Rho_degree-1)
+
+
+
+ def set_solution(self,
+ Rho_degree=1):
+
+ self.set_subsols(
+ Rho_degree=Rho_degree)
+ self.set_solution_finite_element()
+ self.set_solution_function_space()
+ self.set_solution_functions()
+
+ quadrature_degree = max(1, 2*(Rho_degree-1))
+ self.set_quadrature_degree(
+ quadrature_degree=quadrature_degree)
+
+ self.set_foi_finite_elements_DG(
+ degree=0)
+ self.set_foi_function_spaces()
+
+
+
+ def set_kinematics(self):
+
+ self.kinematics = dcm.SphericalKinematics(
+ R=self.R,
+ Rho=self.subsols["Rho"].subfunc,
+ Rho_old=self.subsols["Rho"].func_old)
+
+ # self.add_foi(expr=self.kinematics.Fe, fs=self.mfoi_fs, name="F")
+ # self.add_foi(expr=self.kinematics.Je, fs=self.sfoi_fs, name="J")
+ # self.add_foi(expr=self.kinematics.Ce, fs=self.mfoi_fs, name="C")
+ # self.add_foi(expr=self.kinematics.Ee, fs=self.mfoi_fs, name="E")
+
+
+
+ def set_materials(self,
+ elastic_behavior=None,
+ elastic_behavior_dev=None,
+ elastic_behavior_bulk=None,
+ subdomain_id=None):
+
+ if (self.w_incompressibility):
+ assert (elastic_behavior is None)
+ assert (elastic_behavior_dev is not None)
+ assert (elastic_behavior_bulk is None)
+ else:
+ assert ((elastic_behavior is not None)
+ or ((elastic_behavior_dev is not None)
+ and (elastic_behavior_bulk is not None)))
+
+ subdomain = dcm.SubDomain(
+ problem=self,
+ elastic_behavior=elastic_behavior,
+ elastic_behavior_dev=elastic_behavior_dev,
+ elastic_behavior_bulk=elastic_behavior_bulk,
+ id=subdomain_id)
+
+ self.subdomains += [subdomain]
+
+ # self.add_foi(expr=subdomain.Sigma, fs=self.mfoi_fs, name="Sigma")
+ # self.add_foi(expr=subdomain.PK1 , fs=self.mfoi_fs, name="PK1" )
+ # self.add_foi(expr=subdomain.sigma, fs=self.mfoi_fs, name="sigma")
+
+
+
+ def set_variational_formulation(self,
+ normal_penalties=[],
+ directional_penalties=[],
+ surface_tensions=[],
+ surface0_loadings=[],
+ pressure0_loadings=[],
+ volume0_loadings=[],
+ surface_loadings=[],
+ pressure_loadings=[],
+ volume_loadings=[],
+ dt=None):
+
+ self.Pi = sum([subdomain.Psi * 4 * math.pi * self.R**2 * self.dR(subdomain.id) for subdomain in self.subdomains])
+ # print (self.Pi)
+
+ self.res_form = dolfin.derivative(
+ self.Pi,
+ self.sol_func,
+ self.dsol_test)
+
+ for loading in pressure_loadings:
+ self.res_form -= loading.val * 4 * math.pi * self.R**2 * dolfin.inv(self.kinematics.Ft)[0,0] * self.R * self.subsols["Rho"].dsubtest * loading.measure
+
+ self.jac_form = dolfin.derivative(
+ self.res_form,
+ self.sol_func,
+ self.dsol_tria)
+
+
+
+ def add_strain_qois(self):
+
+ self.add_qoi(
+ name="E_RR",
+ expr=self.kinematics.Ee[0,0] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name="E_TT",
+ expr=self.kinematics.Ee[1,1] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name="E_PP",
+ expr=self.kinematics.Ee[2,2] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name="E_RT",
+ expr=self.kinematics.Ee[0,1] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name="E_TP",
+ expr=self.kinematics.Ee[1,2] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name="E_PR",
+ expr=self.kinematics.Ee[2,0] * 4 * math.pi * self.R**2 * self.dR)
+
+
+
+ def add_J_qois(self):
+
+ self.add_qoi(
+ name="J",
+ expr=self.kinematics.Je * 4 * math.pi * self.R**2 * self.dR)
+
+
+
+ def add_stress_qois(self,
+ stress_type="cauchy"):
+
+ nb_subdomain = 0
+ for subdomain in self.subdomains:
+ nb_subdomain += 1
+
+ if nb_subdomain == 0:
+ if (stress_type in ("cauchy", "sigma")):
+ basename = "s_"
+ stress = self.sigma
+ elif (stress_type in ("piola", "PK2", "Sigma")):
+ basename = "S_"
+ stress = self.Sigma
+ elif (stress_type in ("PK1", "P")):
+ basename = "P_"
+ stress = self.PK1
+
+ elif nb_subdomain == 1:
+ if (stress_type in ("cauchy", "sigma")):
+ basename = "s_"
+ stress = self.subdomains[0].sigma
+ elif (stress_type in ("piola", "PK2", "Sigma")):
+ basename = "S_"
+ stress = self.subdomains[0].Sigma
+ elif (stress_type in ("PK1", "P")):
+ basename = "P_"
+ stress = self.subdomains[0].PK1
+
+ self.add_qoi(
+ name=basename+"XX",
+ expr=stress[0,0] * 4 * math.pi * self.R**2 * self.dR)
+ if (self.dim >= 2):
+ self.add_qoi(
+ name=basename+"YY",
+ expr=stress[1,1] * 4 * math.pi * self.R**2 * self.dR)
+ if (self.dim >= 3):
+ self.add_qoi(
+ name=basename+"ZZ",
+ expr=stress[2,2] * 4 * math.pi * self.R**2 * self.dR)
+ if (self.dim >= 2):
+ self.add_qoi(
+ name=basename+"XY",
+ expr=stress[0,1] * 4 * math.pi * self.R**2 * self.dR)
+ if (self.dim >= 3):
+ self.add_qoi(
+ name=basename+"YZ",
+ expr=stress[1,2] * 4 * math.pi * self.R**2 * self.dR)
+ self.add_qoi(
+ name=basename+"ZX",
+ expr=stress[2,0] * 4 * math.pi * self.R**2 * self.dR)
+
+
+
+ def add_P_qois(self):
+
+ nb_subdomain = 0
+ for subdomain in self.subdomains:
+ nb_subdomain += 1
+ # print nb_subdomain
+
+ if nb_subdomain == 0:
+ basename = "P_"
+ P = -1./3. * dolfin.tr(self.sigma)
+ elif nb_subdomain == 1:
+ basename = "P_"
+ P = -1./3. * dolfin.tr(self.subdomains[0].sigma)
+
+ self.add_qoi(
+ name=basename,
+ expr=P * 4 * math.pi * self.R**2 * self.dR)
diff --git a/Step.py b/Step.py
index cf2cd424b504faff5721e4a1b9f3a6fdb18366c6..061ceacee99687e5ee3bc0d614823f3e46dcc2f1 100644
--- a/Step.py
+++ b/Step.py
@@ -26,7 +26,7 @@ class Step():
dt_ini=None,
dt_min=None,
dt_max=None,
- constraints=None, #MG20180508: Do not use list as default value because it is static
+ constraints=None, # MG20180508: Do not use list as default value because it is static
normal_penalties=None,
directional_penalties=None,
surface_tensions=None,
diff --git a/TimeIntegrator.py b/TimeIntegrator.py
index 213fc3ff46b4e80bb9b191ee59a9fced88693445..a81191c95683db5d2cd5d3baa12fd8c273a5680e 100644
--- a/TimeIntegrator.py
+++ b/TimeIntegrator.py
@@ -28,7 +28,9 @@ class TimeIntegrator():
print_out=True,
print_sta=True,
write_qois=True,
- write_sol=True):
+ write_sol=True,
+ write_vtus=False,
+ write_xmls=False):
self.problem = problem
@@ -88,11 +90,22 @@ class TimeIntegrator():
self.functions_to_write += self.problem.get_fois_func_lst()
self.xdmf_file_sol = dcm.XDMFFile(
- filename=self.write_sol_filebasename+".xdmf",
- functions=self.functions_to_write)
+ filename=self.write_sol_filebasename+".xdmf",
+ functions=self.functions_to_write)
self.problem.update_fois()
self.xdmf_file_sol.write(0.)
+ self.write_vtus = bool(write_vtus)
+ if (self.write_vtus):
+ dcm.write_VTU_file(
+ filebasename=self.write_sol_filebasename,
+ function=self.problem.subsols["U"].subfunc,
+ time=0)
+
+ self.write_xmls = bool(write_xmls)
+ if (self.write_xmls):
+ dolfin.File(self.write_sol_filebasename+"_"+str(0).zfill(3)+".xml") << self.problem.subsols["U"].subfunc
+
def close(self):
@@ -111,6 +124,7 @@ class TimeIntegrator():
def integrate(self):
k_step = 0
+ k_t_tot = 0
n_iter_tot = 0
self.printer.inc()
for step in self.problem.steps:
@@ -173,6 +187,7 @@ class TimeIntegrator():
self.printer.inc()
while (True):
k_t += 1
+ k_t_tot += 1
self.printer.print_var("k_t",k_t,-1)
if (t+dt > step.t_fin):
@@ -247,6 +262,15 @@ class TimeIntegrator():
self.problem.update_fois()
self.xdmf_file_sol.write(t)
+ if (self.write_vtus):
+ dcm.write_VTU_file(
+ filebasename=self.write_sol_filebasename,
+ function=self.problem.subsols["U"].subfunc,
+ time=k_t_tot)
+
+ if (self.write_xmls):
+ dolfin.File(self.write_sol_filebasename+"_"+str(k_t_tot).zfill(3)+".xml") << self.problem.subsols["U"].subfunc
+
if (self.write_qois):
self.problem.update_qois()
self.qoi_printer.write_line([t]+[qoi.value for qoi in self.problem.qois])
@@ -277,6 +301,7 @@ class TimeIntegrator():
constraint.restore_old_value()
k_t -= 1
+ k_t_tot -= 1
t -= dt
dt /= self.decel_coeff
diff --git a/__init__.py b/__init__.py
index 5ac4b35bd38da3be6c7f7d21fabfb807ac5200b8..32445671c8c31bdf9c47e828a80434fe82c8a85f 100644
--- a/__init__.py
+++ b/__init__.py
@@ -13,6 +13,8 @@ from .Material_Elastic_Bulk_Pneumo import *
from .Material_Elastic_Bulk_Poro_Skeleton import *
from .Material_Elastic_Poro_Wpor import *
from .Material_Elastic_Bulk_CiarletGeymonat import *
+from .Material_Elastic_Dev_MooneyRivlin import *
+from .Kinematics_Spherical import *
from .Loading import *
from .Problem_Porosity_Wpor import *
from .Material_Elastic_Porous import *
@@ -26,8 +28,10 @@ from .Material_Elastic_Lung import *
from .Material_Elastic import *
from .Kinematics import *
from .TimeVaryingConstant import *
+from .write_VTU_file import *
from .Step import *
from .SubDomain_Linearized import *
+from .Material_Elastic_CiarletGeymonatNeoHookeanMooneyRivlin import *
from .Problem_RelaxedGrowth import *
from .Problem_Hyperelasticity import *
from .Material_Elastic_Dev_NeoHookean import *
@@ -45,6 +49,8 @@ from .Material_Elastic_Dev_MooneyRivlin import *
# from .Material_Elastic_Dev_MooneyRivlin_reduced import *
from .Constraint import *
from .FOI import *
+from .Problem_Spherical import *
+from .Material_Elastic_Dev_NeoHookeanMooneyRivlin import *
from .Material_Elastic_Bulk import *
from .Material_Elastic_Hooke import *
from .Material_Inelastic_Growth_TimeDriven import *
diff --git a/write_VTU_file.py b/write_VTU_file.py
new file mode 100644
index 0000000000000000000000000000000000000000..b8617efeb7b9f94ad597e817f156fb769f0b3c9e
--- /dev/null
+++ b/write_VTU_file.py
@@ -0,0 +1,31 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2018-2020 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import os
+import shutil
+
+import dolfin_cm as dcm
+
+################################################################################
+
+def write_VTU_file(
+ filebasename,
+ function,
+ time,
+ zfill=3):
+
+ file_pvd = dolfin.File(filebasename+"__.pvd")
+ file_pvd << (function, float(time))
+ os.remove(
+ filebasename+"__.pvd")
+ shutil.move(
+ filebasename+"__"+"".zfill(6)+".vtu",
+ filebasename+"_"+str(time).zfill(zfill)+".vtu")