From 18edf99ce0de8c41b40226e9ab65378998f05ac3 Mon Sep 17 00:00:00 2001 From: Baptiste Durand Date: Mon, 11 Jan 2021 17:39:45 +0100 Subject: [PATCH] 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. --- ho_homog/part.py | 341 ++++++++++++------------------------- ho_homog/toolbox_FEniCS.py | 179 ++++++++++++------- test/homog2d_test.py | 20 +-- test/part_test.py | 303 +++++++++++++++++--------------- 4 files changed, 409 insertions(+), 434 deletions(-) diff --git a/ho_homog/part.py b/ho_homog/part.py index 645c804..225af5a 100644 --- a/ho_homog/part.py +++ b/ho_homog/part.py @@ -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 diff --git a/ho_homog/toolbox_FEniCS.py b/ho_homog/toolbox_FEniCS.py index d3c617d..d8f9396 100644 --- a/ho_homog/toolbox_FEniCS.py +++ b/ho_homog/toolbox_FEniCS.py @@ -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 : - - "_facet_region.xdmf" and - - "_physical_region.xdmf" (for subdomains) + + Returns + ------- + dolfin.Mesh + """ + if not isinstance(mesh_file, PurePath): + mesh_file = Path(mesh_file) + if not mesh_file.suffix == ".xdmf": + raise TypeError("Wrong suffix for the path to the mesh.") + mesh = fe.Mesh() + with fe.XDMFFile(str(mesh_file)) as f_in: + f_in.read(mesh) + return mesh + + +_warning_message_ = "The {} file is missing for the mesh {}." + + +def import_subdomain_data_xml(mesh, mesh_file): + """ Import data from the _physical_region.xml and _facet_region.xml files + that result from the mesh conversion to .xml with dolfin-convert + + Parameters + ---------- + mesh: dolfin.Mesh + mesh already imported from the `.xml` file (mesh_path) + mesh_file : pathlib.Path + Path to the mesh file (.xml). + The files that contain the data about the subdomains and the facet regions + are assumed to be in the same folder, with the names: + `[mesh_name]_physical_region.xml` and `[mesh_name]_facet_region.xml`. + """ + name = mesh_file.stem + subdo_path = mesh_file.with_name(name + "_physical_region.xml") + facet_path = mesh_file.with_name(name + "_facet_region.xml") + + if subdo_path.exists(): + subdomains = fe.MeshFunction("size_t", mesh, subdo_path.as_posix()) + subdo_val = get_mesh_function_value(subdomains) + logger.info(f"{subdo_val[0]} physical regions imported. Tags : {subdo_val[1]}") + else: + + logger.warning(_warning_message_.format("_physical_region.xml", mesh_file.name)) + + subdomains = None + + if facet_path.exists(): + facets = fe.MeshFunction("size_t", mesh, facet_path.as_posix()) + facets_val = get_mesh_function_value(facets) + logger.info(f"{facets_val[0]} facet regions imported. Tags : {facets_val[1]}") + else: + logger.warning(_warning_message_.format("_facet_region.xml", mesh_file.name)) + + facets = None + return subdomains, facets + + +def import_subdomain_data_xdmf(mesh, mesh_file, facet_file="", physical_file=""): + + """Import information about subdomains from .xdmf files from + obtained with the mesh conversion .msh -> .xdmf with meshio. + + 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 : + - "_facet_region.xdmf" and + - "_physical_region.xdmf" (for subdomains) + + + Parameters + ---------- + mesh: dolfin.Mesh + mesh already imported from the `.xdmf` file (mesh_path) + mesh_file : pathlib.Path + Path to the mesh file (.xdmf). facet_file : str or Path, optional Path to the mesh auxiliary file that contains subdomains data. Defaults to "" i.e. the default path will be used. physical_file : str or Path, optional Path to the mesh auxiliary file that contains facet regions data. Defaults to "" i.e. the default path will be used. - Returns ------- - Mesh / tuple - If subdomains are not requested : - - The Mesh instance - If subdomains are requested : - - The Mesh instance, - - The MeshFunction for subdomains if it exists else None, - - The MeshFunction for facets if it exists else None, + Tuple (length: 2) + - The MeshFunction for subdomains if it exists else None; + - The MeshFunction for facets if it exists else None. Source ------ Gist meshtagging_mvc.py, June 2018, Michal Habera https://gist.github.com/michalhabera/bbe8a17f788192e53fd758a67cbf3bed - """ - if not isinstance(mesh_file, PurePath): - mesh_file = Path(mesh_file) - if not mesh_file.suffix == ".xdmf": - raise TypeError("Wrong suffix for the path to the mesh.") - mesh = fe.Mesh() - with fe.XDMFFile(str(mesh_file)) as f_in: - f_in.read(mesh) - if not import_subdomains: - return mesh + """ dim = mesh.geometric_dimension() - if not facet_file: - facet_file = mesh_file.with_name(f"{mesh_file.stem}_facet_region.xdmf") - else: - facet_file = Path(facet_file) - if not facet_file.suffix == ".xdmf": - raise TypeError("Wrong suffix for the path to facet regions.") if not physical_file: physical_file = mesh_file.with_name(f"{mesh_file.stem}_physical_region.xdmf") + if not facet_file: + facet_file = mesh_file.with_name(f"{mesh_file.stem}_facet_region.xdmf") + physical_file, facet_file = Path(physical_file), Path(facet_file) + for p, n in zip((physical_file, facet_file), ("subdomains", "facet regions")): + if not p.suffix == ".xdmf": + raise TypeError(f"Wrong suffix for the path to {n}.") + + subdomains, facets = None, None + if physical_file.exists(): + cell_vc = fe.MeshValueCollection("size_t", mesh, dim) + with fe.XDMFFile(str(physical_file)) as f_in: + f_in.read(cell_vc, "cell_data") + subdomains = fe.cpp.mesh.MeshFunctionSizet(mesh, cell_vc) else: - physical_file = Path(physical_file) - if not physical_file.suffix == ".xdmf": - raise TypeError("Wrong suffix for the path to subdomains.") + logger.warning( + _warning_message_.format("_physical_region.xdmf", mesh_file.name) + ) - if not facet_file.exists(): - facet_regions = None - else: + if facet_file.exists(): facet_vc = fe.MeshValueCollection("size_t", mesh, dim - 1) with fe.XDMFFile(str(facet_file)) as f_in: f_in.read(facet_vc, "facet_data") - facet_regions = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc) - if not physical_file.exists(): - subdomains = None + facets = fe.cpp.mesh.MeshFunctionSizet(mesh, facet_vc) else: - cell_vc = fe.MeshValueCollection("size_t", mesh, dim - 1) - with fe.XDMFFile(str(physical_file)) as f_in: - f_in.read(cell_vc, "cell_data") - subdomains = fe.cpp.mesh.MeshFunctionSizet(mesh, cell_vc) - - return mesh, subdomains, facet_regions + logger.warning(_warning_message_.format("_facet_region.xdmf", mesh_file.name)) + return subdomains, facets diff --git a/test/homog2d_test.py b/test/homog2d_test.py index 400676b..f50fa35 100644 --- a/test/homog2d_test.py +++ b/test/homog2d_test.py @@ -36,7 +36,7 @@ def test_homog_EGG_pantograph_1x1(generate_mesh=False): geometry.init_geo_tools() geometry.set_gmsh_option("Mesh.MshFileVersion", 4.1) a = 1 - b, k = a, a/3 + b, k = a, a / 3 r = 0.02 panto_test = mesh_generate.pantograph.pantograph_RVE( a, b, k, r, nb_cells=(1, 1), soft_mat=False, name="panto_rve_1x1" @@ -56,11 +56,10 @@ def test_homog_EGG_pantograph_1x1(generate_mesh=False): meshio.write("panto_rve_1x1.xdmf", geo_only) geometry.reset() - E, nu = 1., 0.3 - material = materials.Material(E, nu, 'cp') - gen_vect = np.array([[4., 0.], [0., 8.]]) - rve = part.Fenics2DRVE.file_2_Fenics_2DRVE( - "panto_rve_1x1.xdmf", gen_vect, material) + E, nu = 1.0, 0.3 + material = materials.Material(E, nu, "cp") + gen_vect = np.array([[4.0, 0.0], [0.0, 8.0]]) + rve = part.Fenics2DRVE.rve_from_file("panto_rve_1x1.xdmf", gen_vect, material) hom_model = homog2d.Fenics2DHomogenization(rve) *localzt_dicts, constit_tensors = hom_model.homogenizationScheme("EGG") @@ -138,11 +137,10 @@ def test_homogeneous_pantograph(generate_mesh=False, save_results=False): meshio.write("homogeneous_panto.xdmf", geo_only) geometry.reset() - E, nu = 1., 0.3 - material = materials.Material(E, nu, 'cp') - gen_vect = np.array([[4., 0.], [0., 8.]]) - rve = part.Fenics2DRVE.file_2_Fenics_2DRVE( - "homogeneous_panto.xdmf", gen_vect, material) + E, nu = 1.0, 0.3 + material = materials.Material(E, nu, "cp") + gen_vect = np.array([[4.0, 0.0], [0.0, 8.0]]) + rve = part.Fenics2DRVE.rve_from_file("homogeneous_panto.xdmf", gen_vect, material) hom_model = homog2d.Fenics2DHomogenization(rve) *localzt_dicts, constit_tensors = hom_model.homogenizationScheme("EGG") Chom_ref = np.array( diff --git a/test/part_test.py b/test/part_test.py index eeaa51e..7b8d3cc 100644 --- a/test/part_test.py +++ b/test/part_test.py @@ -22,21 +22,15 @@ from ho_homog import part model = gmsh.model factory = model.occ -logger = logging.getLogger(__name__) #http://sametmax.com/ecrire-des-logs-en-python/ +logger = logging.getLogger(__name__) # http://sametmax.com/ecrire-des-logs-en-python/ logger.setLevel(logging.INFO) + if __name__ == "__main__": logger_root = logging.getLogger() logger_root.setLevel(logging.INFO) - formatter = logging.Formatter('%(asctime)s :: %(levelname)s :: %(message)s') # Afficher le temps à chaque message - file_handler = RotatingFileHandler(f'activity_{__name__}.log', 'a', 1000000) - file_handler.setLevel(logging.DEBUG) - file_handler.setFormatter(formatter) - logger_root.addHandler(file_handler) #Pour écriture d'un fichier log - formatter = logging.Formatter('%(levelname)s :: %(message)s') - stream_handler = logging.StreamHandler() - stream_handler.setLevel(logging.INFO) - stream_handler.setFormatter(formatter) - logger_root.addHandler(stream_handler) + logging.basicConfig( + format="%(asctime)s :: %(levelname)s :: %(message)s", level=logging.INFO + ) geo.init_geo_tools() geo.set_gmsh_option("General.Verbosity", 3) @@ -50,20 +44,28 @@ def test_rve_from_gmsh2drve(): ) panto_test.main_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) + 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 + 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") -#? Problème de maillage lorsqu'un domaine mou existe + +# test_rve_from_gmsh2drve() +# * Test ok + +# ? Problème de maillage lorsqu'un domaine mou existe # 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') @@ -85,10 +87,10 @@ def test_rve_from_gmsh2drve(): # 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) -#? Fin de la résolution +# ? Fin de la résolution -#? Test d'un mesh avec un domaine mou, gmsh2DRVE puis import dans FEniCS +# ? Test d'un mesh avec un domaine mou, gmsh2DRVE puis import dans FEniCS # 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') @@ -123,177 +125,206 @@ def test_rve_from_gmsh2drve(): # material_dict[subdo.tag] = mat.Material(coeff[0], coeff[1], 'cp') # rve = prt.Fenics2DRVE.gmsh_2_Fenics_2DRVE(aux_test, material_dict) + def test_global_area_2D(): """FenicsPart method global_aera, for a 2D part""" - L_x, L_y = 10., 4. + L_x, L_y = 10.0, 4.0 mesh = fe.RectangleMesh(fe.Point(0.0, 0.0), fe.Point(L_x, L_y), 10, 10) - material = {'0':mat.Material(1., 0.3, 'cp')} - dimensions = np.array(((L_x, 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.,rel=1e-10) + 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) + def test_mat_area(): """FenicsPart method global_aera, for 2D parts created from a gmsh mesh and from a FEniCS mesh""" - L_x, L_y = 4., 5. - H = 1. + L_x, L_y = 4.0, 5.0 + H = 1.0 size = 0.5 - rectangle = mshr.Rectangle(fe.Point(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)) + 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), + ) domain = rectangle - hole domain.set_subdomain(1, rectangle) mesh = mshr.generate_mesh(domain, size) - dimensions = np.array(((L_x, 0.),(0., L_y))) - material = {'0':mat.Material(1., 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) + 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) name = "test_mat_area" local_dir = Path(__file__).parent - mesh_file = local_dir.joinpath(name+'.msh') + mesh_file = local_dir.joinpath(name + ".msh") gmsh.model.add(name) - vertices = [(0., 0.), (0., L_y), (L_x, L_y), (L_x, 0.)] + vertices = [(0.0, 0.0), (0.0, L_y), (L_x, L_y), (L_x, 0.0)] contour = geo.LineLoop([geo.Point(np.array(c)) for c in vertices], False) surface = geo.PlaneSurface(contour) cut_vertices = list() - for local_coord in [(H, 0., 0.), (0., H, 0.), (-H, 0., 0.), (0., -H, 0.)]: - vertex = geo.translation(contour.vertices[2], - np.array(local_coord)) + 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)) cut_vertices.append(vertex) - cut_surface = geo.PlaneSurface(geo.LineLoop(cut_vertices,False)) + cut_surface = geo.PlaneSurface(geo.LineLoop(cut_vertices, False)) for s in [surface, cut_surface]: s.add_gmsh() factory.synchronize() - surface, = geo.surface_bool_cut(surface, cut_surface) + (surface,) = geo.surface_bool_cut(surface, cut_surface) factory.synchronize() for dim_tag in model.getEntities(2): if not dim_tag[1] == surface.tag: model.removeEntities(dim_tag, True) - charact_field = ho_homog.mesh_tools.MathEvalField("0.1") - ho_homog.mesh_tools.set_background_mesh(charact_field) - geo.set_gmsh_option('Mesh.SaveAll', 1) + charact_field = mesh_tools.MathEvalField("0.1") + mesh_tools.set_background_mesh(charact_field) + geo.set_gmsh_option("Mesh.SaveAll", 1) 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) - mesh = fe.Mesh(str(mesh_file.with_suffix('.xml'))) - dimensions = np.array(((L_x, 0.),(0., L_y))) - material = {'0':mat.Material(1., 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) + 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) geo.reset() + 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) - dimensions = np.array(((L_x, 0.),(0., L_y))) + dimensions = np.array(((L_x, 0.0), (0.0, L_y))) E, nu = 1, 0.3 - material = mat.Material(E, nu, 'cp') - rect_part = part.FenicsPart(mesh, materials=material, subdomains=None, - global_dimensions=dimensions, facet_regions=None) - elem_type = 'CG' + material = mat.Material(E, nu, "cp") + rect_part = part.FenicsPart( + mesh, + materials=material, + subdomains=None, + global_dimensions=dimensions, + facet_regions=None, + ) + elem_type = "CG" degree = 2 - 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 = rect_part.sigma(strain) + 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) energy = fe.assemble(fe.inner(stress, strain) * fe.dx(rect_part.mesh)) - energy_theo = E/(1 + nu) * (1 + 28/(15*(1-nu))) + energy_theo = E / (1 + nu) * (1 + 28 / (15 * (1 - nu))) assert energy == approx(energy_theo, rel=1e-13) + def test_2_materials(): """FenicsPart instance initialized with only one instance of Material in the materials dictionnary""" L_x, L_y = 1, 1 - size = 0.1 mesh = fe.RectangleMesh(fe.Point(-L_x, -L_y), fe.Point(L_x, L_y), 20, 20) - dimensions = np.array(((2*L_x, 0.),(0., 2*L_y))) + dimensions = np.array(((2 * L_x, 0.0), (0.0, 2 * L_y))) subdomains = fe.MeshFunction("size_t", mesh, 2) + class Right_part(fe.SubDomain): def inside(self, x, on_boundary): return x[0] >= 0 - fe.DOLFIN_EPS + subdomain_right = Right_part() subdomains.set_all(0) subdomain_right.mark(subdomains, 1) E_1, E_2, nu = 1, 3, 0.3 - materials = { - 0: mat.Material(1, 0.3, 'cp'), - 1: mat.Material(3, 0.3, 'cp') - } + materials = {0: mat.Material(1, 0.3, "cp"), 1: mat.Material(3, 0.3, "cp")} rect_part = part.FenicsPart(mesh, materials, subdomains, dimensions) - elem_type = 'CG' + elem_type = "CG" degree = 2 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 = rect_part.sigma(strain) + 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) energy = fe.assemble(fe.inner(stress, strain) * fe.dx(rect_part.mesh)) - energy_theo = 2*((E_1+E_2)/(1 + nu) * (1 + 28/(15*(1-nu)))) + energy_theo = 2 * ((E_1 + E_2) / (1 + nu) * (1 + 28 / (15 * (1 - nu)))) assert energy == approx(energy_theo, rel=1e-13) + 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 - mesh_file = local_dir.joinpath(name+'.msh') + mesh_file = local_dir.joinpath(name + ".msh") gmsh.model.add(name) - L_x, L_y = 2., 2. - H = 1. - vertices = [(0., 0.), (0., L_y), (L_x, L_y), (L_x, 0.)] + 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)] contour = geo.LineLoop([geo.Point(np.array(c)) for c in vertices], False) surface = geo.PlaneSurface(contour) inclusion_vertices = list() - for coord in [(H/2, -H/2, 0.), (H/2, H/2, 0.), (-H/2, H/2, 0.), (-H/2, -H/2, 0.)]: - vertex = geo.translation(geo.Point((L_x/2, L_y/2)), coord) + 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) inclusion_vertices.append(vertex) inclusion = geo.PlaneSurface(geo.LineLoop(inclusion_vertices, False)) for s in [surface, inclusion]: s.add_gmsh() factory.synchronize() - stiff_s, = geo.surface_bool_cut(surface, inclusion) + (stiff_s,) = geo.surface_bool_cut(surface, inclusion) factory.synchronize() - soft_s, = geo.surface_bool_cut(surface, stiff_s) + (soft_s,) = geo.surface_bool_cut(surface, stiff_s) factory.synchronize() domains = { - 'stiff': geo.PhysicalGroup(stiff_s, 2), - 'soft': geo.PhysicalGroup(soft_s, 2) + "stiff": geo.PhysicalGroup(stiff_s, 2), + "soft": geo.PhysicalGroup(soft_s, 2), } boundaries = { - '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), + "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), } for group in domains.values(): group.add_gmsh() for group in boundaries.values(): group.add_gmsh() - charact_field = ho_homog.mesh_tools.MathEvalField("0.05") - ho_homog.mesh_tools.set_background_mesh(charact_field) - geo.set_gmsh_option('Mesh.SaveAll', 0) + charact_field = mesh_tools.MathEvalField("0.05") + mesh_tools.set_background_mesh(charact_field) + geo.set_gmsh_option("Mesh.SaveAll", 0) 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 = { - domains['soft'].tag: mat.Material(E_1, nu, 'cp'), - domains['stiff'].tag: mat.Material(E_1, nu, 'cp') + domains["soft"].tag: mat.Material(E_1, nu, "cp"), + domains["stiff"].tag: mat.Material(E_1, nu, "cp"), } - test_part = part.FenicsPart.file_2_FenicsPart( - str(mesh_file), materials, subdomains_import=True) - assert test_part.mat_area() == approx(L_x*L_y) - elem_type = 'CG' + 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" degree = 2 V = fe.VectorFunctionSpace(test_part.mesh, elem_type, degree) W = fe.FunctionSpace( @@ -301,30 +332,28 @@ def test_get_domains_gmsh(plots=False): fe.VectorElement(elem_type, test_part.mesh.ufl_cell(), degree, dim=3), ) boundary_conditions = { - 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) - } + 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), + } bcs = list() for tag, val in boundary_conditions.items(): bcs.append(fe.DirichletBC(V, val, test_part.facet_regions, tag)) - ds = fe.Measure( - 'ds', domain=test_part.mesh, - subdomain_data=test_part.facet_regions - ) + ds = fe.Measure("ds", domain=test_part.mesh, subdomain_data=test_part.facet_regions) v = fe.TestFunctions(V) u = fe.TrialFunctions(V) - F = fe.inner( - test_part.sigma(test_part.epsilon(u)), - test_part.epsilon(v) - ) * fe.dx + F = ( + fe.inner(mat.sigma(test_part.elasticity_tensor, mat.epsilon(u)), mat.epsilon(v)) + * fe.dx + ) a, L = fe.lhs(F), fe.rhs(F) u_sol = fe.Function(V) fe.solve(a == L, u_sol, bcs) - strain = fe.project(test_part.epsilon(u_sol),W) + strain = fe.project(mat.epsilon(u_sol), W) if plots: import matplotlib.pyplot as plt + plt.figure() plot = fe.plot(u_sol) plt.colorbar(plot) @@ -338,15 +367,21 @@ def test_get_domains_gmsh(plots=False): plot = fe.plot(strain[2]) plt.colorbar(plot) plt.show() - 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) + 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) materials = { - domains['soft'].tag: mat.Material(E_1, nu, 'cp'), - domains['stiff'].tag: mat.Material(E_2, nu, 'cp') + domains["soft"].tag: mat.Material(E_1, nu, "cp"), + domains["stiff"].tag: mat.Material(E_2, nu, "cp"), } - test_part = part.FenicsPart.file_2_FenicsPart( - str(mesh_file), materials, subdomains_import=True) + test_part = part.FenicsPart.part_from_file( + mesh_file, materials, subdomains_import=True + ) V = fe.VectorFunctionSpace(test_part.mesh, elem_type, degree) W = fe.FunctionSpace( test_part.mesh, @@ -357,15 +392,15 @@ def test_get_domains_gmsh(plots=False): bcs.append(fe.DirichletBC(V, val, test_part.facet_regions, tag)) v = fe.TestFunctions(V) u = fe.TrialFunctions(V) - F = fe.inner( - test_part.sigma(test_part.epsilon(u)), - test_part.epsilon(v) - ) * fe.dx + F = ( + fe.inner(mat.sigma(test_part.elasticity_tensor, mat.epsilon(u)), mat.epsilon(v)) + * fe.dx + ) a, L = fe.lhs(F), fe.rhs(F) u_sol = fe.Function(V) fe.solve(a == L, u_sol, bcs) - strain = test_part.epsilon(u_sol) - stress = test_part.sigma(strain) + strain = mat.epsilon(u_sol) + stress = mat.sigma(test_part.elasticity_tensor, strain) 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) -- GitLab