part_test.py 14.7 KB
Newer Older
1 2 3 4 5 6 7
# coding: utf-8
"""
Created on 09/01/2019
@author: baptiste

"""
import logging
8 9 10 11
from pathlib import Path
from subprocess import run

import dolfin as fe
12
import gmsh
13
import mshr
14
import numpy as np
15
from pytest import approx
16

17 18 19 20 21
from ho_homog import geometry as geo
from ho_homog import materials as mat
from ho_homog import mesh_generate as mg
from ho_homog import mesh_tools
from ho_homog import part
22 23 24

model = gmsh.model
factory = model.occ
25
logger = logging.getLogger(__name__)  # http://sametmax.com/ecrire-des-logs-en-python/
26
logger.setLevel(logging.INFO)
27

28 29 30
if __name__ == "__main__":
    logger_root = logging.getLogger()
    logger_root.setLevel(logging.INFO)
31 32 33
    logging.basicConfig(
        format="%(asctime)s :: %(levelname)s :: %(message)s", level=logging.INFO
    )
34

35 36 37 38 39 40 41 42 43 44 45
geo.init_geo_tools()
geo.set_gmsh_option("General.Verbosity", 3)

# ? Test du constructeur gmsh_2_Fenics_2DRVE
def test_rve_from_gmsh2drve():
    a = 1
    b, k = a, a / 3
    panto_test = mg.pantograph.pantograph_RVE(
        a, b, k, 0.1, nb_cells=(2, 3), soft_mat=False, name="panto_test"
    )
    panto_test.main_mesh_refinement((0.1, 0.5), (0.03, 0.3), False)
46

47 48
    panto_test.mesh_generate()
    run(f"gmsh {panto_test.name}.msh &", shell=True, check=True)
49

50 51 52 53 54 55 56 57 58 59 60 61 62
    E1, nu1 = 1.0, 0.3
    E2, nu2 = E1 / 100.0, nu1
    E_nu_tuples = [(E1, nu1), (E2, nu2)]
    phy_subdomains = panto_test.phy_surf
    material_dict = dict()
    for coeff, subdo in zip(E_nu_tuples, phy_subdomains):
        material_dict[subdo.tag] = mat.Material(coeff[0], coeff[1], "cp")
    rve = part.Fenics2DRVE.rve_from_gmsh2drve(panto_test, material_dict)
    assert isinstance(rve, part.Fenics2DRVE)
    assert hasattr(rve, "mesh")
    assert hasattr(rve, "gen_vect")
    assert hasattr(rve, "C_per")
    assert hasattr(rve, "rve_area")
63

64 65 66 67 68

# test_rve_from_gmsh2drve()
# * Test ok

# ? Problème de maillage lorsqu'un domaine mou existe
69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
# a = 1
# b, k, r = a, a/3, a/10
# panto_test = prt.Gmsh2DRVE.pantograph(a, b, k, r, soft_mat=True, name='panto_test_soft')
# run(f"gmsh {panto_test.name}.brep &", shell=True, check=True)
# panto_test.main_mesh_refinement((r, 5*r), (r/2, a/3), False)
# panto_test.soft_mesh_refinement((r, 5*r), (r/2, a/3), False)
# panto_test.mesh_generate() #!Sans l'opération de removeDuplicateNodes() dans mesh_generate(). Avec domaine mou obtenu par intersection.
# gmsh.model.mesh.removeDuplicateNodes()
# gmsh.write(f"{panto_test.name}_clean.msh")
# run(f"gmsh {panto_test.name}.msh &", shell=True, check=True)
# run(f"gmsh {panto_test.name}_clean.msh &", shell=True, check=True)
# a = 1
# b, k, r = a, a/3, a/10
# panto_test = prt.Gmsh2DRVE.pantograph(a, b, k, r, soft_mat=True, name='panto_test_soft')
# run(f"gmsh {panto_test.name}.brep &", shell=True, check=True)
# panto_test.main_mesh_refinement((r, 5*r), (r/2, a/3), False)
# panto_test.soft_mesh_refinement((r, 5*r), (2*r, a/3), False)
# panto_test.mesh_generate() #!Sans l'opération de removeDuplicateNodes() dans mesh_generate(). Avec domaine mou obtenu par intersection.
# gmsh.write(f"{panto_test.name}_clean.msh")
# run(f"gmsh {panto_test.name}.msh &", shell=True, check=True)
# run(f"gmsh {panto_test.name}_clean.msh &", shell=True, check=True)
90
# ? Fin de la résolution
91 92


