Commit a73b0450 authored by Baptiste Durand's avatar Baptiste Durand

context.py file for tests

parent 956778f1
# coding: utf8
"""
Created on 12/07/2019
@author: baptiste
Give the individual tests import contex.
Source : https://docs.python-guide.org/writing/structure/#structure-of-the-repository
"""
import site
from pathlib import Path
cur_dir = Path(__file__).resolve().parent
repository_dir = cur_dir.parent
site.addsitedir(repository_dir)
import ho_homog
# coding: utf8
"""
Created on 12/07/2019
@author: baptiste
Give the individual tests import contex.
Source : https://docs.python-guide.org/writing/structure/#structure-of-the-repository
"""
import site
from pathlib import Path
cur_dir = Path(__file__).resolve().parent
repository_dir = cur_dir.parent
site.addsitedir(repository_dir)
import ho_homog
...@@ -6,97 +6,107 @@ Created on 08/04/2019 ...@@ -6,97 +6,107 @@ Created on 08/04/2019
import logging import logging
import dolfin as fe import dolfin as fe
import ho_homog.full_scale_pb as homog_full from context.ho_homog import full_scale_pb as homog_full
from pathlib import Path from pathlib import Path
from pytest import approx from pytest import approx
fe.set_log_level(30) fe.set_log_level(30)
def test_reconstruction_vectors(): def test_reconstruction_vectors():
logger = logging.getLogger("test_reconstruction") logger = logging.getLogger("test_reconstruction")
nb_x = nb_y = 30 nb_x = nb_y = 30
L_x = 2*fe.pi L_x = 2 * fe.pi
L_y = 2*fe.pi L_y = 2 * fe.pi
mesh = fe.RectangleMesh(fe.Point(0, -L_y/2), fe.Point(L_x, L_y/2), nb_x, nb_y) mesh = fe.RectangleMesh(fe.Point(0, -L_y / 2), fe.Point(L_x, L_y / 2), nb_x, nb_y)
elem_type = 'CG' elem_type = "CG"
degree = 2 degree = 2
displ_fspace = fe.VectorFunctionSpace(mesh, elem_type, degree) displ_fspace = fe.VectorFunctionSpace(mesh, elem_type, degree)
strain_fspace = fe.FunctionSpace(mesh, strain_fspace = fe.FunctionSpace(
fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3), mesh, fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3)
) )
scalar_fspace = fe.FunctionSpace(mesh, scalar_fspace = fe.FunctionSpace(
fe.FiniteElement(elem_type, mesh.ufl_cell(), degree), mesh, fe.FiniteElement(elem_type, mesh.ufl_cell(), degree)
) )
macro_fields = { macro_fields = {
'U': [fe.Expression('0.5*x[0]', degree=1), 0], "U": [fe.Expression("0.5*x[0]", degree=1), 0],
'E': [fe.Expression('0.5', degree=0), 0, 0] "E": [fe.Expression("0.5", degree=0), 0, 0],
} }
localization_tensors = { localization_tensors = {
'U': { "U": {
'u':[ "u": [
fe.interpolate(fe.Constant((1.,0.)), displ_fspace).split(), fe.interpolate(fe.Constant((1.0, 0.0)), displ_fspace).split(),
fe.interpolate(fe.Constant((0.,1.)), displ_fspace).split(), fe.interpolate(fe.Constant((0.0, 1.0)), displ_fspace).split(),
], ]
}, },
'E': { "E": {
'u': [ "u": [
[ [
fe.Expression('cos(x[0])', degree=2), fe.Expression("cos(x[0])", degree=2),
fe.Expression('sin(x[1])', degree=2), fe.Expression("sin(x[1])", degree=2),
], ],
[], [],
[], [],
], ],
'eps': [ "eps": [
[ [
fe.Expression('1-sin(x[0])', degree=2), fe.Expression("1-sin(x[0])", degree=2),
fe.Expression('cos(x[1])', degree=2), fe.Expression("cos(x[1])", degree=2),
fe.Constant(value=0), fe.Constant(value=0),
], ],
[], [],
[], [],
] ],
} },
}
function_spaces = {
'u': displ_fspace,
'eps': strain_fspace,
} }
function_spaces = {"u": displ_fspace, "eps": strain_fspace}
exact_sol = { exact_sol = {
'u': fe.project(fe.Expression( "u": fe.project(
("0.5*x[0] + 0.5*cos(x[0])", "0.5*sin(x[1])"), fe.Expression(("0.5*x[0] + 0.5*cos(x[0])", "0.5*sin(x[1])"), degree=2),
degree=2), displ_fspace), displ_fspace,
'eps': fe.project(fe.Expression( ),
(" 0.5 - 0.5*sin(x[0])", "0.5*cos(x[1])", "0"), "eps": fe.project(
degree=2), strain_fspace) fe.Expression((" 0.5 - 0.5*sin(x[0])", "0.5*cos(x[1])", "0"), degree=2),
strain_fspace,
),
} }
reconstr_sol = homog_full.reconstruction( reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces, trunc_order=1) localization_tensors, macro_fields, function_spaces, trunc_order=1
results_file_path = Path(__file__).with_name('test_reconstruction.xdmf') )
results_file_path = Path(__file__).with_name("test_reconstruction.xdmf")
with fe.XDMFFile(str(results_file_path)) as results_file: with fe.XDMFFile(str(results_file_path)) as results_file:
data = [ data = [
(exact_sol['u'], 'disp_exact', 'exact displacement'), (exact_sol["u"], "disp_exact", "exact displacement"),
(exact_sol['eps'], 'strain_exact', 'exact strain'), (exact_sol["eps"], "strain_exact", "exact strain"),
(reconstr_sol['u'], 'disp_reconstruction', 'displacement reconstruction'), (reconstr_sol["u"], "disp_reconstruction", "displacement reconstruction"),
(reconstr_sol['eps'], 'strain_reconstruction', 'strain reconstruction'), (reconstr_sol["eps"], "strain_reconstruction", "strain reconstruction"),
(fe.project(macro_fields['U'][0], scalar_fspace), 'disp_macro', 'disp_macro'), (
(fe.project(macro_fields['E'][0], scalar_fspace), 'strain_macro', 'strain_macro'), fe.project(macro_fields["U"][0], scalar_fspace),
"disp_macro",
"disp_macro",
),
(
fe.project(macro_fields["E"][0], scalar_fspace),
"strain_macro",
"strain_macro",
),
] ]
results_file.parameters["flush_output"] = False results_file.parameters["flush_output"] = False
results_file.parameters["functions_share_mesh"] = True results_file.parameters["functions_share_mesh"] = True
for field, name, descrpt in data: for field, name, descrpt in data:
field.rename(name, descrpt) field.rename(name, descrpt)
results_file.write(field, 0.) results_file.write(field, 0.0)
for f_name in reconstr_sol.keys(): for f_name in reconstr_sol.keys():
dim = exact_sol[f_name].ufl_shape[0] dim = exact_sol[f_name].ufl_shape[0]
exact_norm = fe.norm(exact_sol[f_name], 'L2', mesh) exact_norm = fe.norm(exact_sol[f_name], "L2", mesh)
difference = fe.errornorm(exact_sol[f_name], reconstr_sol[f_name], difference = fe.errornorm(
'L2', mesh=mesh) exact_sol[f_name], reconstr_sol[f_name], "L2", mesh=mesh
error = difference/exact_norm )
logger.debug(f'Relative error for {f_name} = {error}') error = difference / exact_norm
logger.debug(f"Relative error for {f_name} = {error}")
# * ref: Relative errors for u = 3.504e-15; for eps = 2.844e-16 # * ref: Relative errors for u = 3.504e-15; for eps = 2.844e-16
assert error == approx(0., abs=1e-13) assert error == approx(0.0, abs=1e-13)
def test_reconstruction_with_constraint(): def test_reconstruction_with_constraint():
...@@ -105,84 +115,80 @@ def test_reconstruction_with_constraint(): ...@@ -105,84 +115,80 @@ def test_reconstruction_with_constraint():
L_x = 2 L_x = 2
L_y = 2 L_y = 2
mesh = fe.RectangleMesh( mesh = fe.RectangleMesh(
fe.Point(-L_x/2, -L_y/2), fe.Point(L_x/2, L_y/2), fe.Point(-L_x / 2, -L_y / 2), fe.Point(L_x / 2, L_y / 2), nb_x, nb_y
nb_x, nb_y )
)
class PeriodicDomain(fe.SubDomain): class PeriodicDomain(fe.SubDomain):
#? Source : https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html # ? Source : https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html
def __init__(self, tolerance=fe.DOLFIN_EPS): def __init__(self, tolerance=fe.DOLFIN_EPS):
""" vertices stores the coordinates of the 4 unit cell corners""" """ vertices stores the coordinates of the 4 unit cell corners"""
fe.SubDomain.__init__(self, tolerance) fe.SubDomain.__init__(self, tolerance)
self.tol = tolerance self.tol = tolerance
def inside(self, x, on_boundary): def inside(self, x, on_boundary):
Bottom = fe.near(x[1], -L_y/2, self.tol) Bottom = fe.near(x[1], -L_y / 2, self.tol)
return Bottom and on_boundary return Bottom and on_boundary
def map(self, x, y): def map(self, x, y):
SlaveT = fe.near(x[1], L_y/2., self.tol) SlaveT = fe.near(x[1], L_y / 2.0, self.tol)
if SlaveT: if SlaveT:
for i in range(len(x[:])): for i in range(len(x[:])):
y[i] = (x[i] - L_y) if i==1 else x[i] y[i] = (x[i] - L_y) if i == 1 else x[i]
else: else:
for i in range(len(x[:])): for i in range(len(x[:])):
y[i] = 1000. * (L_x + L_y) y[i] = 1000.0 * (L_x + L_y)
pbc = PeriodicDomain() pbc = PeriodicDomain()
elem_type = 'CG' elem_type = "CG"
degree = 2 degree = 2
scalar_fspace = fe.FunctionSpace( scalar_fspace = fe.FunctionSpace(
mesh, mesh,
fe.FiniteElement(elem_type, mesh.ufl_cell(), degree), fe.FiniteElement(elem_type, mesh.ufl_cell(), degree),
constrained_domain=pbc) constrained_domain=pbc,
)
macro_fields = { macro_fields = {
'U': [fe.Expression('0.5*x[0]+x[1]', degree=1)], "U": [fe.Expression("0.5*x[0]+x[1]", degree=1)],
'E': [fe.Expression('0.5', degree=0), fe.Expression('1', degree=0)] "E": [fe.Expression("0.5", degree=0), fe.Expression("1", degree=0)],
} }
localization_tensors = { localization_tensors = {
'U': { "U": {"u": [[fe.interpolate(fe.Constant(1.0), scalar_fspace)]]},
'u': [ "E": {
[fe.interpolate(fe.Constant(1.), scalar_fspace)] "u": [[fe.Expression("x[1]", degree=2)], [fe.Expression("x[0]", degree=2)]]
],
}, },
'E': {
'u': [
[fe.Expression('x[1]', degree=2)],
[fe.Expression('x[0]', degree=2)],
],
}
}
function_spaces = {
'u': scalar_fspace,
} }
localization_rules = {'u': [('U','U'), ('E','E')]} function_spaces = {"u": scalar_fspace}
localization_rules = {"u": [("U", "U"), ("E", "E")]}
exact_sol = fe.project(fe.Expression( exact_sol = fe.project(
"0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", fe.Expression("0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", degree=2), scalar_fspace
degree=2), scalar_fspace) )
reconstr_sol = homog_full.reconstruction( reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces, localization_tensors,
localization_rules=localization_rules, output_request=('u',)) macro_fields,
reconstr_sol = reconstr_sol['u'] function_spaces,
results_file_path = Path(__file__).with_name('test_reconstruction_constraint.xdmf') localization_rules=localization_rules,
output_request=("u",),
)
reconstr_sol = reconstr_sol["u"]
results_file_path = Path(__file__).with_name("test_reconstruction_constraint.xdmf")
with fe.XDMFFile(str(results_file_path)) as results_file: with fe.XDMFFile(str(results_file_path)) as results_file:
data = [ data = [
(exact_sol, 'sol_exact', 'exact solution field'), (exact_sol, "sol_exact", "exact solution field"),
(reconstr_sol, 'sol_reconstructed', 'reconstructed solution field'), (reconstr_sol, "sol_reconstructed", "reconstructed solution field"),
] ]
results_file.parameters["flush_output"] = False results_file.parameters["flush_output"] = False
results_file.parameters["functions_share_mesh"] = True results_file.parameters["functions_share_mesh"] = True
for field, name, descrpt in data: for field, name, descrpt in data:
field.rename(name, descrpt) field.rename(name, descrpt)
results_file.write(field, 0.) results_file.write(field, 0.0)
exact_norm = fe.norm(exact_sol, 'L2', mesh) exact_norm = fe.norm(exact_sol, "L2", mesh)
difference = fe.errornorm(exact_sol, reconstr_sol, 'L2', mesh=mesh) difference = fe.errornorm(exact_sol, reconstr_sol, "L2", mesh=mesh)
error = difference/exact_norm error = difference / exact_norm
logger.debug(f'Relative error = {error}') # * ref: 3.2625e-15 logger.debug(f"Relative error = {error}") # * ref: 3.2625e-15
assert error == approx(0., abs=1e-14) assert error == approx(0.0, abs=1e-14)
def test_select_solver(): def test_select_solver():
"The solver Mumps is selected." "The solver Mumps is selected."
...@@ -191,9 +197,8 @@ def test_select_solver(): ...@@ -191,9 +197,8 @@ def test_select_solver():
L_x = 2 L_x = 2
L_y = 2 L_y = 2
mesh = fe.RectangleMesh( mesh = fe.RectangleMesh(
fe.Point(-L_x/2, -L_y/2), fe.Point(L_x/2, L_y/2), fe.Point(-L_x / 2, -L_y / 2), fe.Point(L_x / 2, L_y / 2), nb_x, nb_y
nb_x, nb_y )
)
class PeriodicDomain(fe.SubDomain): class PeriodicDomain(fe.SubDomain):
def __init__(self, tolerance=fe.DOLFIN_EPS): def __init__(self, tolerance=fe.DOLFIN_EPS):
...@@ -202,65 +207,64 @@ def test_select_solver(): ...@@ -202,65 +207,64 @@ def test_select_solver():
self.tol = tolerance self.tol = tolerance
def inside(self, x, on_boundary): def inside(self, x, on_boundary):
Bottom = fe.near(x[1], -L_y/2, self.tol) Bottom = fe.near(x[1], -L_y / 2, self.tol)
return Bottom and on_boundary return Bottom and on_boundary
def map(self, x, y): def map(self, x, y):
SlaveT = fe.near(x[1], L_y/2., self.tol) SlaveT = fe.near(x[1], L_y / 2.0, self.tol)
if SlaveT: if SlaveT:
for i in range(len(x[:])): for i in range(len(x[:])):
y[i] = (x[i] - L_y) if i==1 else x[i] y[i] = (x[i] - L_y) if i == 1 else x[i]
else: else:
for i in range(len(x[:])): for i in range(len(x[:])):
y[i] = 1000. * (L_x + L_y) y[i] = 1000.0 * (L_x + L_y)
pbc = PeriodicDomain() pbc = PeriodicDomain()
scalar_fspace = fe.FunctionSpace( scalar_fspace = fe.FunctionSpace(
mesh, mesh, fe.FiniteElement("CG", mesh.ufl_cell(), 2), constrained_domain=pbc
fe.FiniteElement('CG', mesh.ufl_cell(), 2), )
constrained_domain=pbc)
macro_fields = { macro_fields = {
'U': [fe.Expression('0.5*x[0]+x[1]', degree=1)], "U": [fe.Expression("0.5*x[0]+x[1]", degree=1)],
'E': [fe.Expression('0.5', degree=0), fe.Expression('1', degree=0)] "E": [fe.Expression("0.5", degree=0), fe.Expression("1", degree=0)],
} }
localization_tensors = { localization_tensors = {
'U': { "U": {"u": [[fe.interpolate(fe.Constant(1.0), scalar_fspace)]]},
'u': [ "E": {
[fe.interpolate(fe.Constant(1.), scalar_fspace)] "u": [[fe.Expression("x[1]", degree=2)], [fe.Expression("x[0]", degree=2)]]
],
}, },
'E': {
'u': [
[fe.Expression('x[1]', degree=2)],
[fe.Expression('x[0]', degree=2)],
],
}
} }
function_spaces = {'u': scalar_fspace} function_spaces = {"u": scalar_fspace}
localization_rules = {'u': [('U', 'U'), ('E', 'E')]} localization_rules = {"u": [("U", "U"), ("E", "E")]}
exact_sol = fe.project(fe.Expression( exact_sol = fe.project(
"0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", fe.Expression("0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", degree=2), scalar_fspace
degree=2), scalar_fspace) )
reconstr_sol = homog_full.reconstruction( reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces, localization_tensors,
localization_rules=localization_rules, output_request=('u',), macro_fields,
proj_solver="mumps") function_spaces,
reconstr_sol = reconstr_sol['u'] localization_rules=localization_rules,
exact_norm = fe.norm(exact_sol, 'L2', mesh) output_request=("u",),
difference = fe.errornorm(exact_sol, reconstr_sol, 'L2', mesh=mesh) proj_solver="mumps",
error = difference/exact_norm )
logger.debug(f'Relative error = {error}') #* ref: 3.2625e-15 reconstr_sol = reconstr_sol["u"]
assert error == approx(0., abs=1e-14) exact_norm = fe.norm(exact_sol, "L2", mesh)
difference = fe.errornorm(exact_sol, reconstr_sol, "L2", mesh=mesh)
error = difference / exact_norm
logger.debug(f"Relative error = {error}") # * ref: 3.2625e-15
assert error == approx(0.0, abs=1e-14)
if __name__ == "__main__": if __name__ == "__main__":
logger_root = logging.getLogger() logger_root = logging.getLogger()
logger_root.setLevel(logging.DEBUG) logger_root.setLevel(logging.DEBUG)
formatter = logging.Formatter('%(asctime)s :: %(levelname)s :: %(message)s',"%H:%M") formatter = logging.Formatter(
"%(asctime)s :: %(levelname)s :: %(message)s", "%H:%M"
)
stream_handler = logging.StreamHandler() stream_handler = logging.StreamHandler()
stream_handler.setLevel(logging.DEBUG) stream_handler.setLevel(logging.DEBUG)
stream_handler.setFormatter(formatter) stream_handler.setFormatter(formatter)
logger_root.addHandler(stream_handler) logger_root.addHandler(stream_handler)
test_reconstruction_vectors() test_reconstruction_vectors()
test_reconstruction_with_constraint() test_reconstruction_with_constraint()
\ No newline at end of file
...@@ -6,7 +6,7 @@ Created on Mon Oct 15 11:00:18 2018 ...@@ -6,7 +6,7 @@ Created on Mon Oct 15 11:00:18 2018
""" """
import os import os
import ho_homog.geometry as geo from .context import ho_homog
import math import math
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -16,6 +16,8 @@ from logging.handlers import RotatingFileHandler ...@@ -16,6 +16,8 @@ from logging.handlers import RotatingFileHandler
import gmsh import gmsh
geo = ho_homog.geometry
# nice shortcuts # nice shortcuts
model = gmsh.model model = gmsh.model
factory = model.occ factory = model.occ
...@@ -502,7 +504,6 @@ def test_ll_modif(): ...@@ -502,7 +504,6 @@ def test_ll_modif():
os.system("gmsh %s.msh &" % name) os.system("gmsh %s.msh &" % name)
def test_gather_line(): def test_gather_line():
""" """
Test of the macro_line_fragments function. Test of the macro_line_fragments function.
......
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