Commit 18edf99c authored by Baptiste Durand's avatar Baptiste Durand

Reorganize FenicsPart and Fenics2DRVE

Use of FenicsPart __init__ in Fenics2DRVE;
Define new functions for plotting the mesh and the subdomains and importing information about subdomains;
Split import of a mesh from a .xdmf and import of the subdomain data from the complementary .xdmf file;
Use of msh_conversion systematically when a conversion is needed.
parent ab886331
...@@ -6,17 +6,16 @@ Created on 17/11/2018 ...@@ -6,17 +6,16 @@ Created on 17/11/2018
""" """
import logging import logging
from pathlib import Path
import dolfin as fe
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import ho_homog.toolbox_FEniCS as fetools
import dolfin as fe
import ho_homog.materials as mat import ho_homog.materials as mat
from subprocess import run import ho_homog.toolbox_FEniCS as fetools
from pathlib import Path
from ho_homog.toolbox_gmsh import msh_conversion from ho_homog.toolbox_gmsh import msh_conversion
plt.ioff()
logger = logging.getLogger(__name__) # http://sametmax.com/ecrire-des-logs-en-python/ logger = logging.getLogger(__name__) # http://sametmax.com/ecrire-des-logs-en-python/
...@@ -50,22 +49,24 @@ class FenicsPart: ...@@ -50,22 +49,24 @@ class FenicsPart:
self.elasticity_tensor = fe.as_matrix(materials.get_C()) self.elasticity_tensor = fe.as_matrix(materials.get_C())
elif isinstance(materials, dict): elif isinstance(materials, dict):
self.elasticity_tensor = mat.mat_per_subdomains( self.elasticity_tensor = mat.mat_per_subdomains(
self.subdomains, self.materials, self.dim subdomains, materials, self.dim
) )
else: else:
raise TypeError( raise TypeError("Expected a (dict of) Material as materials argument")
"materials parameter must be an instance of Material or a dictionnary that contains Material instances."
)
if "bottom_left_corner" in kwargs.keys(): if "bottom_left_corner" in kwargs.keys():
self.bottom_left_corner = np.asarray(kwargs["bottom_left_corner"]) 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): def mat_area(self):
try: # Pattern : "lazy property" https://stevenloria.com/lazy-properties/
return self._mat_area # http://sametmax.com/creer-un-decorateur-a-la-volee/
except AttributeError: if self._mat_area is None:
self._mat_area = fe.assemble(fe.Constant(1) * fe.dx(self.mesh)) 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): def global_area(self):
if self.global_dimensions is not None: if self.global_dimensions is not None:
return np.linalg.det(self.global_dimensions) return np.linalg.det(self.global_dimensions)
...@@ -73,14 +74,38 @@ class FenicsPart: ...@@ -73,14 +74,38 @@ class FenicsPart:
msg = f"global_size information is lacking for FenicsPart {self}." msg = f"global_size information is lacking for FenicsPart {self}."
raise AttributeError(msg) 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 @staticmethod
def file_2_FenicsPart( def part_from_file(
mesh_path, mesh_path,
materials, materials,
global_dimensions=None, global_dimensions=None,
subdomains_import=False, subdomains_import=False,
plots=True, plots=False,
explicit_subdo_val=0,
**kwargs, **kwargs,
): ):
"""Generate an instance of FenicsPart from a .xml or .msh file that contains the mesh. """Generate an instance of FenicsPart from a .xml or .msh file that contains the mesh.
...@@ -104,8 +129,7 @@ class FenicsPart: ...@@ -104,8 +129,7 @@ class FenicsPart:
FenicsPart instance FenicsPart instance
""" """
if not isinstance(mesh_path, Path): mesh_path = Path(mesh_path)
mesh_path = Path(mesh_path)
name = mesh_path.stem name = mesh_path.stem
suffix = mesh_path.suffix suffix = mesh_path.suffix
...@@ -125,73 +149,48 @@ class FenicsPart: ...@@ -125,73 +149,48 @@ class FenicsPart:
if suffix == ".xml": if suffix == ".xml":
mesh = fe.Mesh(str(mesh_path)) mesh = fe.Mesh(str(mesh_path))
if subdomains_import: if subdomains_import:
subdo_path = mesh_path.with_name(f"{name}_physical_region.xml") subdomains, facets = fetools.import_subdomain_data_xml(mesh, mesh_path)
facet_path = mesh_path.with_name(f"{name}_facet_region.xml") elif suffix == ".xdmf":
if subdo_path.exists(): mesh = fetools.xdmf_mesh(mesh_path)
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":
if subdomains_import: if subdomains_import:
mesh, subdomains, facets = fetools.xdmf_mesh(mesh_path, True) subdomains, facets = fetools.import_subdomain_data_xdmf(mesh, mesh_path)
logger.info( msg = f"Import of a mesh from {mesh_path} file, with subdomains data"
f"Import of a mesh from {mesh_path} file, with subdomains data"
)
else: else:
mesh = fetools.xdmf_mesh(mesh_path, False) msg = f"Import of a mesh from {mesh_path} file, without subdomains data"
logger.info( logger.info(msg)
f"Import of a mesh from {mesh_path} file, without subdomains data" 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: if plots:
plt.figure() part.plot_mesh()
fe.plot(mesh, title=f"Mesh of {name}") if subdomains is not None:
if subdomains_import: part.plot_subdomains()
plt.figure() if facets is not None:
subdo_plt = fe.plot( part.plot_facet_regions()
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])
plt.show() plt.show()
logger.info(f"Import of the mesh : DONE")
return FenicsPart(mesh, materials, subdomains, global_dimensions, facets) return part
class Fenics2DRVE(FenicsPart): 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__( def __init__(
self, self, mesh, generating_vectors, materials, subdomains, facet_regions, **kwargs,
mesh,
generating_vectors,
material_dict,
subdomains,
facet_regions,
**kwargs,
): ):
""" """
Parameters Parameters
...@@ -200,88 +199,44 @@ class Fenics2DRVE(FenicsPart): ...@@ -200,88 +199,44 @@ class Fenics2DRVE(FenicsPart):
Indicates the subdomains that have been defined in the mesh. 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 #! 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): super().__init__(
self.C_per = fe.as_matrix(self.materials.get_C()) mesh, materials, subdomains, generating_vectors, facet_regions, **kwargs
else: )
self.C_per = mat.mat_per_subdomains( self.gen_vect = generating_vectors
self.subdomains, self.materials, self.mesh_dim self.rve_area = self.global_area
) self.C_per = self.elasticity_tensor
if "bottom_left_corner" in kwargs.keys():
self.bottom_left_corner = np.asarray(kwargs["bottom_left_corner"])
def epsilon(self, u): self.epsilon = mat.epsilon
return mat.epsilon(u)
def sigma(self, eps): def sigma(self, eps):
return mat.sigma(self.C_per, 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) return mat.strain_cross_energy(sig, eps, self.mesh, self.rve_area)
@staticmethod @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. Generate an instance of Fenics2DRVE from a instance of the Gmsh2DRVE class.
""" """
xml_path = gmsh_2D_RVE.mesh_abs_path.with_suffix(".xml") msh_conversion(gmsh_2d_rve.mesh_abs_path, format_=".xml")
cmd = ( mesh_path = gmsh_2d_rve.mesh_abs_path.with_suffix(".xml")
"dolfin-convert " gen_vect = gmsh_2d_rve.gen_vect
+ 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
kwargs = dict() kwargs = dict()
try: try:
kwargs["bottom_left_corner"] = gmsh_2D_RVE.bottom_left_corner kwargs["bottom_left_corner"] = gmsh_2d_rve.bottom_left_corner
except AttributeError: except AttributeError:
pass pass
return Fenics2DRVE( fenics_rve = Fenics2DRVE.rve_from_file(
mesh, generating_vectors, material_dict, subdomains, facets, **kwargs mesh_path, gen_vect, materials, plots, **kwargs
) )
return fenics_rve
@staticmethod @staticmethod
def file_2_Fenics_2DRVE( def rve_from_file(
mesh_path, generating_vectors, material_dict, plots=True, **kwargs mesh_path, generating_vectors, materials, plots=True, **kwargs
): ):
"""Generate an instance of Fenics2DRVE from a .xml, .msh or .xdmf """Generate an instance of Fenics2DRVE from a .xml, .msh or .xdmf
file that contains the mesh. file that contains the mesh.
...@@ -293,10 +248,11 @@ class Fenics2DRVE(FenicsPart): ...@@ -293,10 +248,11 @@ class Fenics2DRVE(FenicsPart):
(format Dolfin XML, XDMF or MSH version 2) (format Dolfin XML, XDMF or MSH version 2)
generating_vectors : 2D-array generating_vectors : 2D-array
dimensions of the RVE dimensions of the RVE
material_dict : dictionnary materials : dictionnary
[description] #TODO : à compléter [description] #TODO : à compléter
plots : bool, optional 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 Returns
------- -------
...@@ -314,105 +270,28 @@ class Fenics2DRVE(FenicsPart): ...@@ -314,105 +270,28 @@ class Fenics2DRVE(FenicsPart):
with fe.XDMFFile(str(mesh_path)) as file_in: with fe.XDMFFile(str(mesh_path)) as file_in:
file_in.read(mesh) file_in.read(mesh)
else: else:
cmd = f"dolfin-convert {mesh_path} {mesh_path.with_suffix('.xml')}" msh_conversion(mesh_path, format_=".xml")
run(cmd, shell=True, check=True)
mesh_path = mesh_path.with_suffix(".xml") mesh_path = mesh_path.with_suffix(".xml")
mesh = fe.Mesh(str(mesh_path)) mesh = fe.Mesh(str(mesh_path))
subdo_path = mesh_path.with_name(name + "_physical_region.xml") logger.info("Import of the mesh : DONE")
facet_path = mesh_path.with_name(name + "_facet_region.xml") # TODO : import des subdomains
if subdo_path.exists(): subdomains, facets = fetools.import_subdomain_data_xml(mesh, mesh_path)
subdomains = fe.MeshFunction("size_t", mesh, subdo_path.as_posix()) logger.info("Import of the subdomain data and facet-region data: DONE")
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
if plots: if "name" not in kwargs:
# TODO : une seule méthode plot, appelée à chaque fois kwargs["name"] = name
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")
return Fenics2DRVE( fenics_rve = Fenics2DRVE(
mesh, generating_vectors, material_dict, subdomains, facets, **kwargs 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__": return fenics_rve
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 &")
...@@ -37,7 +37,7 @@ DOLFIN_LU_METHODS = { ...@@ -37,7 +37,7 @@ DOLFIN_LU_METHODS = {
SUPPORTED_MESH_SUFFIX = (".xml", ".xdmf") 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. Get information about the set of values that a given MeshFunction outputs.
...@@ -56,7 +56,7 @@ def get_MeshFunction_val(msh_fctn): ...@@ -56,7 +56,7 @@ def get_MeshFunction_val(msh_fctn):
val = np.unique(msh_fctn.array()) val = np.unique(msh_fctn.array())
nb = len(val) nb = len(val)
return (nb, val) return nb, val
def facet_plot2d( def facet_plot2d(
...@@ -116,13 +116,21 @@ def function_from_xdmf(function_space, function_name="", xdmf_path=""): ...@@ -116,13 +116,21 @@ def function_from_xdmf(function_space, function_name="", xdmf_path=""):
If multiple functions are requested: #TODO : test If multiple functions are requested: #TODO : test
the first argument must be a list of tuples (FunctionSpace, function name) the first argument must be a list of tuples (FunctionSpace, function name)
the second argument must be empty. the second argument must be empty.
Ex: function_from_xdmf([(V,"u"),(W,"eps")], xdmf_path=checkpoint_path)
Returns Returns
------- -------
dolfin.Function, list of dolfin.Function dolfin.Function, list of dolfin.Function
New function(s) that is(are) identical to the previously saved function(s). 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): if isinstance(function_space, list):
all_f = list() all_f = list()
with fe.XDMFFile(str(xdmf_path)) as f_in: 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): ...@@ -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): 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 Parameters
---------- ----------
v : [type] v : [type]
...@@ -207,14 +213,20 @@ def local_project(v, fspace, solver_method: str = "", **kwargs): ...@@ -207,14 +213,20 @@ def local_project(v, fspace, solver_method: str = "", **kwargs):
Returns Returns
------- -------
Function 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", {}) metadata = kwargs.get("metadata", {})
dX = kwargs.get("dx", fe.dx(metadata=metadata)) dx = kwargs.get("dx", fe.dx(metadata=metadata))
dv = fe.TrialFunction(fspace) dv = fe.TrialFunction(fspace)
v_ = fe.TestFunction(fspace) v_ = fe.TestFunction(fspace)
a_proj = fe.inner(dv, v_) * dX a_proj = fe.inner(dv, v_) * dx
b_proj = fe.inner(v, v_) * dX