93
# ? Test d'un mesh avec un domaine mou, gmsh2DRVE puis import dans FEniCS
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
# a = 1
# b, k = a, a/3
# panto_test = mesh_generate_2D.Gmsh2DRVE.pantograph(a, b, k, 0.1, nb_cells=(2, 3), soft_mat=True, name='panto_test_mou')
# panto_test.main_mesh_refinement((0.1,0.5),(0.03,0.3),False)
# panto_test.soft_mesh_refinement((0.1,0.5),(0.03,0.3),False)
# panto_test.mesh_generate()
# run(f"gmsh {panto_test.name}.msh &",shell=True, check=True)
# E1, nu1 = 1., 0.3
# E2, nu2 = E1/100., nu1
# E_nu_tuples = [(E1, nu1), (E2, nu2)]
# phy_subdomains = panto_test.phy_surf
# material_dict = dict()
# for coeff, subdo in zip(E_nu_tuples, phy_subdomains):
#     material_dict[subdo.tag] = mat.Material(coeff[0], coeff[1], 'cp')
# rve = prt.Fenics2DRVE.gmsh_2_Fenics_2DRVE(panto_test, material_dict)
# # #* Test ok


# L, t = 1, 0.05
# a = L-3*t
# aux_test = mesh_generate_2D.Gmsh2DRVE.auxetic_square(L, a, t, nb_cells=(2, 3), soft_mat=True, name='aux_test_mou')
# aux_test.main_mesh_refinement((0.1,0.5),(0.03,0.3),False)
# aux_test.soft_mesh_refinement((0.1,0.5),(0.03,0.3),False)
# aux_test.mesh_generate()
# run(f"gmsh {aux_test.name}.msh &",shell=True, check=True)
# E1, nu1 = 1., 0.3
# E2, nu2 = E1/100., nu1
# E_nu_tuples = [(E1, nu1), (E2, nu2)]
# phy_subdomains = aux_test.phy_surf
# material_dict = dict()
# for coeff, subdo in zip(E_nu_tuples, phy_subdomains):
#     material_dict[subdo.tag] = mat.Material(coeff[0], coeff[1], 'cp')
# rve = prt.Fenics2DRVE.gmsh_2_Fenics_2DRVE(aux_test, material_dict)

128

129 130
def test_global_area_2D():
    """FenicsPart method global_aera, for a 2D part"""
131
    L_x, L_y = 10.0, 4.0
132
    mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(L_x, L_y), 10, 10)
133 134 135 136 137 138 139 140 141 142 143
    material = {"0": mat.Material(1.0, 0.3, "cp")}
    dimensions = np.array(((L_x, 0.0), (0.0, L_y)))
    rect_part = part.FenicsPart(
        mesh,
        materials=material,
        subdomains=None,
        global_dimensions=dimensions,
        facet_regions=None,
    )
    assert rect_part.global_area == approx(40.0, rel=1e-10)

144 145 146

def test_mat_area():
    """FenicsPart method global_aera, for 2D parts created from a gmsh mesh and from a FEniCS mesh"""
147 148
    L_x, L_y = 4.0, 5.0
    H = 1.0
149
    size = 0.5
150 151 152 153 154
    rectangle = mshr.Rectangle(fe.Point(0.0, 0), fe.Point(L_x, L_y))
    hole = mshr.Rectangle(
        fe.Point(L_x / 2 - H / 2, L_y / 2 - H / 2),
        fe.Point(L_x / 2 + H / 2, L_y / 2 + H / 2),
    )
155 156 157
    domain = rectangle - hole
    domain.set_subdomain(1, rectangle)
    mesh = mshr.generate_mesh(domain, size)
