Skip to content

GitLab

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

Expression designed for periodic extensionPeriodicExpr is now part of the... 

Expression designed for periodic extensionPeriodicExpr is now part of the ho_homog package.Tested with : `demo_pantograph_same_mesh.py`+ Refactoring reconstruction fuction+ minor corrections in geometry sub-package
parent 711c0c2e
Hide whitespace changes
Inline Side-by-side
Showing with 172 additions and 144 deletions
demo/demo_pantograph_compare_same_mesh.py
View file @ d4147ef0
......@@ -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}")
ho_homog/full_scale_pb.py
View file @ d4147ef0
......@@ -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,49 @@ 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
=============== ===== =============================================
# ! La construction de champs de localisation périodiques doit être faite en dehors de cette fonction.
proj_solver str The solver for the projections
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.
=============== ====== ============================================
"""
# 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 +374,12 @@ 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}
proj_solver = kwargs.pop("proj_solver", None)
solver_param = {"solver_type": proj_solver} if proj_solver else dict()
# 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 +396,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 +449,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 +462,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]
......
ho_homog/geometry/__init__.py
View file @ d4147ef0
......@@ -117,6 +117,7 @@ __all__ = [
# * points
"Point",
# * curves
"Curve",
"AbstractCurve",
"Arc",
"Line",
......
ho_homog/geometry/tools.py
View file @ d4147ef0
......@@ -12,7 +12,6 @@ from . import np, plt, logger, math
from .point import Point
from .curves import Line, Arc
import copy
E3 = np.array((0.0, 0.0, 1.0))
......@@ -109,6 +108,7 @@ def round_corner(inp_pt, pt_amt, pt_avl, r, junction_raduis=False, plot=False):
"""
from .transformations import translation
# Direction en amont et en aval
v_amt = unit_vect(pt_amt.coord - inp_pt.coord)
v_avl = unit_vect(pt_avl.coord - inp_pt.coord)
......
ho_homog/mesh_tools.py
View file @ d4147ef0
......@@ -30,8 +30,9 @@ logger = logging.getLogger(__name__)
class Field(object):
"""
Class mère (=générique) pour représenter les champs scalaires utilisés pour prescrire la taile caractéristique des éléments dans gmsh.
Info pour l'héritage : https://openclassrooms.com/fr/courses/235344-apprenez-a-programmer-en-python/233164-lheritage
Représentation des champs scalaires utilisés pour prescrire
la taile caractéristique des éléments dans gmsh.
Tutorial : gmsh-4.0.2-Linux64/share/doc/gmsh/demos/api/t10.py
"""
......@@ -41,27 +42,23 @@ class Field(object):
self.parents = parent_fields
def set_params(self):
"""
Cette partie sera précisée pour chaque type de champ de scalaire.
"""
# ? C'est a priori pas nécessaire de définir une méthode set_params dans cette classe.
"""Cette méthode sera précisée pour chaque type de champ de scalaire."""
pass
def add_gmsh(self):
"""
Partie générique des intructions nécessaires pour ajouter un field de contrôle de la taille des éléments.
Fonction générique pour l'ajout d'un characteristic length field au modèle gmsh.
"""
if (
self.tag
): # That means that the field has already been instantiated in the gmsh model.
return self.tag # The tag is returned for information purposes only.
if self.tag:
return None
if self.parents:
for p_field in self.parents:
if not p_field.tag:
p_field.add_gmsh()
self.tag = api_field.add(self.f_type)
self.set_params()
return self.tag
class AttractorField(Field):
......
ho_homog/periodicity.py
View file @ d4147ef0
......@@ -165,3 +165,87 @@ class PeriodicDomain(fe.SubDomain):
per_vectors.append(basis[1])
return PeriodicDomain(per_vectors, master_tests, slave_tests, dim, tol)
# * Periodic localization fields
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);
}
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);
}
"""
periodic_cppcode = fe.compile_cpp_code(per_scalar_fnct_cpp_code)
class PeriodicExpr(fe.UserExpression):
"""Periodic extension of a (periodic) function. Represented by an expression"""
def __init__(self, base_function, per_vectors, topo_dim=2, **kwargs):
"""Create a periodic extension of a function which is defined only on a cell.
For now, only the case of periodic vectors aligned with the global basis
in 2D is supported.
Parameters
----------
base_function : dolfin Function
Function that will be extended by periodicity.
Must be at least defined on a cell delimited by the 4 points :
O, vp1, vp2, vp1+vp2 where vp1 and vp2 are the periodicity vectors.
per_vectors : ndarray
periodicity vectors in column
topo_dim : int, optional
[description], by default 2
"""
self.base_function = base_function
self.per_vectors = np.asarray(per_vectors)
assert self.per_vectors.shape == (2, 2), "Invalid shape. Only 2D supported"
diag_test = np.allclose(self.per_vectors[[0, 1], [1, 0]], [0.0, 0.0], atol=1e-12)
assert diag_test, "Periodicity vectors must be align with global basis."
# TODO : généralisation à des vecteurs de périodicité non alignés avec le repère.
super().__init__(degree=kwargs["degree"])
per_cpp_expr = periodic_cppcode.PeriodicExpr()
per_cpp_expr.per_x = self.per_vectors[0, 0]
per_cpp_expr.per_y = self.per_vectors[1, 1]
base_function.set_allow_extrapolation(True)
per_cpp_expr.f = base_function.cpp_object()
self._cpp_object = per_cpp_expr
def value_shape(self):
return ()
test/test_mesh_generate_2D.py
View file @ d4147ef0
......@@ -5,6 +5,7 @@ Created on 12/06/2019
"""
from ho_homog import geometry as geo
from ho_homog import mesh_tools
from ho_homog import mesh_generate_2D as mesh2D
import gmsh
......@@ -29,5 +30,5 @@ def test_rve_2_part():
gmsh.model.mesh.renumberElements()
gmsh.write(str(panto_rve.mesh_abs_path))
panto_part = mesh2D.Gmsh2DPartFromRVE(panto_rve, (2, 3))
mesh2D.msh_conversion(panto_rve.mesh_abs_path, ".xdmf")
mesh2D.msh_conversion(panto_part.mesh_abs_path, ".xdmf")
mesh_tools.msh_conversion(panto_rve.mesh_abs_path, ".xdmf")
mesh_tools.msh_conversion(panto_part.mesh_abs_path, ".xdmf")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment