Page MenuHomec4science
Files F91927906

fenics_norms.py
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 15, 19:46

fenics_norms.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import fenics as fen
from rrompy.solver.fenics import (L2NormMatrix, H1NormMatrix,
Hminus1NormMatrix, elasticNormMatrix,
elasticDualNormMatrix)
def test_fenics_L2():
V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 3)
u = fen.interpolate(fen.Constant(3.), V)
v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
mass = L2NormMatrix(V)
inner = fen.assemble(fen.dot(u, v) * fen.dx)
assert np.isclose(uv.T.dot(mass.dot(vv)), 2.5, rtol = 1e-8)
assert np.isclose(inner, 2.5, rtol = 1e-8)
def test_fenics_H1():
V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 2)
u = fen.interpolate(fen.Expression("x[0]+exp(x[1])", degree = 4), V)
v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
stiffness = H1NormMatrix(V)
helmholtz = H1NormMatrix(V, 12)
inners = fen.assemble(fen.dot(fen.grad(u), fen.grad(v)) * fen.dx)
innerh = 12 * fen.assemble(fen.dot(u, v) * fen.dx)
assert np.isclose(uv.T.dot(stiffness.dot(vv)), np.exp(1), rtol = 1e-6)
assert np.isclose(inners, np.exp(1), rtol = 1e-6)
assert np.isclose(uv.T.dot(helmholtz.dot(vv)), 5 * np.exp(1) + 14,
rtol = 1e-3)
assert np.isclose(inners + innerh, 5 * np.exp(1) + 14, rtol = 1e-3)
def test_fenics_Hminus1():
V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 2)
u = fen.interpolate(fen.Expression("x[0]+exp(x[1])", degree = 4), V)
v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
energyFull = Hminus1NormMatrix(V, 12)
energyLR = Hminus1NormMatrix(V, 12, compressRank = 20)
assert np.isclose(uv.T.dot(energyFull.dot(vv)),
uv.T.dot(energyLR.dot(vv)), rtol = 1e-2)
assert np.isclose(uv.T.dot(energyFull.dot(vv)), .1641618355, rtol = 1e-6)
def test_fenics_elastic():
V = fen.VectorFunctionSpace(fen.UnitCubeMesh(5, 5, 5), "P", 1)
l_ = 1.
m_ = fen.Expression(".5*x[0]+1.", degree = 1)
u = fen.interpolate(fen.Expression(("exp(x[1])", "x[0]-x[2]", "3."),
degree = 4), V)
v = fen.interpolate(fen.Expression(("x[0]*x[0]+x[2]", "1.", "-1. * x[1]"),
degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
energyMat = elasticNormMatrix(V, l_, m_, 10)
epsilon = lambda f: 0.5 * (fen.grad(f) + fen.nabla_grad(f))
sigma = (l_ * fen.div(u) * fen.Identity(3) + 2. * m_ * epsilon(u))
energy = fen.assemble((10 * fen.dot(u, v)
+ fen.inner(sigma, epsilon(v))) * fen.dx)
assert np.isclose(uv.T.dot(energyMat.dot(vv)), energy, rtol = 1e-8)
def test_fenics_elastic_dual():
V = fen.VectorFunctionSpace(fen.UnitCubeMesh(5, 5, 5), "P", 1)
l_ = 1.
m_ = fen.Expression(".5*x[0]+1.", degree = 1)
u = fen.interpolate(fen.Expression(("exp(x[1])", "x[0]-x[2]", "3."),
degree = 4), V)
v = fen.interpolate(fen.Expression(("x[0]*x[0]+x[2]", "1.", "-1. * x[1]"),
degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
energyFull = elasticDualNormMatrix(V, l_, m_, 10)
energyLR = elasticDualNormMatrix(V, l_, m_, 10, compressRank = 50)
assert np.isclose(uv.T.dot(energyFull.dot(vv)),
uv.T.dot(energyLR.dot(vv)), rtol = 1e-1)
assert np.isclose(uv.T.dot(energyFull.dot(vv)), -.00804628936, rtol = 1e-6)

Event Timeline

c4science · Help