158 159 160 161 162 163 164 165 166 167
    dimensions = np.array(((L_x, 0.0), (0.0, L_y)))
    material = {"0": mat.Material(1.0, 0.3, "cp")}
    rect_part = part.FenicsPart(
        mesh,
        materials=material,
        subdomains=None,
        global_dimensions=dimensions,
        facet_regions=None,
    )
    assert rect_part.mat_area == (L_x * L_y - H ** 2)
168 169 170

    name = "test_mat_area"
    local_dir = Path(__file__).parent
171
    mesh_file = local_dir.joinpath(name + ".msh")
172
    gmsh.model.add(name)
173
    vertices = [(0.0, 0.0), (0.0, L_y), (L_x, L_y), (L_x, 0.0)]
174 175 176
    contour = geo.LineLoop([geo.Point(np.array(c)) for c in vertices], False)
    surface = geo.PlaneSurface(contour)
    cut_vertices = list()
177 178
    for local_coord in [(H, 0.0, 0.0), (0.0, H, 0.0), (-H, 0.0, 0.0), (0.0, -H, 0.0)]:
        vertex = geo.translation(contour.vertices[2], np.array(local_coord))
179
        cut_vertices.append(vertex)
180
    cut_surface = geo.PlaneSurface(geo.LineLoop(cut_vertices, False))
181 182 183
    for s in [surface, cut_surface]:
        s.add_gmsh()
    factory.synchronize()
184
    (surface,) = geo.surface_bool_cut(surface, cut_surface)
185 186 187 188
    factory.synchronize()
    for dim_tag in model.getEntities(2):
        if not dim_tag[1] == surface.tag:
            model.removeEntities(dim_tag, True)
189 190 191
    charact_field = mesh_tools.MathEvalField("0.1")
    mesh_tools.set_background_mesh(charact_field)
    geo.set_gmsh_option("Mesh.SaveAll", 1)
192 193 194 195
    model.mesh.generate(2)
    gmsh.write(str(mesh_file))
    cmd = f"dolfin-convert {mesh_file} {mesh_file.with_suffix('.xml')}"
    run(cmd, shell=True, check=True)
196 197 198 199 200 201 202 203 204 205 206
    mesh = fe.Mesh(str(mesh_file.with_suffix(".xml")))
    dimensions = np.array(((L_x, 0.0), (0.0, L_y)))
    material = {"0": mat.Material(1.0, 0.3, "cp")}
    rect_part = part.FenicsPart(
        mesh,
        materials=material,
        subdomains=None,
        global_dimensions=dimensions,
        facet_regions=None,
    )
    assert rect_part.mat_area == approx(L_x * L_y - H * H / 2)
207 208
    geo.reset()

209

210 211 212 213
def test_mat_without_dictionnary():
    """FenicsPart instance initialized with only one instance of Material"""
    L_x, L_y = 1, 1
    mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(L_x, L_y), 10, 10)
214
    dimensions = np.array(((L_x, 0.0), (0.0, L_y)))
215
    E, nu = 1, 0.3
216 217 218 219 220 221 222 223 224
    material = mat.Material(E, nu, "cp")
    rect_part = part.FenicsPart(
        mesh,
        materials=material,
        subdomains=None,
        global_dimensions=dimensions,
        facet_regions=None,
    )
    elem_type = "CG"
225
    degree = 2
226 227 228 229 230 231 232
    strain_fspace = fe.FunctionSpace(
        mesh, fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3),
    )
    strain = fe.project(
        fe.Expression(("1.0+x[0]*x[0]", "0", "1.0"), degree=2), strain_fspace
    )
    stress = mat.sigma(rect_part.elasticity_tensor, strain)
233
    energy = fe.assemble(fe.inner(stress, strain) * fe.dx(rect_part.mesh))
234
    energy_theo = E / (1 + nu) * (1 + 28 / (15 * (1 - nu)))
235 236
    assert energy == approx(energy_theo, rel=1e-13)

237

238 239 240 241
def test_2_materials():
    """FenicsPart instance initialized with only one instance of Material in the materials dictionnary"""
    L_x, L_y = 1, 1
    mesh = fe.RectangleMesh(fe.Point(-L_x, -L_y), fe.Point(L_x, L_y), 20, 20)
242
    dimensions = np.array(((2 * L_x, 0.0), (0.0, 2 * L_y)))
