Skip to content

GitLab

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

Merge branch 'improve_structure'

parents 46f86d8f a73b0450
Hide whitespace changes
Inline Side-by-side
Showing with 176 additions and 133 deletions
test/__init__.py 0 → 100644
View file @ ea814ebe
# 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
test/context.py 0 → 100644
View file @ ea814ebe
# 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
test/full_scale_pb_test.py
View file @ ea814ebe
......@@ -6,97 +6,107 @@ Created on 08/04/2019
import logging
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 pytest import approx
fe.set_log_level(30)
def test_reconstruction_vectors():
logger = logging.getLogger("test_reconstruction")
nb_x = nb_y = 30
L_x = 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)
elem_type = 'CG'
L_x = 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)
elem_type = "CG"
degree = 2
displ_fspace = fe.VectorFunctionSpace(mesh, elem_type, degree)
strain_fspace = fe.FunctionSpace(mesh,
fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3),
)
scalar_fspace = fe.FunctionSpace(mesh,
fe.FiniteElement(elem_type, mesh.ufl_cell(), degree),
)
strain_fspace = fe.FunctionSpace(
mesh, fe.VectorElement(elem_type, mesh.ufl_cell(), degree, dim=3)
)
scalar_fspace = fe.FunctionSpace(
mesh, fe.FiniteElement(elem_type, mesh.ufl_cell(), degree)
)
macro_fields = {
'U': [fe.Expression('0.5*x[0]', degree=1), 0],
'E': [fe.Expression('0.5', degree=0), 0, 0]
"U": [fe.Expression("0.5*x[0]", degree=1), 0],
"E": [fe.Expression("0.5", degree=0), 0, 0],
}
localization_tensors = {
'U': {
'u':[
fe.interpolate(fe.Constant((1.,0.)), displ_fspace).split(),
fe.interpolate(fe.Constant((0.,1.)), displ_fspace).split(),
],
"U": {
"u": [
fe.interpolate(fe.Constant((1.0, 0.0)), displ_fspace).split(),
fe.interpolate(fe.Constant((0.0, 1.0)), displ_fspace).split(),
]
},
'E': {
'u': [
"E": {
"u": [
[
fe.Expression('cos(x[0])', degree=2),
fe.Expression('sin(x[1])', degree=2),
fe.Expression("cos(x[0])", degree=2),
fe.Expression("sin(x[1])", degree=2),
],
[],
[],
],
'eps': [
"eps": [
[
fe.Expression('1-sin(x[0])', degree=2),
fe.Expression('cos(x[1])', degree=2),
fe.Expression("1-sin(x[0])", degree=2),
fe.Expression("cos(x[1])", degree=2),
fe.Constant(value=0),
],
[],
[],
]
}
}
function_spaces = {
'u': displ_fspace,
'eps': strain_fspace,
],
},
}
function_spaces = {"u": displ_fspace, "eps": strain_fspace}
exact_sol = {
'u': fe.project(fe.Expression(
("0.5*x[0] + 0.5*cos(x[0])", "0.5*sin(x[1])"),
degree=2), displ_fspace),
'eps': fe.project(fe.Expression(
(" 0.5 - 0.5*sin(x[0])", "0.5*cos(x[1])", "0"),
degree=2), strain_fspace)
"u": fe.project(
fe.Expression(("0.5*x[0] + 0.5*cos(x[0])", "0.5*sin(x[1])"), degree=2),
displ_fspace,
),
"eps": fe.project(
fe.Expression((" 0.5 - 0.5*sin(x[0])", "0.5*cos(x[1])", "0"), degree=2),
strain_fspace,
),
}
reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces, trunc_order=1)
results_file_path = Path(__file__).with_name('test_reconstruction.xdmf')
localization_tensors, macro_fields, function_spaces, trunc_order=1
)
results_file_path = Path(__file__).with_name("test_reconstruction.xdmf")
with fe.XDMFFile(str(results_file_path)) as results_file:
data = [
(exact_sol['u'], 'disp_exact', 'exact displacement'),
(exact_sol['eps'], 'strain_exact', 'exact strain'),
(reconstr_sol['u'], 'disp_reconstruction', 'displacement 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'),
(exact_sol["u"], "disp_exact", "exact displacement"),
(exact_sol["eps"], "strain_exact", "exact strain"),
(reconstr_sol["u"], "disp_reconstruction", "displacement 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",
),
]
results_file.parameters["flush_output"] = False
results_file.parameters["functions_share_mesh"] = True
for field, name, descrpt in data:
field.rename(name, descrpt)
results_file.write(field, 0.)
results_file.write(field, 0.0)
for f_name in reconstr_sol.keys():
dim = exact_sol[f_name].ufl_shape[0]
exact_norm = fe.norm(exact_sol[f_name], 'L2', mesh)
difference = fe.errornorm(exact_sol[f_name], reconstr_sol[f_name],
'L2', mesh=mesh)
error = difference/exact_norm
logger.debug(f'Relative error for {f_name} = {error}')
exact_norm = fe.norm(exact_sol[f_name], "L2", mesh)
difference = fe.errornorm(
exact_sol[f_name], reconstr_sol[f_name], "L2", mesh=mesh
)
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
assert error == approx(0., abs=1e-13)
assert error == approx(0.0, abs=1e-13)
def test_reconstruction_with_constraint():
......@@ -105,84 +115,80 @@ def test_reconstruction_with_constraint():
L_x = 2
L_y = 2
mesh = fe.RectangleMesh(
fe.Point(-L_x/2, -L_y/2), fe.Point(L_x/2, L_y/2),
nb_x, nb_y
)
fe.Point(-L_x / 2, -L_y / 2), fe.Point(L_x / 2, L_y / 2), nb_x, nb_y
)
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):
""" vertices stores the coordinates of the 4 unit cell corners"""
fe.SubDomain.__init__(self, tolerance)
self.tol = tolerance
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
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:
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:
for i in range(len(x[:])):
y[i] = 1000. * (L_x + L_y)
y[i] = 1000.0 * (L_x + L_y)
pbc = PeriodicDomain()
elem_type = 'CG'
elem_type = "CG"
degree = 2
scalar_fspace = fe.FunctionSpace(
mesh,
fe.FiniteElement(elem_type, mesh.ufl_cell(), degree),
constrained_domain=pbc)
constrained_domain=pbc,
)
macro_fields = {
'U': [fe.Expression('0.5*x[0]+x[1]', degree=1)],
'E': [fe.Expression('0.5', degree=0), fe.Expression('1', degree=0)]
"U": [fe.Expression("0.5*x[0]+x[1]", degree=1)],
"E": [fe.Expression("0.5", degree=0), fe.Expression("1", degree=0)],
}
localization_tensors = {
'U': {
'u': [
[fe.interpolate(fe.Constant(1.), scalar_fspace)]
],
"U": {"u": [[fe.interpolate(fe.Constant(1.0), scalar_fspace)]]},
"E": {
"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(
"0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]",
degree=2), scalar_fspace)
exact_sol = fe.project(
fe.Expression("0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", degree=2), scalar_fspace
)
reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces,
localization_rules=localization_rules, output_request=('u',))
reconstr_sol = reconstr_sol['u']
results_file_path = Path(__file__).with_name('test_reconstruction_constraint.xdmf')
localization_tensors,
macro_fields,
function_spaces,
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:
data = [
(exact_sol, 'sol_exact', 'exact solution field'),
(reconstr_sol, 'sol_reconstructed', 'reconstructed solution field'),
(exact_sol, "sol_exact", "exact solution field"),
(reconstr_sol, "sol_reconstructed", "reconstructed solution field"),
]
results_file.parameters["flush_output"] = False
results_file.parameters["functions_share_mesh"] = True
for field, name, descrpt in data:
field.rename(name, descrpt)
results_file.write(field, 0.)
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., abs=1e-14)
results_file.write(field, 0.0)
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)
def test_select_solver():
"The solver Mumps is selected."
......@@ -191,9 +197,8 @@ def test_select_solver():
L_x = 2
L_y = 2
mesh = fe.RectangleMesh(
fe.Point(-L_x/2, -L_y/2), fe.Point(L_x/2, L_y/2),
nb_x, nb_y
)
fe.Point(-L_x / 2, -L_y / 2), fe.Point(L_x / 2, L_y / 2), nb_x, nb_y
)
class PeriodicDomain(fe.SubDomain):
def __init__(self, tolerance=fe.DOLFIN_EPS):
......@@ -202,65 +207,64 @@ def test_select_solver():
self.tol = tolerance
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
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:
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:
for i in range(len(x[:])):
y[i] = 1000. * (L_x + L_y)
y[i] = 1000.0 * (L_x + L_y)
pbc = PeriodicDomain()
scalar_fspace = fe.FunctionSpace(
mesh,
fe.FiniteElement('CG', mesh.ufl_cell(), 2),
constrained_domain=pbc)
mesh, fe.FiniteElement("CG", mesh.ufl_cell(), 2), constrained_domain=pbc
)
macro_fields = {
'U': [fe.Expression('0.5*x[0]+x[1]', degree=1)],
'E': [fe.Expression('0.5', degree=0), fe.Expression('1', degree=0)]
"U": [fe.Expression("0.5*x[0]+x[1]", degree=1)],
"E": [fe.Expression("0.5", degree=0), fe.Expression("1", degree=0)],
}
localization_tensors = {
'U': {
'u': [
[fe.interpolate(fe.Constant(1.), scalar_fspace)]
],
"U": {"u": [[fe.interpolate(fe.Constant(1.0), scalar_fspace)]]},
"E": {
"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')]}
exact_sol = fe.project(fe.Expression(
"0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]",
degree=2), scalar_fspace)
function_spaces = {"u": scalar_fspace}
localization_rules = {"u": [("U", "U"), ("E", "E")]}
exact_sol = fe.project(
fe.Expression("0.5*x[0] + x[1] + 0.5*x[1] + 1.*x[0]", degree=2), scalar_fspace
)
reconstr_sol = homog_full.reconstruction(
localization_tensors, macro_fields, function_spaces,
localization_rules=localization_rules, output_request=('u',),
proj_solver="mumps")
reconstr_sol = reconstr_sol['u']
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., abs=1e-14)
localization_tensors,
macro_fields,
function_spaces,
localization_rules=localization_rules,
output_request=("u",),
proj_solver="mumps",
)
reconstr_sol = reconstr_sol["u"]
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__":
logger_root = logging.getLogger()
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.setLevel(logging.DEBUG)
stream_handler.setFormatter(formatter)
logger_root.addHandler(stream_handler)
test_reconstruction_vectors()
test_reconstruction_with_constraint()
\ No newline at end of file
test_reconstruction_with_constraint()
test/geometry_test.py
View file @ ea814ebe
......@@ -6,7 +6,7 @@ Created on Mon Oct 15 11:00:18 2018
"""
import os
import ho_homog.geometry as geo
from .context import ho_homog
import math
import numpy as np
import matplotlib.pyplot as plt
......@@ -16,6 +16,8 @@ from logging.handlers import RotatingFileHandler
import gmsh
geo = ho_homog.geometry
# nice shortcuts
model = gmsh.model
factory = model.occ
......@@ -502,7 +504,6 @@ def test_ll_modif():
os.system("gmsh %s.msh &" % name)
def test_gather_line():
"""
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