Commit 7330b4d9 authored by Baptiste Durand's avatar Baptiste Durand

Merge branch 'Improve_readability'

parents fad7449f 5332a022
......@@ -23,11 +23,13 @@ from ho_homog import (
homog2d,
materials,
mesh_generate_2D,
mesh_tools,
part,
pckg_logger,
toolbox_FEniCS,
periodicity,
)
from ho_homog.toolbox_FEniCS import function_errornorm
from ho_homog.toolbox_gmsh import process_gmsh_log
logger = logging.getLogger("demo_full_compare")
logger_root = logging.getLogger()
......@@ -54,40 +56,11 @@ logging.getLogger("UFL").setLevel(logging.DEBUG)
logging.getLogger("FFC").setLevel(logging.DEBUG)
fe.set_log_level(20)
def process_gmsh_log(gmsh_log: list, detect_error=True):
"""Treatment of log messages gathered with gmsh.logger.get()"""
err_msg, warn_msg = list(), list()
for line in gmsh_log:
if "error" in line.lower():
err_msg.append(line)
if "warning" in line.lower():
warn_msg.append(line)
gmsh_logger.info("**********")
gmsh_logger.info(
f"{len(gmsh_log)} logging messages got from Gmsh : "
f"{len(err_msg)} errors, {len(warn_msg)} warnings."
)
if err_msg:
gmsh_logger.error("Gmsh errors details :")
for line in err_msg:
gmsh_logger.error(line)
if warn_msg:
gmsh_logger.warning("Gmsh warnings details :")
for line in warn_msg:
gmsh_logger.warning(line)
gmsh_logger.debug("All gmsh messages :")
gmsh_logger.debug(gmsh_log)
gmsh_logger.info("**********")
if detect_error and err_msg:
raise AssertionError("Gmsh logging messages signal errors.")
# * Step 1 : Modeling the geometry of the RVE
geometry.init_geo_tools()
geometry.set_gmsh_option("Mesh.Algorithm", 6)
geometry.set_gmsh_option("Mesh.MshFileVersion", 4.1)
gmsh.option.setNumber("Geometry.Tolerance", 1e-16)
gmsh.option.setNumber("Geometry.Tolerance", 1e-15)
gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.1)
gmsh.option.setNumber("Mesh.LcIntegrationPrecision", 1e-9)
......@@ -113,7 +86,7 @@ gmsh.logger.stop()
gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements()
gmsh.write(str(rve_geo.mesh_abs_path))
rve_path, *_ = mesh_generate_2D.msh_conversion(rve_geo.mesh_abs_path, ".xdmf")
rve_path, *_ = mesh_tools.msh_conversion(rve_geo.mesh_abs_path, ".xdmf")
# * Step 3 : Build the mesh of the part from the mesh of the RVE
......@@ -121,7 +94,7 @@ gmsh.logger.start()
part_geo = mesh_generate_2D.Gmsh2DPartFromRVE(rve_geo, (75, 1))
process_gmsh_log(gmsh.logger.get())
gmsh.logger.stop()
part_path, *_ = mesh_generate_2D.msh_conversion(part_geo.mesh_abs_path, ".xdmf")
part_path, *_ = mesh_tools.msh_conversion(part_geo.mesh_abs_path, ".xdmf")
# * Step 4 : Defining the material properties
E, nu = 1.0, 0.3
......@@ -245,54 +218,6 @@ macro_fields = {
"EGG": [U_hom[3]] + [0] * 11,
}
per_scalar_fnct_cpp_code = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/function/Function.h>
class PeriodicExpr : public dolfin::Expression
{
public:
PeriodicExpr() : dolfin::Expression() {}
void eval(
Eigen::Ref<Eigen::VectorXd> values,
Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell
) const override
{
Eigen::Vector2d coord_equi;
coord_equi[0] = x[0] -per_x*floor(x[0]/per_x);
coord_equi[1] = x[1] -per_y*floor(x[1]/per_y);
f->eval(values, coord_equi);
//passer par un pointeur ? *f->eval(values, coord_equi);
}
std::shared_ptr<dolfin::Function> f;
double per_x;
double per_y;
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<PeriodicExpr, std::shared_ptr<PeriodicExpr>, dolfin::Expression>
(m, "PeriodicExpr", py::dynamic_attr())
.def(py::init<>())
.def_readwrite("f", &PeriodicExpr::f)
.def_readwrite("per_x", &PeriodicExpr::per_x)
.def_readwrite("per_y", &PeriodicExpr::per_y);
}
"""
class PeriodicExpr(fe.UserExpression):
def value_shape(self):
return ()
key_conversion = {"U": "u", "Epsilon": "eps", "Sigma": "sigma"}
localztn_expr = dict()
for key, scd_dict in hom_model.localization.items():
......@@ -307,16 +232,9 @@ for key, scd_dict in hom_model.localization.items():
for i, field in enumerate(localztn_fields):
new_fields = list() # 1 field per component of U, Sigma and Epsilon
for component in field.split():
per_field = PeriodicExpr(degree=3)
periodic_cppcode = fe.compile_cpp_code(
per_scalar_fnct_cpp_code
).PeriodicExpr()
periodic_cppcode.per_x = rve_geo.gen_vect[0, 0]
periodic_cppcode.per_y = rve_geo.gen_vect[1, 1]
component.set_allow_extrapolation(True)
periodic_cppcode.f = component.cpp_object()
# https://bitbucket.org/fenics-project/dolfin/issues/1041/compiledexpression-cant-be-initialized
per_field._cpp_object = periodic_cppcode
per_field = periodicity.PeriodicExpr(
component, rve_geo.gen_vect, degree=3
)
new_fields.append(per_field)
localztn_expr[key][updated_key2].append(new_fields)
function_spaces = {"u": model.displ_fspace, "eps": model.strain_fspace}
......@@ -356,11 +274,8 @@ errors = dict()
for f_name in reconstr_sol.keys():
dim = exact_sol[f_name].ufl_shape[0]
exact_norm = fe.norm(exact_sol[f_name], "L2")
difference = toolbox_FEniCS.function_errornorm(
reconstr_sol[f_name], exact_sol[f_name], "L2"
)
difference = function_errornorm(reconstr_sol[f_name], exact_sol[f_name], "L2")
error = difference / exact_norm
errors[f_name] = error
print(f_name, error, difference, exact_norm)
logger.info(f"Relative error for {f_name} = {error}")
......@@ -112,7 +112,7 @@ macro_s = geo.PlaneSurface(macro_ll)
logger.info("Start boolean operations on surfaces")
pattern_s = [geo.PlaneSurface(ll) for ll in pattern_ll]
rve_s = geo.AbstractSurface.bool_cut(macro_s, pattern_s)
rve_s = geo.surface_bool_cut(macro_s, pattern_s)
rve_s = rve_s[0]
logger.info("Done boolean operations on surfaces")
rve_s_phy = geo.PhysicalGroup([rve_s], 2, "partition_plein")
......
......@@ -5,11 +5,13 @@ Created on 08/04/2019
"""
import logging
import warnings
from pathlib import Path
import dolfin as fe
import numpy as np
import ho_homog
from pathlib import Path
from fenicstools import interpolate_nonmatching_mesh_any
import ho_homog.materials as mat
from ho_homog.toolbox_FEniCS import (
DOLFIN_KRYLOV_METHODS,
DOLFIN_LU_METHODS,
......@@ -18,8 +20,12 @@ from ho_homog.toolbox_FEniCS import (
logger = logging.getLogger(__name__)
GEO_TOLERANCE = ho_homog.GEO_TOLERANCE
mat = ho_homog.materials
try:
from fenicstools import interpolate_nonmatching_mesh_any
except ImportError:
warnings.warn("Import of fenicstools has failed.", ImportWarning)
logger.warning("Import : fenicstools cannot be imported.")
# * For mechanical fields reconstruction
MACRO_FIELDS_NAMES = ["U", "E", "EG", "EGG"]
......@@ -295,9 +301,8 @@ def reconstruction(
macro_kinematic: dict,
function_spaces: dict,
localization_rules: dict = {},
trunc_order: int = 0,
output_request=("u", "eps"),
proj_solver=None,
**kwargs,
):
"""
One argument among localization_rules and trunc_order must be used.
......@@ -318,31 +323,47 @@ def reconstruction(
None can be used to indicate a component that is equal to 0 or irrelevant.
function_spaces : dictionnary
Function space into which each mechanical field have to be built.
keys : 'u', 'eps' or 'sigma'
values : FEniCS function space
- keys : 'u', 'eps' or 'sigma'
- values : FEniCS function space
localization_rules : dict, optional
Indicates the rules that have to be followed for the construction of the fluctuations.
(the default is {}, which means that the trunc_order argument will be used)
trunc_order : int, optional
Order of truncation for the reconstruction of the displacement field
following the higher order homogenization scheme defined in ???.
The rules that have to be followed for the construction of the fluctuations.
Defaults to {} i.e. the trunc_order parameter will be used.
output_request : tuple of strings, optional
Which fields have to be calculated.
This can contain : 'u', eps' and 'sigma'
(the default is ('u', 'eps'), displacement and strain fields will be reconstructed)
proj_solver : string
impose the use of a desired solver for the projections.
Fields that have to be calculated.
The request must be consistent with the keys of other parameters :
- function_spaces
- localization_rules
outputs can be :
======= ===================
name Description
======= ===================
'u' displacement field
'eps' strain field
'sigma' stress field
======= ===================
Defaults to ('u', 'eps').
Return
------
Dictionnary
Mechanical fields with microscopic fluctuations. Keys are "eps", "sigma" and "u" respectively for the strain, stress and displacement fields.
"""
# ? Changer les inputs ?
# ? Remplacer le dictionnaire de functionspace et output_request par un seul argument : list de tuples qui contiennent nom + function_space ?
# ? Comme ça on peut aussi imaginer reconstruire le même champs dans différents espaces ?
Mechanical fields with microscopic fluctuations.
Keys are "eps", "sigma" and "u" respectively for the strain, stress and displacement fields.
Other Parameters
----------------
**kwargs :
Valid kwargs are
=============== ===== =============================================
Key Type Description
=============== ===== =============================================
interp_fnct str The name of the desired function for the interpolations. Allowed values are : "dolfin.interpolate" and "interpolate_nonmatching_mesh_any"
trunc_order int Order of truncation for the reconstruction of the displacement according to the notations used in ???. Override localization_rules parameter.
=============== ====== ============================================
# ! La construction de champs de localisation périodiques doit être faite en dehors de cette fonction.
"""
# TODO : récupérer topo_dim à partir des tenseurs de localisation, ou mieux, à partir des espaces fonctionnels
# TODO : choisir comment on fixe la len des listes correspondantes aux composantes de u et de epsilon.
......@@ -351,14 +372,9 @@ def reconstruction(
# TODO : translation_microstructure: np.array, optional
# TODO : Vector, 1D array (shape (2,) or (3,)), position origin used for the description of the RVE with regards to the macroscopic origin.
solver_param = {}
if proj_solver:
solver_param = {"solver_type": proj_solver}
# Au choix, utilisation de trunc_order ou localization_rules dans les kargs
if localization_rules:
pass
elif trunc_order:
trunc_order = kwargs.pop("trunc_order", None)
if trunc_order:
localization_rules = {
"u": [
(MACRO_FIELDS_NAMES[i], MACRO_FIELDS_NAMES[i])
......@@ -375,6 +391,17 @@ def reconstruction(
# * 'eps': [('E','E'), ('EG', 'EG')]
# *}
interpolate = kwargs.pop("interp_fnct", None)
if interpolate:
if interpolate == "dolfin.interpolate":
interpolate = fe.interpolate
elif interpolate == "interpolate_nonmatching_mesh_any":
interpolate = interpolate_nonmatching_mesh_any
else:
interpolate = fe.interpolate
else:
interpolate = fe.interpolate
reconstructed_fields = dict()
for mecha_field in output_request:
......@@ -417,8 +444,7 @@ def reconstruction(
macro_kin_funct[key] = list()
for comp in field:
if comp:
# macro_kin_funct[key].append(fe.interpolate(comp, scalar_fspace))
macro_kin_funct[key].append(interpolate_nonmatching_mesh_any(comp, scalar_fspace))
macro_kin_funct[key].append(interpolate(comp, scalar_fspace))
else:
macro_kin_funct[key].append(0)
......@@ -431,8 +457,7 @@ def reconstruction(
if not macro_comp:
continue
for i in range(vector_dim):
# loc_comp = fe.interpolate(loc_tens_comps[i], scalar_fspace)
loc_comp = interpolate_nonmatching_mesh_any(loc_tens_comps[i], scalar_fspace)
loc_comp = interpolate(loc_tens_comps[i], scalar_fspace)
contributions[i].append((macro_comp, loc_comp))
# components = [sum(compnt_contrib) for compnt_contrib in contributions]
......@@ -452,21 +477,8 @@ def reconstruction(
# * Components -> vector field
field = fe.Function(fspace)
components_proj = list()
logger.info(f"Projection of reconstructed {mecha_field}...")
components_proj = components
# for scl_field in components:
# # if element_family in ("Discontinuous Lagrange", "Quadrature"):
# # comp = local_project(scl_field, scalar_fspace)
# # else:
# # comp = fe.project(scl_field, scalar_fspace, **solver_param)
# scl_field = scl_field.cpp_object()
# f = fe.Function(scalar_fspace)
# f.interpolate(scl_field)
# components_proj.append(comp)
# logger.info(f"Projection of reconstructed {mecha_field} DONE")
if len(value_shape) == 0:
components_proj = components_proj[0]
assigner.assign(field, components_proj)
components = components[0]
assigner.assign(field, components)
reconstructed_fields[mecha_field] = field
return reconstructed_fields
# coding: utf8
"""
Created on 09/10/2018
@author: baptiste
Définition de classes d'objets géométriques et de fonctions permettant de créer un modèle géométrique de RVE dans gmsh.
sources :
- https://deptinfo-ensip.univ-poitiers.fr/ENS/doku/doku.php/stu:python:pypoo
"""
import copy
import logging
from logging.handlers import RotatingFileHandler
import math
from pathlib import Path
import warnings
import gmsh
import matplotlib.pyplot as plt
import numpy as np
logger = logging.getLogger(__name__)
bndry_logger = logging.getLogger("bndry")
bndry_logger.setLevel(logging.DEBUG)
bndry_logger.propagate = False
log_path = Path('~/ho_homog_log/gmsh_getBoundary_output.log').expanduser()
if not log_path.parent.exists():
log_path.parent.mkdir(mode=0o777, parents=True)
file_handler = RotatingFileHandler(str(log_path), 'a', 1000000, 1)
file_handler.setLevel(logging.DEBUG)
file_handler.setFormatter(
logging.Formatter('%(asctime)s :: %(levelname)s :: %(message)s'))
bndry_logger.addHandler(file_handler)
bndry_logger.warning("***Output of the getBoundary function of the gmsh API***")
bndry_logger.warning("""
*****
This logger and the associated handler are a hacky workaround to avoid errors
related to gmsh.model.occ.getBoundary(). Therefore they should not be removed
or deactivated.
*****""")
# nice shortcuts
model = gmsh.model
factory = model.occ
warnings.simplefilter("always") # ? Doc: https://docs.python.org/3.6/library/warnings.html
# info about gmsh module
logger.info(f"gmsh module path : {Path(gmsh.__file__).resolve()}")
#TODO : placer un asarray dans la def de __init__ pour Point
#TODO : docstring à faire
def unit_vect(v):
""" Renvoie le vecteur normé. Nécessite un vecteur non nul"""
return v / np.linalg.norm(v)
#TODO : docstring à faire
def angle_between(v1, v2, orient=True):
""" Renvoie l'angle en radian, orienté ou non entre les deux vecteurs.
Valeur comprise entre -pi (strictement) et pi.
"""
v1u = unit_vect(v1)
v2u = unit_vect(v2)
if orient:
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.arctan2.html
# EN 2D SEULEMENT
angle = np.arctan2(v2u[1], v2u[0]) - np.arctan2(v1u[1], v1u[0])
if angle <= -math.pi:
angle += 2 * math.pi
if angle > math.pi:
angle -= 2 * math.pi
return angle
else:
return np.arccos(np.clip(np.dot(v1u, v2u), -1.0, 1.0))
#TODO : docstring à faire
#TODO : Refactoring
def bissect(v1, v2):
"""Renvoie la direction de la bissectrice d'un angle défini par deux vecteurs, sous forme d'un vecteur unitaire.
Seulement en 2D pour l'instant.
"""
uv = list(map(unit_vect, [v1, v2]))
e3 = np.array((0., 0., 1.))
# if all(np.shape(x) == (2,) for x in uv): #* travail directement avec vecteurs de dimension 3.
#TODO : voir si d'autres répercussions dans le code
# if uv[0].shape == (2,) and uv[1].shape == (2,) :
# vu=[np.hstack((v,np.zeros((1,)))) for v in vu]
# Complèter les vecteurs 2D par une troisième coordonnée nulle n'est pas nécessaire.
# Fait automatiquement avec numpy.cross
biss = np.cross(uv[1] - uv[0], e3)
#biss = np.delete(biss, 2, 0) # retour à un vecteur 2D. Maintenant on travaille en dimension 3
biss = unit_vect(biss)
return biss
def dual_base(basis):
"""
Calculates the dual basis associated with a given basis.
Parameters
----------
basis : numpy array_like, square matrix
The components of the basis vectors in an orthogonal coordinate system.
2-D array, Dimensions of the matrix : 2×2 or 3×3
Return
------
dual_b : np.array
The components of the vectors that composed the dual base, in the same orthogonal coordinate system.
"""
return np.linalg.inv(basis).T
def set_gmsh_option(option, val):
"""
Set a gmsh option to the given value and print a log message.
Parameters
----------
option : string
val : string or number
Type of valid val depend of the option.
See the gmsh reference manual for more information.
"""
if isinstance(val, (int, float)):
setter = gmsh.option.setNumber
getter = gmsh.option.getNumber
elif isinstance(val, str):
setter = gmsh.option.setString
getter = gmsh.option.getString
else:
raise TypeError("Wrong type of parameter for a gmsh option.")
preval = getter(option)
setter(option, val)
logger.info(f"Gmsh option {option} set to {val} (previously : {preval}).")
def init_geo_tools():
"""
The Gmsh Python API must be initialized before using any functions.
In addition, some options are set to custom values.
"""
gmsh.initialize() #? Utiliser l'argument sys.argv ? cf script boolean.py
set_gmsh_option("General.Terminal", 1)
set_gmsh_option("General.Verbosity", 5)
set_gmsh_option("Geometry.AutoCoherence", 0)
set_gmsh_option("Mesh.ColorCarousel", 2) #0=by element type, 1=by elementary entity, 2=by physical entity, 3=by partition
gmsh.option.setNumber("Mesh.MeshOnlyVisible", 0) #TODO : Should be in the init file of the mesh_tools module.
set_gmsh_option('Mesh.CharacteristicLengthExtendFromBoundary', 0)
set_gmsh_option("Mesh.SaveAll", 0)
set_gmsh_option("Mesh.Binary", 0)
set_gmsh_option('Mesh.MshFileVersion', 2.2)
set_gmsh_option('Mesh.Algorithm', 5) #* 2D mesh algorithm (1=MeshAdapt, 2=Automatic,...)
logger.info("Gmsh SDK version : %s", gmsh.option.getString("General.Version"))
def reset():
"""Throw out all information about the created geometry and remove all gmsh models."""
PhysicalGroup.all_groups = dict()
gmsh.clear()
class Point(object):
"""Classe définissant un point caractérisée par :
- ses coordonnées en 2D ou 3D ;
- (à rajouter ?) la dimension de l'espace de travail;
- (à rajouter ?) son identifiant unique;
- (à rajouter ?) sa longueur caractéristique pour le maillage
Cette classe possède un attribut de classe qui s'incrémente à chaque fois que l'on crée un Point
"""
# TODO faire docstring
def __init__(self, coord=np.array((0., 0.)), mesh_size=0):
"""Constructeur de notre classe. Coordonnées importéées sous forme d'un np.array"""
coord = np.asarray(coord) #*existing arrays are not copied
dim = coord.shape[0]
self.coord = coord
if dim == 2:
self.coord = np.append(self.coord, [0.])
#? Choix : on utilise toujours des coordonnés en 3D. Les points définis en 2D sont mis dans le plan z=0.
self.tag = None
#* Nouveau ! Pour identifier les points "importants" autour des quels il faut raffiner le maillage
self.fine_msh = False #TODO A définir et voir son utilisation...
self.mesh_size = mesh_size
def __repr__(self):
"""Represent a Point object with the coordinates of this point."""
return f"Pt {self.tag} ({str(self.coord)}) "
def __eq__(self, other):
"""
Opérateur de comparaison == redéfini pour les objets de la classe Point
Renvoie True ssi les coordonnées sont égales, sans prendre en compte la signature de l'objet.
"""
if not isinstance(other, Point):
return False
if np.array_equal(self.coord, other.coord):
return True
elif np.allclose(self.coord, other.coord):
return True
else:
return False
def __ne__(self, other):
return not self.__eq__(other)
#TODO Adapter à la 3D avec matplotlib option 3D ?
#? Passer la couleur la taille et tout le reste en kargs, puis passer les kargs en parametres du plot() ?
def plot(self, color="red", size=5):
plt.plot(self.coord[0], self.coord[1], marker='o', markersize=size, color=color)
def add_gmsh(self):
"""
"""
if self.tag: # That means that the geometrical entity has already been instantiated in the gmsh model.
return self.tag #The tag is returned for information purposes only.
self.tag = factory.addPoint(self.coord[0], self.coord[1], self.coord[2], self.mesh_size)
#* Pas utilisé pour le moment. Associer à la fonction dilate de l'API.
# def homothetie(self, centre, k):
# """homothétie de rapport k du point par rapport à un centre de type Point précisé en argument """
# self.coord = centre.coord + k * (self.coord - centre.coord)
#* Pas utilisé pour le moment
# def distance(self, other):
# """Return the distance between this point and a reference point"""
# return np.linalg.norm(self.coord - other.coord)
#### Fonctions permettant des opérations géométriques de bases sur des objets de type Point ####
def centroSym(pt, center):
warnings.warn("Deprecated. Should use the point_reflection function instead.", DeprecationWarning)
""" Renvoie un nouveau point obtenu par symétrie centrale."""
new_coord = -(pt.coord - center.coord) + center.coord
return Point(new_coord)
def mirrorSym(pt, centre, axe):
"""
Renvoie un nouveau point obtenu par symétrie mirroir. axe de symétrie décrit par un point et un vecteur
SEULEMENT EN 2D POUR LE MOMENT
To do : généraliser 2D et 3D
"""
warnings.warn("Deprecated. Should use the plane_reflection function instead.", DeprecationWarning)
new_coord = (2 * (pt.coord - centre.coord).dot(axe) * axe - (pt.coord - centre.coord) + centre.coord)
return Point(new_coord)