243
    subdomains = fe.MeshFunction("size_t", mesh, 2)
244

245 246 247
    class Right_part(fe.SubDomain):
        def inside(self, x, on_boundary):
            return x[0] >= 0 - fe.DOLFIN_EPS
248

249 250 251 252
    subdomain_right = Right_part()
    subdomains.set_all(0)
    subdomain_right.mark(subdomains, 1)
    E_1, E_2, nu = 1, 3, 0.3
253
    materials = {0: mat.Material(1, 0.3, "cp"), 1: mat.Material(3, 0.3, "cp")}
254
    rect_part = part.FenicsPart(mesh, materials, subdomains, dimensions)
255
    elem_type = "CG"
256 257
    degree = 2
    strain_fspace = fe.FunctionSpace(
258 259 260 261 262 263
        mesh, fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3),
    )
    strain = fe.project(
        fe.Expression(("1.0+x[0]*x[0]", "0", "1.0"), degree=2), strain_fspace
    )
    stress = mat.sigma(rect_part.elasticity_tensor, strain)
264
    energy = fe.assemble(fe.inner(stress, strain) * fe.dx(rect_part.mesh))
265
    energy_theo = 2 * ((E_1 + E_2) / (1 + nu) * (1 + 28 / (15 * (1 - nu))))
266 267
    assert energy == approx(energy_theo, rel=1e-13)

268

269 270 271 272
def test_get_domains_gmsh(plots=False):
    """ Get subdomains and partition of the boundary from a .msh file. """
    name = "test_domains"
    local_dir = Path(__file__).parent
273
    mesh_file = local_dir.joinpath(name + ".msh")
274
    gmsh.model.add(name)
275 276 277
    L_x, L_y = 2.0, 2.0
    H = 1.0
    vertices = [(0.0, 0.0), (0.0, L_y), (L_x, L_y), (L_x, 0.0)]
278 279 280
    contour = geo.LineLoop([geo.Point(np.array(c)) for c in vertices], False)
    surface = geo.PlaneSurface(contour)
    inclusion_vertices = list()
281 282 283 284 285 286 287
    for coord in [
        (H / 2, -H / 2, 0.0),
        (H / 2, H / 2, 0.0),
        (-H / 2, H / 2, 0.0),
        (-H / 2, -H / 2, 0.0),
    ]:
        vertex = geo.translation(geo.Point((L_x / 2, L_y / 2)), coord)
288
        inclusion_vertices.append(vertex)
289
    inclusion = geo.PlaneSurface(geo.LineLoop(inclusion_vertices, False))
290 291 292
    for s in [surface, inclusion]:
        s.add_gmsh()
    factory.synchronize()
293
    (stiff_s,) = geo.surface_bool_cut(surface, inclusion)
294
    factory.synchronize()
295
    (soft_s,) = geo.surface_bool_cut(surface, stiff_s)
296 297
    factory.synchronize()
    domains = {
298 299
        "stiff": geo.PhysicalGroup(stiff_s, 2),
        "soft": geo.PhysicalGroup(soft_s, 2),
300 301
    }
    boundaries = {
302 303 304 305
        "S": geo.PhysicalGroup(surface.ext_contour.sides[0], 1),
        "W": geo.PhysicalGroup(surface.ext_contour.sides[1], 1),
        "N": geo.PhysicalGroup(surface.ext_contour.sides[2], 1),
        "E": geo.PhysicalGroup(surface.ext_contour.sides[3], 1),
306 307 308 309 310
    }
    for group in domains.values():
        group.add_gmsh()
    for group in boundaries.values():
        group.add_gmsh()
311 312 313
    charact_field = mesh_tools.MathEvalField("0.05")
    mesh_tools.set_background_mesh(charact_field)
    geo.set_gmsh_option("Mesh.SaveAll", 0)
314 315 316 317 318 319
    model.mesh.generate(1)
    model.mesh.generate(2)
    gmsh.model.mesh.removeDuplicateNodes()
    gmsh.write(str(mesh_file))
    E_1, E_2, nu = 1, 3, 0.3
    materials = {
320 321
        domains["soft"].tag: mat.Material(E_1, nu, "cp"),
        domains["stiff"].tag: mat.Material(E_1, nu, "cp"),
322
    }
