Skip to content

GitLab

HO_homog
HO_homog
Commit d350fef5 authored by Baptiste Durand's avatar Baptiste Durand
Browse files
Options

Merge branch 'reorganize_part_classes' 

Improve structure of the code.
Better separation between the class definitions on one hand and the imports of the mesh and the subdomain data on the other hand.
parents 02a6ef7a 99c53448
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 646 additions and 3100 deletions
.gitignore
View file @ d350fef5
......@@ -8,4 +8,5 @@
demo/loc_E12_u.pdf
*.pyc.*
*.egg-info/
ho_homog/.idea
\ No newline at end of file
ho_homog/.idea
.idea/HO_homog.iml
demo/demo_pantograph_mesh.py
View file @ d350fef5
......@@ -153,7 +153,7 @@ msh.set_periodicity_pairs(micro_bndry[1], micro_bndry[3])
logger.info("Periodicity constraint : Done")
logger.info("Cleaning model")
factory.remove([(1, l.tag) for l in macro_ll.sides])
factory.remove([(1, s.tag) for s in macro_ll.sides])
factory.synchronize()
factory.removeAllDuplicates()
factory.synchronize()
......@@ -164,4 +164,3 @@ geo.PhysicalGroup.set_group_mesh(1)
gmsh.model.mesh.generate(2)
gmsh.write(str(mesh_file))
run(f"gmsh {str(mesh_file)} &", shell=True, check=True)
gmsh.fltk.run()
ho_homog/geometry/__init__.py
View file @ d350fef5
......@@ -79,7 +79,7 @@ from .transformations import plane_reflection, point_reflection, rotation, trans
DEFAULT_GMSH_OPTIONS = {
"General.Terminal": 1,
"General.Verbosity": 5,
"General.Verbosity": 3,
"Geometry.AutoCoherence": 0,
"Mesh.ColorCarousel": 2,
# * 0=by element type, 1=by elementary entity, 2=by physical entity,3=by partition
......
ho_homog/homog2d.py
View file @ d350fef5
......@@ -61,7 +61,7 @@ class Fenics2DHomogenization(object):
"""
self.rve = fenics_2d_rve
self.topo_dim = topo_dim = fenics_2d_rve.mesh_dim
self.topo_dim = topo_dim = fenics_2d_rve.dim
try:
bottom_left_corner = fenics_2d_rve.bottom_left_corner
except AttributeError:
......
ho_homog/homog_beam.py
View file @ d350fef5
......@@ -58,7 +58,7 @@ class FenicsHomogenizationBeam(object):
"""
self.rve = fenics_2d_rve
self.topo_dim = topo_dim = fenics_2d_rve.mesh_dim
self.topo_dim = topo_dim = fenics_2d_rve.dim
self.pbc = periodicity.PeriodicDomain.pbc_dual_base(
fenics_2d_rve.gen_vect, "X", topo_dim
) # * Périodicité selon X seulement
......
ho_homog/mesh_generate/__init__.py
View file @ d350fef5
......@@ -451,7 +451,6 @@ def Gmsh2DPartFromRVE(cell: Gmsh2DRVE, nb_cells, part_name=None):
tag = cell.phy_surf[i].tag + 1000
name = cell.phy_surf[i].name
phy_surfaces.append(geo.PhysicalGroup(all_surfaces, 2, name, tag))
# gmsh.fltk.run()
for gp in cell.phy_surf:
gp.remove_gmsh()
factory.synchronize()
......
ho_homog/part.py
View file @ d350fef5
......@@ -6,17 +6,16 @@ Created on 17/11/2018
"""
import logging
from pathlib import Path
import dolfin as fe
import matplotlib.pyplot as plt
import numpy as np
import ho_homog.toolbox_FEniCS as fetools
import dolfin as fe
import ho_homog.materials as mat
from subprocess import run
from pathlib import Path
import ho_homog.toolbox_FEniCS as fetools
from ho_homog.toolbox_gmsh import msh_conversion
plt.ioff()
logger = logging.getLogger(__name__) # http://sametmax.com/ecrire-des-logs-en-python/
......@@ -50,22 +49,24 @@ class FenicsPart:
self.elasticity_tensor = fe.as_matrix(materials.get_C())
elif isinstance(materials, dict):
self.elasticity_tensor = mat.mat_per_subdomains(
self.subdomains, self.materials, self.dim
subdomains, materials, self.dim
)
else:
raise TypeError(
"materials parameter must be an instance of Material or a dictionnary that contains Material instances."
)
raise TypeError("Expected a (dict of) Material as materials argument")
if "bottom_left_corner" in kwargs.keys():
self.bottom_left_corner = np.asarray(kwargs["bottom_left_corner"])
self.name = kwargs.get("name", None)
self._mat_area = None
@property
def mat_area(self):
try:
return self._mat_area
except AttributeError:
# Pattern : "lazy property" https://stevenloria.com/lazy-properties/
# http://sametmax.com/creer-un-decorateur-a-la-volee/
if self._mat_area is None:
self._mat_area = fe.assemble(fe.Constant(1) * fe.dx(self.mesh))
return self._mat_area
return self._mat_area
@property
def global_area(self):
if self.global_dimensions is not None:
return np.linalg.det(self.global_dimensions)
......@@ -73,14 +74,38 @@ class FenicsPart:
msg = f"global_size information is lacking for FenicsPart {self}."
raise AttributeError(msg)
def plot_mesh(self):
""" Plot the mesh of a FenicsPart in a matplotlib Figure."""
plt.figure()
try:
fe.plot(self.mesh, title=f"Mesh of {self.name}")
except AttributeError:
fe.plot(self.mesh)
def plot_subdomains(self):
plt.figure()
try:
title = f"Subdomains imported for {self.name}"
except AttributeError:
title = "Imported subdomains"
subdo_plt = fe.plot(self.subdomains, self.mesh, title=title)
plt.colorbar(subdo_plt)
def plot_facet_regions(self):
plt.figure()
facets_val = fetools.get_mesh_function_value(self.facet_regions)
facets_val_range = max(facets_val[1]) - min(facets_val[1])
cmap = plt.cm.get_cmap("viridis", facets_val_range)
facets_plt = fetools.facet_plot2d(self.facet_regions, self.mesh, cmap=cmap)
plt.colorbar(facets_plt[0])
@staticmethod
def file_2_FenicsPart(
def part_from_file(
mesh_path,
materials,
global_dimensions=None,
subdomains_import=False,
plots=True,
explicit_subdo_val=0,
plots=False,
**kwargs,
):
"""Generate an instance of FenicsPart from a .xml or .msh file that contains the mesh.
......@@ -104,8 +129,7 @@ class FenicsPart:
FenicsPart instance
"""
if not isinstance(mesh_path, Path):
mesh_path = Path(mesh_path)
mesh_path = Path(mesh_path)
name = mesh_path.stem
suffix = mesh_path.suffix
......@@ -125,73 +149,48 @@ class FenicsPart:
if suffix == ".xml":
mesh = fe.Mesh(str(mesh_path))
if subdomains_import:
subdo_path = mesh_path.with_name(f"{name}_physical_region.xml")
facet_path = mesh_path.with_name(f"{name}_facet_region.xml")
if subdo_path.exists():
subdomains = fe.MeshFunction("size_t", mesh, str(subdo_path))
subdo_val = fetools.get_MeshFunction_val(subdomains)
logger.info(f"{subdo_val[0]} physical regions imported.")
logger.info(f"The values of their tags are : {subdo_val[1]}")
else:
logger.info(
f"For mesh file {mesh_path.name}, _physical_region.xml file is missing."
)
if facet_path.exists():
facets = fe.MeshFunction("size_t", mesh, str(facet_path))
facets_val = fetools.get_MeshFunction_val(facets)
logger.info(f"{facets_val[0]} facet regions imported.")
logger.info(f"The values of their tags are : {facets_val[1]}")
else:
logger.info(
f"For mesh file {mesh_path.name}, _facet_region.xml file is missing."
)
if suffix == ".xdmf":
subdomains, facets = fetools.import_subdomain_data_xml(mesh, mesh_path)
elif suffix == ".xdmf":
mesh = fetools.xdmf_mesh(mesh_path)
if subdomains_import:
mesh, subdomains, facets = fetools.xdmf_mesh(mesh_path, True)
logger.info(
f"Import of a mesh from {mesh_path} file, with subdomains data"
)
subdomains, facets = fetools.import_subdomain_data_xdmf(mesh, mesh_path)
msg = f"Import of a mesh from {mesh_path} file, with subdomains data"
else:
mesh = fetools.xdmf_mesh(mesh_path, False)
logger.info(
f"Import of a mesh from {mesh_path} file, without subdomains data"
)
msg = f"Import of a mesh from {mesh_path} file, without subdomains data"
logger.info(msg)
else:
raise ValueError(f"expected a mesh path with a suffix `.xml` or `.xdmf`.")
logger.info(f"Import of the mesh : DONE")
if "name" not in kwargs:
kwargs["name"] = name
part = FenicsPart(
mesh, materials, subdomains, global_dimensions, facets, **kwargs
)
if plots:
plt.figure()
fe.plot(mesh, title=f"Mesh of {name}")
if subdomains_import:
plt.figure()
subdo_plt = fe.plot(
subdomains, mesh, title=f"Subdomains imported for {name}"
)
plt.colorbar(subdo_plt)
plt.figure()
facets_val_range = max(facets_val[1]) - min(facets_val[1])
cmap = plt.cm.get_cmap("viridis", facets_val_range)
facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
plt.colorbar(facets_plt[0])
part.plot_mesh()
if subdomains is not None:
part.plot_subdomains()
if facets is not None:
part.plot_facet_regions()
plt.show()
logger.info(f"Import of the mesh : DONE")
return FenicsPart(mesh, materials, subdomains, global_dimensions, facets)
return part
class Fenics2DRVE(FenicsPart):
"""
Contrat : Créer un couple maillage + matériaux pour des RVE 2D, plans, comportant au plus 2 matériaux constitutifs et pouvant contenir plusieurs cellules.
Create a mesh + constituent material(s) pair that represents a 2D RVE.
This class is intended to be used as an input for 2D homogenization schemes, such as Fenics2DHomogenization.
The RVE can be made up of one constitutive material or several constitutive materials.
In the second case, the distribution of materials must be indicated with a dolfin.MeshFunction.
The RVE can contains one or several unit cells.
"""
def __init__(
self,
mesh,
generating_vectors,
material_dict,
subdomains,
facet_regions,
**kwargs,
self, mesh, generating_vectors, materials, subdomains, facet_regions, **kwargs,
):
"""
Parameters
......@@ -200,88 +199,44 @@ class Fenics2DRVE(FenicsPart):
Indicates the subdomains that have been defined in the mesh.
#! L'ordre des facet function à probablement de l'importance pour la suite des opérations
"""
self.mesh = mesh
self.gen_vect = generating_vectors
self.rve_area = np.linalg.det(self.gen_vect)
self.mat_area = fe.assemble(fe.Constant(1) * fe.dx(mesh))
self.mesh_dim = mesh.topology().dim() # dimension d'espace de depart
self.materials = material_dict
self.subdomains = subdomains
self.facet_regions = facet_regions
if isinstance(self.materials, mat.Material):
self.C_per = fe.as_matrix(self.materials.get_C())
else:
self.C_per = mat.mat_per_subdomains(
self.subdomains, self.materials, self.mesh_dim
)
if "bottom_left_corner" in kwargs.keys():
self.bottom_left_corner = np.asarray(kwargs["bottom_left_corner"])
super().__init__(
mesh, materials, subdomains, generating_vectors, facet_regions, **kwargs
)
self.gen_vect = generating_vectors
self.rve_area = self.global_area
self.C_per = self.elasticity_tensor
def epsilon(self, u):
return mat.epsilon(u)
self.epsilon = mat.epsilon
def sigma(self, eps):
return mat.sigma(self.C_per, eps)
def StrainCrossEnergy(self, sig, eps):
def strain_cross_energy(self, sig, eps):
return mat.strain_cross_energy(sig, eps, self.mesh, self.rve_area)
@staticmethod
def gmsh_2_Fenics_2DRVE(gmsh_2D_RVE, material_dict, plots=True):
def rve_from_gmsh2drve(gmsh_2d_rve, materials, plots=True):
"""
Generate an instance of Fenics2DRVE from a instance of the Gmsh2DRVE class.
"""
xml_path = gmsh_2D_RVE.mesh_abs_path.with_suffix(".xml")
cmd = (
"dolfin-convert "
+ gmsh_2D_RVE.mesh_abs_path.as_posix()
+ " "
+ xml_path.as_posix()
)
run(cmd, shell=True, check=True)
mesh = fe.Mesh(xml_path.as_posix())
name = xml_path.stem
subdomains = fe.MeshFunction(
"size_t", mesh, xml_path.with_name(name + "_physical_region.xml").as_posix()
)
facets = fe.MeshFunction(
"size_t", mesh, xml_path.with_name(name + "_facet_region.xml").as_posix()
)
subdo_val = fetools.get_MeshFunction_val(subdomains)
facets_val = fetools.get_MeshFunction_val(facets)
logger.info(
f"{subdo_val[0]} physical regions imported. The values of their tags are : {subdo_val[1]}"
)
logger.info(
f"{facets_val[0]} facet regions imported. The values of their tags are : {facets_val[1]}"
)
if plots:
plt.figure()
subdo_plt = fe.plot(subdomains)
plt.colorbar(subdo_plt)
plt.figure()
cmap = plt.cm.get_cmap("viridis", max(facets_val[1]) - min(facets_val[1]))
facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
plt.colorbar(facets_plt[0])
# clrbar.set_ticks(facets_val)
plt.draw()
logger.info(f"Import of the mesh : DONE")
generating_vectors = gmsh_2D_RVE.gen_vect
msh_conversion(gmsh_2d_rve.mesh_abs_path, format_=".xml")
mesh_path = gmsh_2d_rve.mesh_abs_path.with_suffix(".xml")
gen_vect = gmsh_2d_rve.gen_vect
kwargs = dict()
try:
kwargs["bottom_left_corner"] = gmsh_2D_RVE.bottom_left_corner
kwargs["bottom_left_corner"] = gmsh_2d_rve.bottom_left_corner
except AttributeError:
pass
return Fenics2DRVE(
mesh, generating_vectors, material_dict, subdomains, facets, **kwargs
fenics_rve = Fenics2DRVE.rve_from_file(
mesh_path, gen_vect, materials, plots, **kwargs
)
return fenics_rve
@staticmethod
def file_2_Fenics_2DRVE(
mesh_path, generating_vectors, material_dict, plots=True, **kwargs
def rve_from_file(
mesh_path, generating_vectors, materials, plots=True, **kwargs
):
"""Generate an instance of Fenics2DRVE from a .xml, .msh or .xdmf
file that contains the mesh.
......@@ -293,10 +248,11 @@ class Fenics2DRVE(FenicsPart):
(format Dolfin XML, XDMF or MSH version 2)
generating_vectors : 2D-array
dimensions of the RVE
material_dict : dictionnary
materials : dictionnary
[description] #TODO : à compléter
plots : bool, optional
If True (default) the physical regions and the facet regions are plotted at the end of the import.
If True (default) the physical regions and the facet regions are plotted
at the end of the import.
Returns
-------
......@@ -314,105 +270,28 @@ class Fenics2DRVE(FenicsPart):
with fe.XDMFFile(str(mesh_path)) as file_in:
file_in.read(mesh)
else:
cmd = f"dolfin-convert {mesh_path} {mesh_path.with_suffix('.xml')}"
run(cmd, shell=True, check=True)
msh_conversion(mesh_path, format_=".xml")
mesh_path = mesh_path.with_suffix(".xml")
mesh = fe.Mesh(str(mesh_path))
subdo_path = mesh_path.with_name(name + "_physical_region.xml")
facet_path = mesh_path.with_name(name + "_facet_region.xml")
if subdo_path.exists():
subdomains = fe.MeshFunction("size_t", mesh, subdo_path.as_posix())
subdo_val = fetools.get_MeshFunction_val(subdomains)
logger.info(
f"{subdo_val[0]} physical regions imported. Tags : {subdo_val[1]}"
)
else:
logger.info(
f"For mesh file {mesh_path.name}, _physical_region.xml file is missing."
)
subdomains = None
if facet_path.exists():
facets = fe.MeshFunction("size_t", mesh, facet_path.as_posix())
facets_val = fetools.get_MeshFunction_val(facets)
logger.info(
f"{facets_val[0]} facet regions imported. Tags : {facets_val[1]}"
)
else:
logger.info(
f"For mesh file {mesh_path.name}, _facet_region.xml file is missing."
)
facets = None
logger.info("Import of the mesh : DONE")
# TODO : import des subdomains
subdomains, facets = fetools.import_subdomain_data_xml(mesh, mesh_path)
logger.info("Import of the subdomain data and facet-region data: DONE")
if plots:
# TODO : une seule méthode plot, appelée à chaque fois
plt.figure()
fe.plot(mesh)
if subdomains is not None:
plt.figure()
subdo_plt = fe.plot(subdomains)
plt.colorbar(subdo_plt)
if facets is not None:
plt.figure()
cmap = plt.cm.get_cmap(
"viridis", max(facets_val[1]) - min(facets_val[1])
)
facets_plt = fetools.facet_plot2d(facets, mesh, cmap=cmap)
plt.colorbar(facets_plt[0])
plt.draw()
logger.info(f"Import of the mesh : DONE")
if "name" not in kwargs:
kwargs["name"] = name
return Fenics2DRVE(
mesh, generating_vectors, material_dict, subdomains, facets, **kwargs
fenics_rve = Fenics2DRVE(
mesh, generating_vectors, materials, subdomains, facets, **kwargs
)
if plots:
fenics_rve.plot_mesh()
if subdomains is not None:
fenics_rve.plot_subdomains()
if facets is not None:
fenics_rve.plot_facet_regions()
plt.show()
if __name__ == "__main__":
pass
# TODO : Nettoyer
# geo.init_geo_tools()
# a = 1
# b, k = a, a/3
# panto_test = Gmsh2DRVE.pantograph(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)
# panto_test.mesh_generate()
# os.system(f"gmsh {panto_test.name}.msh &")
# a = 1
# b, k = a, a/3
# panto_test_offset = Gmsh2DRVE.pantograph(a, b, k, 0.1, nb_cells=(2,3), offset=(0.25,0.25), soft_mat=False, name='panto_test_offset')
# panto_test_offset.main_mesh_refinement((0.1,0.5),(0.03,0.3),False)
# panto_test_offset.mesh_generate()
# os.system(f"gmsh {panto_test_offset.name}.msh &")
# L, t = 1, 0.05
# a = L-3*t
# aux_sqr_test = Gmsh2DRVE.auxetic_square(a, L, t, nb_cells=(4,3), soft_mat=False, name='aux_square_test')
# os.system(f"gmsh {aux_sqr_test.name}.brep &")
# aux_sqr_test.main_mesh_refinement((0.1,0.3), (0.01,0.05), False)
# aux_sqr_test.mesh_generate()
# os.system(f"gmsh {aux_sqr_test.name}.msh &")
# a = 1
# b = a
# w = a/50
# r = 4*w
# beam_panto_test = Gmsh2DRVE.beam_pantograph(a, b, w, r, nb_cells=(1, 1), offset=(0., 0.), soft_mat=False, name='beam_panto_test')
# os.system(f"gmsh {beam_panto_test.name}.brep &")
# beam_panto_test.main_mesh_refinement((5*w, a/2),(w/5, w), True)
# beam_panto_test.mesh_generate()
# os.system(f"gmsh {beam_panto_test.name}.msh &")
# gmsh.option.setNumber('Mesh.SurfaceFaces',1) #Display faces of surface mesh?
# gmsh.fltk.run()
# msh.set_background_mesh(field)
# gmsh.option.setNumber('Mesh.CharacteristicLengthExtendFromBoundary',0)
# geo.PhysicalGroup.set_group_mesh(True)
# model.mesh.generate(1)
# model.mesh.generate(2)
# gmsh.write(f"{self.name}.msh")
# os.system(f"gmsh {self.name}.msh &")
return fenics_rve
ho_homog/toolbox_FEniCS.py
View file @ d350fef5
......@@ -37,7 +37,7 @@ DOLFIN_LU_METHODS = {
SUPPORTED_MESH_SUFFIX = (".xml", ".xdmf")
def get_MeshFunction_val(msh_fctn):
def get_mesh_function_value(msh_fctn):
"""
Get information about the set of values that a given MeshFunction outputs.
......@@ -56,7 +56,7 @@ def get_MeshFunction_val(msh_fctn):
val = np.unique(msh_fctn.array())
nb = len(val)
return (nb, val)
return nb, val
def facet_plot2d(
......@@ -116,13 +116,21 @@ def function_from_xdmf(function_space, function_name="", xdmf_path=""):
If multiple functions are requested: #TODO : test
the first argument must be a list of tuples (FunctionSpace, function name)
the second argument must be empty.
Ex: function_from_xdmf([(V,"u"),(W,"eps")], xdmf_path=checkpoint_path)
Returns
-------
dolfin.Function, list of dolfin.Function
New function(s) that is(are) identical to the previously saved function(s).
Example
--------
>>> import dolfin as fe
>>> mesh = xdmf_mesh("mesh.xdmf")
>>> V = fe.VectorFunctionSpace(mesh, "CG", 2)
>>> W = fe.VectorFunctionSpace(mesh, "DG", 1)
>>> function_from_xdmf([(V,"u"),(W,"eps")], xdmf_path="checkpoint.xdmf")
"""
if isinstance(function_space, list):
all_f = list()
with fe.XDMFFile(str(xdmf_path)) as f_in:
......@@ -181,8 +189,6 @@ def function_errornorm(u, v, norm_type="L2", enable_diff_fspace=False):
def local_project(v, fspace, solver_method: str = "", **kwargs):
"""
Info : https://comet-fenics.readthedocs.io/en/latest/tips_and_tricks.html#efficient-projection-on-dg-or-quadrature-spaces #noqa
Parameters
----------
v : [type]
......@@ -207,14 +213,20 @@ def local_project(v, fspace, solver_method: str = "", **kwargs):
Returns
-------
Function
Notes
-----
Source for this code:
https://comet-fenics.readthedocs.io/en/latest/tips_and_tricks.html#efficient-projection-on-dg-or-quadrature-spaces
"""
metadata = kwargs.get("metadata", {})
dX = kwargs.get("dx", fe.dx(metadata=metadata))
dx = kwargs.get("dx", fe.dx(metadata=metadata))
dv = fe.TrialFunction(fspace)
v_ = fe.TestFunction(fspace)
a_proj = fe.inner(dv, v_) * dX
b_proj = fe.inner(v, v_) * dX
a_proj = fe.inner(dv, v_) * dx
b_proj = fe.inner(v, v_) * dx
if solver_method == "LU":
solver = fe.LocalSolver(
......@@ -249,87 +261,138 @@ def _wrap_in_list(obj, name, types=type):
lst = [obj]
for obj in lst:
if not isinstance(obj, types):
raise TypeError(
"expected a (list of) %s as '%s' argument" % (str(types), name)
)
raise TypeError(f"expected a (list of) {types} as {name} argument")
return lst
def xdmf_mesh(mesh_file, import_subdomains=False, facet_file="", physical_file=""):
def xdmf_mesh(mesh_file):
"""Create a FeniCS mesh from a mesh file with xdmf format.
Parameters
----------
mesh_file : str or Path
Path of the file mesh (xdmf file)
import_subdomains : bool, optional
True if information about subdomains have to be imported.
The paths of the auxiliary files that contains information about subdomains
can be indicated with facet_file and physical_file.
The paths used by default are :
- "<mesh path>_facet_region.xdmf" and
- "<mesh path>_physical_region.xdmf" (for subdomains)
Returns
-------