demo_pantograph_full_scale.py 3.79 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
# coding: utf8
"""
Created on 18/04/2018
@author: baptiste

"""
import logging
from pathlib import Path

import dolfin as fe
import gmsh

13 14
import ho_homog.geometry as geo
from ho_homog import GEO_TOLERANCE, pckg_logger
15
from ho_homog.mesh_generate import pantograph
16 17 18 19 20
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
21

22
logger = logging.getLogger("demo_full_scale")
23 24 25 26 27 28 29 30
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)
31 32 33
formatter = logging.Formatter(
    "%(asctime)s :: %(levelname)s :: %(name)s :: %(message)s", "%H:%M:%S"
)
34 35 36 37 38
stream_handler = logging.StreamHandler()
stream_handler.setLevel(logging.INFO)
stream_handler.setFormatter(formatter)
logger_root.addHandler(stream_handler)

39 40

# * Step 1 : Modeling the geometry of the part
41
geo.init_geo_tools()
42
a = 1
43
b, k = a, a / 3
44 45
r = 0.05
gmsh.logger.start()
46
panto_geo = pantograph.pantograph_RVE(
47 48
    a, b, k, 0.1, nb_cells=(10, 1), soft_mat=True, name="panto_with_soft"
)
49 50 51 52
process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop()


53 54 55 56 57 58
# * 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)
59 60 61 62 63 64
gmsh.logger.start()
panto_geo.mesh_generate()
process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop()


65 66 67
# * Step 2 : Defining the material properties
E1, nu1 = 1.0, 0.3
E2, nu2 = E1 / 100.0, nu1
68 69 70 71 72
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):
73
    material_dict[tag] = Material(coeff[0], coeff[1], "cp")
74

75
# * Step 3 : Create part object
76 77
panto_part = FenicsPart.file_2_FenicsPart(
    panto_geo.mesh_abs_path, material_dict, panto_geo.gen_vect, subdomains_import=True
78
)
79 80 81
LX = panto_part.global_dimensions[0, 0]


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


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


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

101
pbc = PeriodicDomain.pbc_dual_base(panto_part.global_dimensions, "Y")
102
boundary_conditions = [
103 104 105 106 107 108 109 110
    {"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),
111 112
]

113
# * Step 5 : Defining the load
114
indicator_fctn = fe.Expression(
115 116 117 118 119 120 121
    "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
122
s_load = load_magnitude / load_area
123
s_load = fe.Constant((0.0, s_load))  # forme utilisée dans tuto FEniCS linear elasticity
124 125
load = (2, s_load, indicator_fctn)

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

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

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