Merge branch 'master' of gitlab.inria.fr:mgenet/dolfin_cm

Kinematics.py

@@ -14,7 +14,6 @@ import dolfin
 import numpy
 
 import dolfin_cm as dcm
-from .Problem import Problem
 
 ################################################################################
 

Kinematics_Spherical.py

+#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

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)
 

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)
 

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)
 
 
 

Material_Elastic_CiarletGeymonatNeoHookeanMooneyRivlin.py

+#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

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                                   ###
 ###                                                                          ###

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

Material_Elastic_Dev_NeoHookeanMooneyRivlin.py

+#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

Material_Elastic_Jerome.py

@@ -2,7 +2,7 @@
 
 ################################################################################
 ###                                                                          ###
-### Created by Martin Genet, 2018-2019                                       ###
+### Created by Martin Genet, 2018-2020                                       ###
 ###                                                                          ###
 ### École Polytechnique, Palaiseau, France                                   ###
 ###                                                                          ###

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,

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(

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]],

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

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

Problem_Spherical.py

+#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)

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,

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

__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 *

write_VTU_file.py

+#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")