Thussein mechanical bc

README.md

@@ -1,17 +1,66 @@
 # SurroGlas
 
-## Requirements
+## Introduction
 
-- [```FEniCSx```](https://github.com/FEniCS/dolfinx#installation)
+This is the working repository for the publicly funded research project "SurroGlas".
+It contains the necessary files to run the Finite Element simulation that models float glass tempering.
 
-- [```gmsh```](https://gmsh.info/doc/texinfo/gmsh.html#Installing-and-running-Gmsh-on-your-computer)
+## Structure
 
-- [```NumPy```](https://numpy.org/install/)
+### `main.py`
+
+At the top level, inputs are given in `main.py`.
+From there, you can configure various aspects of the model, such as:
+
+- the problem dimension using `problem_dim`,
+- the solution method (Continuous or Discontinuous Galerkin) and order using `fe_config` and
+- model parameters using the `model_params` dict.
+
+### `ThermoViscoProblem`
+
+This class holds the top-level data structure that includes the sub-models, as well as all data structures relevant for Finite Element analysis.
+
+### `ThermalModel`
+
+Here, all parameters are stored that are required to solve the heat equation.
+
+### `ViscoelasticModel`
+
+This class holds the necessary parameters, expressions and functions that are required to compute derived quantities in the viscoelastic material model, i.e. fictive temperature, shift function, scaled time and the various strain and stress increments.
+
+## Installation
+
+### Local
+
+To install all required packages locally, a working Python installation is required.
+In a corresponding environment, run within a shell:
+
+```bash
+pip3 install -r requirements.txt
+```
+
+### Docker
+
+You might also want to use the pre-built docker container `dolfinx/dolfinx:v0.7.3` that the FEniCS project provides.
+To execute the scripts in the container, clone the git repo and run the following command in a shell:
+
+```bash
+docker run -ti -v $(pwd):/root dolfinx/dolfinx:v0.7.3
+```
+
+This mounts the repository in the root folder of the container, in which a terminal session is opened.
 
 ## Running
 
-From within the local Python environment that the required Python packages are installed, execute
+From within the local Python environment where the required Python packages are installed, execute
 
 ```bash
-python main.py
+python3 main.py
 ```
+to run in serial.
+
+For larger models, if you wish to simulate in parallel using `N` cores, execute
+
+```bash
+mpiexec -np N python3 main.py -parallel
+````

ThermalModel.py

@@ -26,4 +26,4 @@ def __init__(self,mesh: Mesh, model_parameters: dict) -> None:
             # Ambient temperature [K]
             self.T_ambient = Constant(mesh, model_parameters["T_ambient"])
 
-            return
+            return

ThermoViscoProblem.py

@@ -6,7 +6,7 @@
 from dolfinx.io import gmshio
 from dolfinx.fem import (FunctionSpace, Function, Constant,)
 from dolfinx.fem.petsc import NonlinearProblem
-from ufl import (TestFunction, FiniteElement, TensorElement,
+from ufl import (TestFunction,TrialFunction, FiniteElement, TensorElement,
                  VectorElement, grad, inner,
                  CellDiameter, avg, jump,
                  Measure, SpatialCoordinate, FacetNormal)#, ds, dx
@@ -18,16 +18,19 @@
 from time import time
 from ViscoelasticModel import ViscoelasticModel
 from ThermalModel import ThermalModel
-
+from dolfinx import default_scalar_type
+from ufl import (inner, tr, sym, Identity,dot, nabla_div)
+from dolfinx.fem import Constant, Expression
 
 class ThermoViscoProblem:
-    def __init__(self,mesh_path: str,time: tuple, dt: float,
-                 config: dict, model_parameters: dict,
+    def __init__(self,mesh_path: str, problem_dim: int, time: tuple, 
+                 dt: float, config: dict, model_parameters: dict,
                  jit_options: (dict|None) = None ) -> None:
+
+        self.dim = problem_dim
         self.mesh, self.cell_tags, self.facet_tags = gmshio.read_from_msh(
-            mesh_path, MPI.COMM_WORLD, 0, gdim=1)
+            mesh_path, MPI.COMM_WORLD, 0, gdim=problem_dim)
 
-        self.dim = self.mesh.topology.dim
         self.dt = dt
         # The time domain
         self.time = time
@@ -100,6 +103,12 @@ def __init_function_spaces(self,config: dict) -> None:
                                  shape=(self.material_model.tableau_size,self.dim,self.dim))
         self.functionSpaces["sigma_partial"] = FunctionSpace(self.mesh,self.finiteElements["sigma_partial"])
 
+        # Displacements
+        self.finiteElements["U"] = VectorElement(config["U"]["element"],
+                                    self.mesh.ufl_cell(),
+                                    degree=config["U"]["degree"])
+        self.functionSpaces["U"] = FunctionSpace(mesh=self.mesh, element=self.finiteElements["U"])
+
         return
 
 
@@ -169,19 +178,29 @@ def __init_functions(self) -> None:
         self.functions_next["sigma_partial"] = Function(self.functionSpaces["sigma_partial"])
 
         self.functions_next["sigma"] = Function(self.functionSpaces["sigma"], name="Stress_tensor")
+
+        self.functions["U"] = Function(self.functionSpaces["U"], name="Displacement")
+        self.u_trial = TrialFunction(self.functionSpaces["U"])
+        self.v_test = TestFunction(self.functionSpaces["U"])
+
+        self.functions["mech_strain"] = Function(self.functionSpaces["sigma"], name="Mechanical_strain")
+
+
 
         return
 
 
-    def setup(self, dirichlet_bc: bool = False,
-              outfile_name: str = "visco",
-              outfile_name1: str = "stresses") -> None:
+    def setup(self, dirichlet_bc: bool = True,
+              outfile_T: str = "visco",
+              outfile_sigma: str = "stresses") -> None:
         self._set_initial_condition(temp_value=self.material_model.T_init)
         if dirichlet_bc:
-            self._set_dirichlet_bc(bc_value=self.material_model.T_ambient)
+            self._set_dirichlet_bc()
         self._write_initial_output(t=self.t)
-        self._setup_weak_form()
-        self._setup_solver()
+        self._setup_weak_form_T()
+        self._setup_solver_T()
+        self._setup_weak_form_u()
+        self._setup_solver_u()
 
 
     def _set_initial_condition(self, temp_value: float) -> None:
@@ -218,9 +237,6 @@ def __set_IC_Tf_partial(self) -> None:
         values.
         For t0, Tf(n) = T (c.f. Nielsen et al., eq. 27)
         """
-        #for (previous,current) in zip(self.functions_previous["Tf_partial"],self.functions_current["Tf_partial"]):
-        #    current.x.array[:] = self.functions_current["T"].x.array[:]
-        #    previous.x.array[:] = self.functions_previous["T"].x.array[:]
         temp_value = self.functions_current["T"].x.array[0]
         dim = self.material_model.tableau_size
         def Tf_init(x):
@@ -233,15 +249,6 @@ def Tf_init(x):
         return
 
 
-    def _set_dirichlet_bc(self, bc_value: float) -> None:
-        fdim = self.mesh.topology.dim - 1
-        boundary_facets = locate_entities_boundary(
-            self.mesh, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
-        self.bc = fem.dirichletbc(PETSc.ScalarType(bc_value), 
-                                  fem.locate_dofs_topological(self.fs, fdim, boundary_facets), self.fs)
-
-        return
-
 
     def _write_initial_output(self,t: float = 0.0) -> None:
         self.vtx_files = [
@@ -257,9 +264,16 @@ def _write_initial_output(self,t: float = 0.0) -> None:
             # BUG: VTXWriter doesn't support mixed elements
             #io.VTXWriter(self.mesh.comm,"output/Tf_partial.bp",
             #             [self.functions_current["Tf_partial"]],engine="BP4"),
+            # thermal strain
+            io.VTXWriter(self.mesh.comm,"output/eth.bp",
+                         [self.functions["thermal_strain"]],engine="BP4"),
             # Shifted time
             io.VTXWriter(self.mesh.comm,"output/xi.bp",
                          [self.functions["xi"]],engine="BP4"),
+            # Displacements
+            io.VTXWriter(self.mesh.comm,"output/u.bp",
+                         [self.functions["U"]],engine="BP4"),
+
         ]
 
         for file in self.vtx_files:
@@ -275,9 +289,9 @@ def _write_initial_output(self,t: float = 0.0) -> None:
 
         return
 
-        
-
-    def _setup_weak_form(self) -> None:
+    # Heat diffusion equation #
+    
+    def _setup_weak_form_T(self) -> None:
         ds = Measure("exterior_facet",domain=self.mesh)
         dx = Measure("dx",domain=self.mesh)
 
@@ -327,7 +341,7 @@ def _setup_weak_form(self) -> None:
         return
 
 
-    def _setup_solver(self) -> None:
+    def _setup_solver_T(self) -> None:
         self.prob = NonlinearProblem(F=self.F,u=self.functions_current["T"],
                                                jit_options=self.jit_options)
 
@@ -344,7 +358,61 @@ def _setup_solver(self) -> None:
         opts[f"{option_prefix}pc_type"] = "gamg"
         opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
         self.ksp.setFromOptions()
+
+    # linear elasticity equation #
+
 
+    def _set_dirichlet_bc(self) -> None:
+
+        u_bc = fem.Constant(self.mesh, PETSc.ScalarType((0.0,0.0)))
+        facet_dim = self.mesh.topology.dim-1
+        facets_l = locate_entities_boundary(self.mesh, dim=facet_dim, marker=lambda x: np.isclose(x[0], 0.0))
+        facets_r = locate_entities_boundary(self.mesh, dim=facet_dim, marker=lambda x: np.isclose(x[0], 50.0))
+        facets_t = locate_entities_boundary(self.mesh, dim=facet_dim, marker=lambda x: np.isclose(x[1], 0.0))
+        facets_b = locate_entities_boundary(self.mesh, dim=facet_dim, marker=lambda x: np.isclose(x[1], 25.0))
+        dofs_l = fem.locate_dofs_topological(V=self.functionSpaces["U"], entity_dim=facet_dim, entities=facets_l)
+        dofs_r = fem.locate_dofs_topological(V=self.functionSpaces["U"], entity_dim=facet_dim, entities=facets_r)
+        dofs_t = fem.locate_dofs_topological(V=self.functionSpaces["U"], entity_dim=facet_dim, entities=facets_t)
+        dofs_b = fem.locate_dofs_topological(V=self.functionSpaces["U"], entity_dim=facet_dim, entities=facets_b)
+        self.bc = [#fem.dirichletbc(u_bc, dofs_l, self.functionSpaces["U"]),
+                   #fem.dirichletbc(u_bc, dofs_r, self.functionSpaces["U"]),
+                   #fem.dirichletbc(u_bc, dofs_t, self.functionSpaces["U"]),
+                   fem.dirichletbc(u_bc, dofs_b, self.functionSpaces["U"])
+                   ]
+
+    def _setup_weak_form_u(self) -> None:
+
+        ds = Measure("exterior_facet",domain=self.mesh)
+        dx = Measure("dx",domain=self.mesh)
+
+        self.ss = Constant(self.mesh,default_scalar_type((0,-3.0)))        # Body force
+        self.traction = Constant(self.mesh,default_scalar_type((0.0,0.0)))   # traction force
+
+        # Lame's elasticity parameters
+
+        mu = self.material_model.mu 
+        lambda_ = self.material_model.lambda_  
+        I = self.material_model.I
+
+
+        def elastic_epsilon(ua):
+            return sym(grad(ua)) 
+
+        def elastic_sigma(ua):
+            return lambda_ * (nabla_div(ua))* I + 2 * mu * elastic_epsilon(ua)
+
+        self.a = inner(elastic_sigma(self.u_trial), elastic_epsilon(self.v_test)) * dx 
+        self.L = dot(self.ss, self.v_test) * dx + dot(self.traction,self.v_test) * ds
+
+        self.expressions = {}
+        self.expressions["U"] = Expression(elastic_epsilon(self.functions["U"]),
+            self.functionSpaces["sigma"].element.interpolation_points()
+        )
+
+
+    def _setup_solver_u(self) -> None:
+
+        self.problem = fem.petsc.LinearProblem(self.a, self.L, u=self.functions["U"], bcs=self.bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
 
     def _update_values(self,current: Function,previous: Function) -> None:
         # Update ghost values across processes, relevant for MPI computations
@@ -358,25 +426,26 @@ def _write_output(self) -> None:
         for file in self.vtx_files:
             file.write(t=self.t)
 
-        self.outfile_sigma.write_function(self.functions_next["sigma"],
-                                          self.t)
+        self.outfile_sigma.write_function(self.functions_next["sigma"], self.t)
 
         return
 
 
     def solve_timestep(self,t) -> None:
         print(f"t={self.t}")
         self._solve_T()
+        self._solve_u()
         self._solve_Tf()
         self._solve_strains()
         self._solve_shifted_time()
         self._solve_stress()
         self._write_output()
 
-        # For some computations, functions_previous["T"] is needed
+        # For some computations, functions_previous["T"] and functions_previous["displacement"] is needed
         # thus, we update only at the end of each timestep
         self._update_values(current=self.functions_current["T"],
                             previous=self.functions_previous["T"])
+        #self._update_values(current=self.functions_current["displacement"],previous=self.functions_previous["displacement"])
 
         return
 
@@ -390,6 +459,15 @@ def _solve_T(self) -> None:
         assert(converged)
         return
 
+    def _solve_u(self) -> None:
+        """
+        Solve the linear elasticity equation for each time step.
+        Update values and write current values to file.
+        """
+        self.problem.solve()
+
+        return
+
     def _solve_Tf(self) -> None:
         """
         Calculate the current fictive temperature,
@@ -448,6 +526,7 @@ def _solve_stress(self) -> None:
         self.__update_deviatoric_stress()
         self.__update_hydrostatic_stress()
         self.__update_total_stress()
+        self.__interpolation_mechanical_strain()
 
         return
 
@@ -528,10 +607,8 @@ def __update_T_next(self) -> None:
         return
 
     def __update_phi(self) -> None:
-        self.functions["phi"].interpolate(
-            self.material_model.expressions["phi"])
-        self.functions_next["phi"].interpolate(
-            self.material_model.expressions["phi_next"]
+        self.functions["phi"].interpolate(self.material_model.expressions["phi"])
+        self.functions_next["phi"].interpolate(self.material_model.expressions["phi_next"]
         )
 
         return
@@ -593,6 +670,13 @@ def __update_total_stress(self) -> None:
         )
 
         return
+
+    def __interpolation_mechanical_strain(self) -> None:
+        self.functions["mech_strain"].interpolate(
+            self.expressions["U"]
+        )
+
+        return
 
 
     def solve(self) -> None: