[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: master
commit ff6fe3768eab5b0ae7393177aa779079c81ac5e5
Author: Konstantinos Poulios <address@hidden>
Date:   Thu Jun 28 15:42:44 2018 +0200

    Add phase field with elasticity example
---
 interface/tests/python/Makefile.am         |   1 +
 interface/tests/python/demo_phase_field.py | 194 +++++++++++++++++++++++++++++
 2 files changed, 195 insertions(+)

diff --git a/interface/tests/python/Makefile.am 
b/interface/tests/python/Makefile.am
index e9af543..bd3c0ba 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -35,6 +35,7 @@ EXTRA_DIST=                                           \
        demo_laplacian_DG.py                            \
        demo_laplacian_aposteriori.py                   \
        demo_mortar.py                                  \
+       demo_phase_field.py                             \
        demo_plasticity.py                              \
        demo_plate.py                                   \
        demo_static_contact.py                          \
diff --git a/interface/tests/python/demo_phase_field.py 
b/interface/tests/python/demo_phase_field.py
new file mode 100644
index 0000000..a223fd4
--- /dev/null
+++ b/interface/tests/python/demo_phase_field.py
@@ -0,0 +1,194 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2018-2018 Konstantinos Poulios.
+#
+# This file is a part of GetFEM++
+#
+# GetFEM++  is  free software;  you  can  redistribute  it  and/or modify it
+# under  the  terms  of the  GNU  Lesser General Public License as published
+# by  the  Free Software Foundation;  either version 2.1 of the License,  or
+# (at your option) any later version.
+# This program  is  distributed  in  the  hope  that it will be useful,  but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+# License for more details.
+# You  should  have received a copy of the GNU Lesser General Public License
+# along  with  this program;  if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
+#
+############################################################################
+
+import getfem as gf
+import numpy as np
+import os, sys
+
+np.set_printoptions(threshold=np.nan)
+gf.util_trace_level(1)
+
+# Input data
+NX = 60         # number of elements in horizontal direction
+NY = 60         # number of elements in verical direction
+
+LX = 1.         # [mm] Length
+LY = 1.         # [mm] Height
+
+E = 210e3       # [N/mm^2]
+nu = 0.3
+
+Gc0 = 2.7       # [N/mm]
+
+ll = 0.01       # [mm] length scale
+
+strain_rate = 2e-4
+init_time_step = 1.
+
+steps = 100
+
+#------------------------------------
+geotrans = "GT_QK(2,2)"  # geometric transformation
+disp_fem_order = 2       # displacements finite element order
+phi_fem_order = 2        # phase field finite element order
+
+integration_degree = 5   # 9 gauss points per quad
+
+#------------------------------------
+
+resultspath = "./demo_phase_field_results"
+if not os.path.exists(resultspath):
+   os.makedirs(resultspath)
+
+# auxiliary constants
+B_BOUNDARY = 4
+T_BOUNDARY = 6
+TB_BOUNDARY = 7
+
+NX_seed1 = np.linspace(-1, 0, NX/3+1)
+NX_seed2 = np.linspace(0, 1, (NX-NX/3)+1)[1:]
+NY_seed = np.linspace(-1, 1, NY+1)
+X_seed1 = 
LX/2*(0.2*NX_seed1+0.8*np.sign(NX_seed1)*np.power(np.abs(NX_seed1),1.5))
+X_seed2 = 
LX/2*(0.6*NX_seed2+0.4*np.sign(NX_seed2)*np.power(np.abs(NX_seed2),1.5))
+X_seed = np.concatenate((X_seed1,X_seed2))
+Y_seed = LY/2*(0.2*NY_seed+0.8*np.sign(NY_seed)*np.power(np.abs(NY_seed),1.5))
+mesh = gf.Mesh("cartesian", X_seed, Y_seed)
+N = mesh.dim()
+
+bottom_faces = 
mesh.outer_faces_in_box([-LX/2-1e-5,-LY/2-1e-5],[LX/2+1e-5,-LY/2+1e-5])
+top_faces = 
mesh.outer_faces_in_box([-LX/2-1e-5,LY/2-1e-5],[LX/2+1e-5,LY/2+1e-5])
+
+mesh.set_region(B_BOUNDARY, bottom_faces)
+mesh.set_region(T_BOUNDARY, top_faces)
+
+mesh.region_merge(TB_BOUNDARY, T_BOUNDARY)
+mesh.region_merge(TB_BOUNDARY, B_BOUNDARY)
+
+mesh.export_to_vtk("%s/mesh.vtk" % resultspath)
+
+# FEM
+mfu = gf.MeshFem(mesh, N)
+mfu.set_classical_fem(disp_fem_order)
+
+mfdir = mfu
+
+mfphi = gf.MeshFem(mesh, 1)
+mfphi.set_classical_fem(phi_fem_order)
+
+mfout = gf.MeshFem(mesh)
+mfout.set_classical_discontinuous_fem(2)
+
+# Integration method
+mim = gf.MeshIm(mesh, integration_degree)
+mimd1 = gf.MeshImData(mim, -1)
+
+# Model
+md = gf.Model("real")
+
+md.add_fem_variable("u", mfu)          # displacements
+md.add_fem_variable("phi", mfphi)      # phase field
+md.set_variable("phi", np.ones(mfphi.nbdof()))
+
+md.add_fem_data("u_stored", mfu)
+
+md.add_im_data("psi0_max", mimd1)
+
+md.add_initialized_data("kappa", E/(3.*(1.-2.*nu)))
+md.add_initialized_data("mu", E/(2*(1+nu)))
+
+md.add_initialized_data("Gc0", Gc0)
+md.add_initialized_data("l", ll)
+
+md.set_variable("psi0_max", 
md.interpolation("100000*min(1,max(0,1-(abs(X(2))+pos_part(X(1)))/(l/10)))", 
mimd1, -1))
+
+md.add_macro("damage", "max(1e-7,sqr(1-phi))")
+md.add_macro("psi0", 
"(0.5*kappa*sqr(Trace(Grad_u))+mu*Norm_sqr(Deviator(Sym(Grad_u))))")
+md.add_macro("Gc", "Gc0")
+
+md.add_nonlinear_term(mim, 
"damage*(kappa*Trace(Grad_u)*Id(2)+2*mu*Deviator(Sym(Grad_u))):Grad_Test_u")
+md.add_nonlinear_term(mim, 
"(-2*(1-phi)*max(psi0_max,psi0)*Test_phi+Gc*(phi/l*Test_phi+l*Grad_phi.Grad_Test_phi))")
+
+md.add_initialized_data("MM", 0)
+md.add_nonlinear_generic_assembly_brick(mim, "MM*(u-u_stored).Test_u")
+
+# Load
+md.add_fem_data("dirichlet_data", mfu);
+md.add_Dirichlet_condition_with_multipliers(mim, "u", mfdir, TB_BOUNDARY, 
"dirichlet_data")
+
+print("Displacement dofs: %i" % mfu.nbdof())
+print("Total model dofs: %i" % md.nbdof())
+
+time_step = init_time_step
+eps = 0.
+pseudodynamic = False
+MM = 0.
+for step in range(steps):
+  eps_old = eps
+  while True:
+    if step > 0:
+      eps = eps_old + strain_rate*time_step
+
+    X = md.from_variables()
+    print("Step %i with eps=%e" % (step, eps))
+    md.set_variable("dirichlet_data", 
md.interpolation("{eps}*[0;X(2)]".format(eps=eps), mfu, -1))
+    if pseudodynamic:
+      print("With damping %e" % MM)
+    md.set_variable("MM",[MM])
+
+    iters = 20
+    nit = \
+    md.solve("noisy", "lsolver", "mumps", "max_iter", iters, "max_res", 1e-7, 
#)[0]
+             "lsearch", "simplest", "alpha max ratio", 1.5, "alpha min", 0.4, 
"alpha mult", 0.6,
+              "alpha threshold res", 5e3)[0]
+    if nit >= iters:
+      if step == 0:
+        break
+      md.to_variables(X)
+      eps = eps_old
+      if time_step > init_time_step/50.:
+        time_step *= 0.5
+      else:
+        if pseudodynamic:
+          MM *= 100.
+        else:
+          pseudodynamic = True
+          MM = 10000.
+    else:
+      md.set_variable("u_stored", md.variable("u"))
+      if pseudodynamic:
+        if MM < 1e-4:
+          MM = 0.
+          pseudodynamic = False
+        else:
+          MM /= 10.
+      else:
+        TMP = md.interpolation("max(psi0_max,psi0)", mimd1, -1)
+        md.set_variable("psi0_max", TMP)
+        break
+
+  mfout.export_to_vtk("%s/demo_phase_field_%i.vtk" % (resultspath, step),
+                      mfu, md.variable("u"), "Displacements",
+                      mfphi, md.variable("phi"), "phi")
+
+
+
+