323 324 325 326 327
    test_part = part.FenicsPart.part_from_file(
        mesh_file, materials, subdomains_import=True
    )
    assert test_part.mat_area == approx(L_x * L_y)
    elem_type = "CG"
328 329 330 331 332
    degree = 2
    V = fe.VectorFunctionSpace(test_part.mesh, elem_type, degree)
    W = fe.FunctionSpace(
        test_part.mesh,
        fe.VectorElement(elem_type, test_part.mesh.ufl_cell(), degree, dim=3),
333
    )
334
    boundary_conditions = {
335 336 337 338 339
        boundaries["N"].tag: fe.Expression(("x[0]-1", "1"), degree=1),
        boundaries["S"].tag: fe.Expression(("x[0]-1", "-1"), degree=1),
        boundaries["E"].tag: fe.Expression(("1", "x[1]-1"), degree=1),
        boundaries["W"].tag: fe.Expression(("-1", "x[1]-1"), degree=1),
    }
340 341 342
    bcs = list()
    for tag, val in boundary_conditions.items():
        bcs.append(fe.DirichletBC(V, val, test_part.facet_regions, tag))
343
    ds = fe.Measure("ds", domain=test_part.mesh, subdomain_data=test_part.facet_regions)
344 345
    v = fe.TestFunctions(V)
    u = fe.TrialFunctions(V)
346 347 348 349
    F = (
        fe.inner(mat.sigma(test_part.elasticity_tensor, mat.epsilon(u)), mat.epsilon(v))
        * fe.dx
    )
350 351 352
    a, L = fe.lhs(F), fe.rhs(F)
    u_sol = fe.Function(V)
    fe.solve(a == L, u_sol, bcs)
353
    strain = fe.project(mat.epsilon(u_sol), W)
354 355
    if plots:
        import matplotlib.pyplot as plt
356

357 358 359 360 361 362 363 364 365 366 367 368 369
        plt.figure()
        plot = fe.plot(u_sol)
        plt.colorbar(plot)
        plt.figure()
        plot = fe.plot(strain[0])
        plt.colorbar(plot)
        plt.figure()
        plot = fe.plot(strain[1])
        plt.colorbar(plot)
        plt.figure()
        plot = fe.plot(strain[2])
        plt.colorbar(plot)
        plt.show()
370 371 372 373 374 375 376
    error = fe.errornorm(
        strain,
        fe.Expression(("1", "1", "0"), degree=0),
        degree_rise=3,
        mesh=test_part.mesh,
    )
    assert error == approx(0, abs=1e-12)
377

378
    materials = {
379 380
        domains["soft"].tag: mat.Material(E_1, nu, "cp"),
        domains["stiff"].tag: mat.Material(E_2, nu, "cp"),
381
    }
382 383 384
    test_part = part.FenicsPart.part_from_file(
        mesh_file, materials, subdomains_import=True
    )
385 386 387 388 389 390 391 392 393 394
    V = fe.VectorFunctionSpace(test_part.mesh, elem_type, degree)
    W = fe.FunctionSpace(
        test_part.mesh,
        fe.VectorElement(elem_type, test_part.mesh.ufl_cell(), degree, dim=3),
    )
    bcs = list()
    for tag, val in boundary_conditions.items():
        bcs.append(fe.DirichletBC(V, val, test_part.facet_regions, tag))
    v = fe.TestFunctions(V)
    u = fe.TrialFunctions(V)
395 396 397 398
    F = (
        fe.inner(mat.sigma(test_part.elasticity_tensor, mat.epsilon(u)), mat.epsilon(v))
        * fe.dx
    )
399 400 401
    a, L = fe.lhs(F), fe.rhs(F)
    u_sol = fe.Function(V)
    fe.solve(a == L, u_sol, bcs)
402 403
    strain = mat.epsilon(u_sol)
    stress = mat.sigma(test_part.elasticity_tensor, strain)
404 405 406
    energy = 0.5 * fe.assemble(fe.inner(stress, strain) * fe.dx(test_part.mesh))
    energy_abaqus = 12.8788939
    assert energy == approx(energy_abaqus, rel=1e-3)
407
    geo.reset()