demo_pantograph_full_scale.py

# coding: utf8
"""
Created on 18/04/2018
@author: baptiste

"""
import logging
from pathlib import Path

import dolfin as fe
import gmsh

import ho_homog.geometry as geo
from ho_homog import GEO_TOLERANCE, pckg_logger
from ho_homog.mesh_generate import pantograph
from ho_homog.full_scale_pb import FullScaleModel
from ho_homog.materials import Material
from ho_homog.part import FenicsPart
from ho_homog.periodicity import PeriodicDomain
from ho_homog.toolbox_gmsh import process_gmsh_log

logger = logging.getLogger("demo_full_scale")
logger_root = logging.getLogger()
logger.setLevel(logging.INFO)
logger_root.setLevel(logging.INFO)

for hdlr in pckg_logger.handlers[:]:
    if isinstance(hdlr, logging.StreamHandler):
        pckg_logger.removeHandler(hdlr)
pckg_logger.setLevel(logging.INFO)
formatter = logging.Formatter(
    "%(asctime)s :: %(levelname)s :: %(name)s :: %(message)s", "%H:%M:%S"
)
stream_handler = logging.StreamHandler()
stream_handler.setLevel(logging.INFO)
stream_handler.setFormatter(formatter)
logger_root.addHandler(stream_handler)


# * Step 1 : Modeling the geometry of the part
geo.init_geo_tools()
a = 1
b, k = a, a / 3
r = 0.05
gmsh.logger.start()
panto_geo = pantograph.pantograph_RVE(
    a, b, k, 0.1, nb_cells=(10, 1), soft_mat=True, name="panto_with_soft"
)
process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop()


# * Step 2 : Generating the mesh
lc_ratio = 1 / 3
d_min_max = (2 * r * a, a)
lc_min_max = (lc_ratio * r * a, lc_ratio * a)
panto_geo.main_mesh_refinement((0.1, 0.5), (0.1, 0.3), False)
panto_geo.soft_mesh_refinement((0.1, 0.5), (0.1, 0.3), False)
gmsh.logger.start()
panto_geo.mesh_generate()
process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop()


# * Step 2 : Defining the material properties
E1, nu1 = 1.0, 0.3
E2, nu2 = E1 / 100.0, nu1
E_nu_tuples = [(E1, nu1), (E2, nu2)]

subdo_tags = tuple([subdo.tag for subdo in panto_geo.phy_surf])
material_dict = dict()
for coeff, tag in zip(E_nu_tuples, subdo_tags):
    material_dict[tag] = Material(coeff[0], coeff[1], "cp")

# * Step 3 : Create part object
panto_part = FenicsPart.file_2_FenicsPart(
    panto_geo.mesh_abs_path, material_dict, panto_geo.gen_vect, subdomains_import=True
)
LX = panto_part.global_dimensions[0, 0]


# * Step 4 : Defining boundary conditions
class LeftBorder(fe.SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] <= GEO_TOLERANCE


class RightBorder(fe.SubDomain):
    def inside(self, x, on_boundary):
        right = x[0] >= LX - GEO_TOLERANCE
        return on_boundary and right


l_border = LeftBorder()
r_border = RightBorder()
boundary_markers = fe.MeshFunction("size_t", panto_part.mesh, dim=panto_part.dim - 1)
boundary_markers.set_all(9999)
l_border.mark(boundary_markers, 11)
r_border.mark(boundary_markers, 99)

pbc = PeriodicDomain.pbc_dual_base(panto_part.global_dimensions, "Y")
boundary_conditions = [
    {"type": "Periodic", "constraint": pbc},
    {
        "type": "Dirichlet",
        "constraint": (0.0, 0.0),
        "facet_function": boundary_markers,
        "facet_idx": 11,
    },
    ("Dirichlet", (0.0, 0.0), boundary_markers, 99),
]

# * Step 5 : Defining the load
indicator_fctn = fe.Expression(
    "x[0] >= (LX-Lx_cell)/2.0 && x[0] <= (LX+Lx_cell)/2.0 ? 1 : 0",
    degree=0,
    LX=LX,
    Lx_cell=4.0 * a,
)
load_area = fe.assemble(indicator_fctn * fe.dx(panto_part.mesh))
load_magnitude = 1.0
s_load = load_magnitude / load_area
s_load = fe.Constant((0.0, s_load))  # forme utilisée dans tuto FEniCS linear elasticity
load = (2, s_load, indicator_fctn)

# * Step 6 : Choosing FE characteristics
element = ("Lagrange", 2)

# * Step 7 : Gathering all data in a model
model = FullScaleModel(panto_part, [load], boundary_conditions, element)
model.set_solver("mumps")

# * Step 8 : Solving problem
model.solve(results_file=Path("demo_results/demo_pantograph_full_scale.xdmf"))