diff --git a/CHANGELOG.md b/CHANGELOG.md
index 89276d1ab..0e9008434 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,16 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+## [Unreleased]
+
+### Changed
+- Renamed skfem.importers to skfem.io.
+
+### Added
+- MortarBasis and MortarPair that support mortar methods in 2D.
+- Serialization of meshes to json-files via skfem.io.json.
+- ElementQuadDG for transforming H1 elements to DG elements.
+
## [0.4.1] - 2020-01-19
### Added
diff --git a/docs/api.rst b/docs/api.rst
index f9503afb3..f9983c6ca 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -169,6 +169,8 @@ Defining a global basis
.. autoclass:: skfem.assembly.InteriorBasis
+.. automethod:: skfem.assembly.InteriorBasis.refinterp
+
.. autoclass:: skfem.assembly.FacetBasis
Assembling matrices
diff --git a/docs/conf.py b/docs/conf.py
index fc5b10195..e8e3b6326 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -19,7 +19,7 @@
# -- Project information -----------------------------------------------------
project = 'scikit-fem'
-copyright = '2018-2019, Tom Gustafsson, Geordie McBain'
+copyright = '2018-20, Tom Gustafsson, Geordie McBain'
author = 'Tom Gustafsson, Geordie McBain'
diff --git a/docs/learning.rst b/docs/examples.rst
similarity index 97%
rename from docs/learning.rst
rename to docs/examples.rst
index f7432a536..9dccffffa 100644
--- a/docs/learning.rst
+++ b/docs/examples.rst
@@ -15,13 +15,14 @@ Poisson equation
examples/ex14.rst
examples/ex12.rst
examples/ex13.rst
- examples/ex15.rst
examples/ex26.rst
examples/ex09.rst
examples/ex07.rst
examples/ex05.rst
examples/ex10.rst
examples/ex22.rst
+ examples/ex15.rst
+ examples/ex06.rst
Heat transfer
#############
@@ -51,11 +52,12 @@ Structural mechanics
.. toctree::
- examples/ex03.rst
- examples/ex11.rst
- examples/ex08.rst
examples/ex02.rst
examples/ex21.rst
+ examples/ex04.rst
+ examples/ex11.rst
+ examples/ex03.rst
+ examples/ex08.rst
Others
######
@@ -64,4 +66,3 @@ Others
examples/ex23.rst
examples/ex16.rst
- examples/ex06.rst
diff --git a/docs/examples/Makefile b/docs/examples/Makefile
index d71a6dcba..4ff37034d 100644
--- a/docs/examples/Makefile
+++ b/docs/examples/Makefile
@@ -1,5 +1,6 @@
all: ex01_solution.png \
ex02_solution.png \
+ ex04_solution.png \
ex05_solution.png \
ex06_solution.png \
ex17_solution.png \
diff --git a/docs/examples/ex01.py b/docs/examples/ex01.py
index ff85d171f..a45a75668 100644
--- a/docs/examples/ex01.py
+++ b/docs/examples/ex01.py
@@ -8,11 +8,11 @@
@bilinear_form
def laplace(u, du, v, dv, w):
- return du[0]*dv[0] + du[1]*dv[1]
+ return du[0] * dv[0] + du[1] * dv[1]
@linear_form
def load(v, dv, w):
- return 1.0*v
+ return 1. * v
A = asm(laplace, basis)
b = asm(load, basis)
diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py
index 00045dd5c..474632caf 100644
--- a/docs/examples/ex04.py
+++ b/docs/examples/ex04.py
@@ -1,133 +1,213 @@
-"""
-Author: kinnala
-
-Solve a linearized contact problem between two linear elastic
-bodies using the Nitsche's method.
-"""
-
from skfem import *
from skfem.models.elasticity import linear_elasticity
import numpy as np
-
-m = MeshTri()
-m.refine(3)
-M = MeshLine(np.linspace(0, 1, 6))
-M = M*M
-M = M._splitquads()
+from skfem.io import from_meshio
+from skfem.io.json import from_file, to_file
+
+
+# create meshes
+try:
+ m = from_file("docs/examples/ex04_mesh.json")
+except FileNotFoundError:
+ from pygmsh import generate_mesh
+ from pygmsh.built_in import Geometry
+ geom = Geometry()
+ points = []
+ lines = []
+ points.append(geom.add_point([0., 0., 0.], .1))
+ points.append(geom.add_point([0., 1., 0.], .1))
+ points.append(geom.add_point([0.,-1., 0.], .1))
+ lines.append(geom.add_circle_arc(points[2], points[0], points[1]))
+ geom.add_physical(lines[-1], 'contact')
+ lines.append(geom.add_line(points[1], points[2]))
+ geom.add_physical(lines[-1], 'dirichlet')
+ geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain')
+ m = from_meshio(generate_mesh(geom, dim=2))
+ to_file(m, "docs/examples/ex04_mesh.json")
+
+M = MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10))
M.translate((1.0, 0.0))
+M.refine()
-map = MappingAffine(m)
-Map = MappingAffine(M)
-e1 = ElementTriP1()
-e = ElementVectorH1(e1)
-ib = InteriorBasis(m, e, map, 2)
-Ib = InteriorBasis(M, e, Map, 2)
+# define elements and bases
+e1 = ElementTriP2()
+e = ElementVectorH1(e1)
-def rule(x, y):
- return (x==1.0)
+E1 = ElementQuad2()
+E = ElementVectorH1(E1)
-def param(x, y):
- return y
+ib = InteriorBasis(m, e, intorder=4)
+Ib = InteriorBasis(M, E, intorder=4)
-mortar = InterfaceMesh1D(m, M, rule, param)
+mappings = MortarPair.init_2D(m, M,
+ m.boundaries['contact'],
+ M.facets_satisfying(lambda x: x[0] == 1.0),
+ np.array([0.0, 1.0]))
-mb = {}
-mortar_map = MappingAffine(mortar)
-mb[0] = FacetBasis(mortar, e, mortar_map, 2, side=0)
-mb[1] = FacetBasis(mortar, e, mortar_map, 2, side=1)
+mb = [
+ MortarBasis(m, e, mapping = mappings[0], intorder=4),
+ MortarBasis(M, E, mapping = mappings[1], intorder=4),
+]
+# define bilinear forms
E1 = 1000.0
E2 = 1000.0
nu1 = 0.3
nu2 = 0.3
-Mu1 = E1/(2.0*(1.0 + nu1))
-Mu2 = E2/(2.0*(1.0 + nu2))
+Mu1 = E1 / (2. * (1. + nu1))
+Mu2 = E2 / (2. * (1. + nu2))
-Lambda1 = E1*nu1/((1.0 + nu1)*(1.0 - 2.0*nu1))
-Lambda2 = E2*nu2/((1.0 + nu2)*(1.0 - 2.0*nu2))
+Lambda1 = E1 * nu1 / ((1. + nu1) * (1. - 2. * nu1))
+Lambda2 = E2 * nu2 / ((1. + nu2) * (1. - 2. * nu2))
Mu = Mu1
Lambda = Lambda1
-weakform1 = linear_elasticity(Lambda=Lambda1, Mu=Mu1)
-weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu2)
+weakform1 = linear_elasticity(Lambda=Lambda, Mu=Mu)
+weakform2 = linear_elasticity(Lambda=Lambda, Mu=Mu)
+
-alpha = 1
+alpha = 1000
+limit = 0.3
+
+# assemble the stiffness matrices
K1 = asm(weakform1, ib)
K2 = asm(weakform2, Ib)
-L = 0
+K = [[K1, 0.], [0., K2]]
+f = [None] * 2
+
+def tr(T):
+ return T[0, 0] + T[1, 1]
+
+def C(T):
+ return np.array([[2. * Mu * T[0, 0] + Lambda * tr(T), 2. * Mu * T[0, 1]],
+ [2. * Mu * T[1, 0], 2. * Mu * T[1, 1] + Lambda * tr(T)]])
+
+def Eps(dw):
+ return np.array([[dw[0][0], .5 * (dw[0][1] + dw[1][0])],
+ [.5 * (dw[1][0] + dw[0][1]), dw[1][1]]])
+
+
+# assemble the mortar matrices
for i in range(2):
for j in range(2):
@bilinear_form
def bilin_penalty(u, du, v, dv, w):
n = w.n
- ju = (-1.0)**i*(u[0]*n[0] + u[1]*n[1])
- jv = (-1.0)**j*(v[0]*n[0] + v[1]*n[1])
-
- def tr(T):
- return T[0, 0] + T[1, 1]
-
- def C(T):
- return np.array([[2*Mu*T[0, 0] + Lambda*tr(T), 2*Mu*T[0, 1]],
- [2*Mu*T[1, 0], 2*Mu*T[1, 1] + Lambda*tr(T)]])
-
- def Eps(dw):
- return np.array([[dw[0][0], 0.5*(dw[0][1] + dw[1][0])],
- [0.5*(dw[1][0] + dw[0][1]), dw[1][1]]])
- mu = 0.5*(n[0]*C(Eps(du))[0, 0]*n[0] +\
- n[0]*C(Eps(du))[0, 1]*n[1] +\
- n[1]*C(Eps(du))[1, 0]*n[0] +\
- n[1]*C(Eps(du))[1, 1]*n[1])
- mv = 0.5*(n[0]*C(Eps(dv))[0, 0]*n[0] +\
- n[0]*C(Eps(dv))[0, 1]*n[1] +\
- n[1]*C(Eps(dv))[1, 0]*n[0] +\
- n[1]*C(Eps(dv))[1, 1]*n[1])
+ ju = (-1.) ** i * (u[0] * n[0] + u[1] * n[1])
+ jv = (-1.) ** j * (v[0] * n[0] + v[1] * n[1])
+ mu = .5 * (n[0] * C(Eps(du))[0, 0] * n[0] +
+ n[0] * C(Eps(du))[0, 1] * n[1] +
+ n[1] * C(Eps(du))[1, 0] * n[0] +
+ n[1] * C(Eps(du))[1, 1] * n[1])
+ mv = .5 * (n[0] * C(Eps(dv))[0, 0] * n[0] +
+ n[0] * C(Eps(dv))[0, 1] * n[1] +
+ n[1] * C(Eps(dv))[1, 0] * n[0] +
+ n[1] * C(Eps(dv))[1, 1] * n[1])
h = w.h
- return 1.0/(alpha*h)*ju*jv - mu*jv - mv*ju
-
- L = asm(bilin_penalty, mb[i], mb[j]) + L
-
-@linear_form
-def load(v, dv, w):
- return -50*v[1]
-
-f1 = asm(load, ib)
-f2 = np.zeros(K2.shape[0])
+ return ((1. / (alpha * h) * ju * jv - mu * jv - mv * ju)
+ * (np.abs(w.x[1]) <= limit))
+
+ K[j][i] += asm(bilin_penalty, mb[i], mb[j])
+
+ @linear_form
+ def lin_penalty(v, dv, w):
+ n = w.n
+ jv = (-1.) ** i * (v[0] * n[0] + v[1] * n[1])
+ mv = .5 * (n[0] * C(Eps(dv))[0, 0] * n[0] +
+ n[0] * C(Eps(dv))[0, 1] * n[1] +
+ n[1] * C(Eps(dv))[1, 0] * n[0] +
+ n[1] * C(Eps(dv))[1, 1] * n[1])
+ h = w.h
+ def gap(x):
+ return (1. - np.sqrt(1. - x[1] ** 2))
+ return ((1. / (alpha * h) * gap(w.x) * jv - gap(w.x) * mv)
+ * (np.abs(w.x[1]) <= limit))
+
+ f[i] = asm(lin_penalty, mb[i])
import scipy.sparse
-K = (scipy.sparse.bmat([[K1, None],[None, K2]]) + L).tocsr()
+K = (scipy.sparse.bmat(K)).tocsr()
+
+# set boundary conditions and solve
i1 = np.arange(K1.shape[0])
i2 = np.arange(K2.shape[0]) + K1.shape[0]
-D1 = ib.get_dofs(lambda x, y: x==0.0).all()
-D2 = Ib.get_dofs(lambda x, y: x==2.0).all()
+D1 = ib.get_dofs(lambda x: x[0] == 0.0).all()
+D2 = Ib.get_dofs(lambda x: x[0] == 2.0).all()
x = np.zeros(K.shape[0])
-f = np.hstack((f1, f2))
+f = np.hstack((f[0], f[1]))
x = np.zeros(K.shape[0])
D = np.concatenate((D1, D2 + ib.N))
I = np.setdiff1d(np.arange(K.shape[0]), D)
-x[I] = solve(*condense(K, f, I=I))
+x[ib.get_dofs(lambda x: x[0] == 0.0).nodal['u^1']] = 0.1
+x[ib.get_dofs(lambda x: x[0] == 0.0).facet['u^1']] = 0.1
+x = solve(*condense(K, f, I=I, x=x))
+
+
+# create a displaced mesh
sf = 1
-m.p[0, :] = m.p[0, :] + sf*x[i1][ib.nodal_dofs[0, :]]
-m.p[1, :] = m.p[1, :] + sf*x[i1][ib.nodal_dofs[1, :]]
+m.p[0, :] = m.p[0, :] + sf * x[i1][ib.nodal_dofs[0, :]]
+m.p[1, :] = m.p[1, :] + sf * x[i1][ib.nodal_dofs[1, :]]
-M.p[0, :] = M.p[0, :] + sf*x[i2][Ib.nodal_dofs[0, :]]
-M.p[1, :] = M.p[1, :] + sf*x[i2][Ib.nodal_dofs[1, :]]
+M.p[0, :] = M.p[0, :] + sf * x[i2][Ib.nodal_dofs[0, :]]
+M.p[1, :] = M.p[1, :] + sf * x[i2][Ib.nodal_dofs[1, :]]
-if __name__ == "__main__":
- from skfem.visuals.matplotlib import draw, show
+# post processing
+s, S = {}, {}
+e_dg = ElementTriDG(ElementTriP1())
+E_dg = ElementQuadDG(ElementQuad1())
+
+for itr in range(2):
+ for jtr in range(2):
+ @bilinear_form
+ def proj_cauchy(u, du, v, dv, w):
+ return C(Eps(du))[itr, jtr] * v
+
+ @bilinear_form
+ def mass(u, du, v, dv, w):
+ return u * v
+
+ ib_dg = InteriorBasis(m, e_dg, intorder=4)
+ Ib_dg = InteriorBasis(M, E_dg, intorder=4)
+
+ s[itr, jtr] = solve(asm(mass, ib_dg),
+ asm(proj_cauchy, ib, ib_dg) @ x[i1])
+ S[itr, jtr] = solve(asm(mass, Ib_dg),
+ asm(proj_cauchy, Ib, Ib_dg) @ x[i2])
+
+s[2, 2] = nu1 * (s[0, 0] + s[1, 1])
+S[2, 2] = nu2 * (S[0, 0] + S[1, 1])
+
+vonmises1 = np.sqrt(.5 * ((s[0, 0] - s[1, 1]) ** 2 +
+ (s[1, 1] - s[2, 2]) ** 2 +
+ (s[2, 2] - s[0, 0]) ** 2 +
+ 6. * s[0, 1]**2))
+
+vonmises2 = np.sqrt(.5 * ((S[0, 0] - S[1, 1]) ** 2 +
+ (S[1, 1] - S[2, 2]) ** 2 +
+ (S[2, 2] - S[0, 0]) ** 2 +
+ 6. * S[0, 1]**2))
+
+
+if __name__ == "__main__":
+ from os.path import splitext
+ from sys import argv
+ from skfem.visuals.matplotlib import *
- ax = draw(m)
+ ax = plot(ib_dg, vonmises1, shading='gouraud')
+ draw(m, ax=ax)
+ plot(Ib_dg, vonmises2, ax=ax, Nrefs=3, shading='gouraud')
draw(M, ax=ax)
- show()
+ savefig(splitext(argv[0])[0] + '_solution.png')
diff --git a/docs/examples/ex04.rst b/docs/examples/ex04.rst
new file mode 100644
index 000000000..0a7b3fb82
--- /dev/null
+++ b/docs/examples/ex04.rst
@@ -0,0 +1,70 @@
+.. _mortaring:
+
+Contact problem
+---------------
+
+.. note::
+
+ This example requires the external package `pygmsh
+ `_.
+
+Mortar methods allow setting interface conditions on non-matching meshes.
+They are useful also when solving variational inequalities such as
+`elastic contact problems `_.
+
+This example solves the first contact iteration for the following prototype
+contact problem: find :math:`\boldsymbol{u}_i : \Omega_i \rightarrow
+\mathbb{R}^2`, :math:`i = 1,2`, such that
+
+
+.. math::
+
+ \begin{aligned}
+ \boldsymbol{\mathrm{div}}\,\boldsymbol{\sigma}_i(\boldsymbol{u}_i)&=\boldsymbol{0} \quad && \text{in $\Omega_i$,} \\
+ \boldsymbol{u}_1&=(0.1, 0) \quad && \text{on $\Gamma_{D,1}$,} \\
+ \boldsymbol{u}_2&=\boldsymbol{0} \quad && \text{on $\Gamma_{D,2}$,} \\
+ \boldsymbol{\sigma}_2(\boldsymbol{u}_2) \boldsymbol{n}_2 &=\boldsymbol{0} \quad && \text{on $\Gamma_{N,2}$,} \\
+ \boldsymbol{\sigma}_{i,t}(\boldsymbol{u}_i) &= \boldsymbol{0} && \text{on $\Gamma$,} \\
+ \sigma_{1,n}(\boldsymbol{u}_1(\boldsymbol{\gamma}(\boldsymbol{x})) - \sigma_{2,n}(\boldsymbol{u}_2)&=0 && \text{on $\Gamma$,} \\
+ [[u_n]] - g &\geq 0 && \text{on $\Gamma$,} \\
+ \sigma_{2,n}(\boldsymbol{u}_2)&\leq 0 && \text{on $\Gamma$,} \\
+ ([[u_n]] - g) \sigma_{2,n}(\boldsymbol{u}_2) &= 0 && \text{on $\Gamma$,}
+ \end{aligned}
+
+
+where
+
+* :math:`\Omega_1 = \{ (x, y) : x^2 + y^2 < 1 \} \setminus \{ (x, y) : x < 0\}`,
+* :math:`\Omega_2 = (1, 2) \times (-1, 1)`,
+* :math:`\Gamma_{D,1} = \{ (x, y) \in \Omega_1 : x=0 \}`,
+* :math:`\Gamma_{D,2} = \{ (x, y) \in \Omega_2 : x=2 \}`,
+* :math:`\Gamma = \{ (x, y) \in \Omega_2 : x=1 \}`,
+* :math:`g((x,y)) = 1 - \sqrt{1 - y^2}`,
+* :math:`\boldsymbol{\gamma} : \Gamma \rightarrow \{ (x, y) \in \partial
+ \Omega_1 : x > 0 \}`, :math:`\boldsymbol{\gamma}((x,y)) = (g(x-1)+1, y)`,
+* :math:`[[u_n]] = \boldsymbol{u}_1(\boldsymbol{\gamma}(\boldsymbol{x})) \cdot \boldsymbol{n} - \boldsymbol{u}_2(\boldsymbol{x}) \cdot \boldsymbol{n}`,
+
+and the directions for evaluating :math:`[[u_n]]`, :math:`\sigma_{1,n}` and
+:math:`\boldsymbol{\sigma}_{1,t}` are defined as :math:`\boldsymbol{n}=(1,0)`
+and :math:`\boldsymbol{t}=(0,1)`.
+This is a nonlinear problem since we do not know a priori which subset
+:math:`\Gamma_C \subset \Gamma` satisfies :math:`([[u_n]] - g)|_{\Gamma_C} = 0`.
+
+.. note::
+
+ The example solves a simplified prototype problem.
+ Instead of iterating for the true contact boundary,
+ we solve a single contact iteration (a linear problem) with the initial
+ guess :math:`\{ (x, y) \in \Gamma : |y| < 0.1 \}`.
+ Solving a real contact problem involves repeatedly solving and guessing a new
+ candidate boundary :math:`\Gamma_C` until convergence.
+ Extending this example should be straightforward.
+
+.. figure:: ex04_solution.png
+
+ The von Mises stress in a displaced geometry.
+
+The complete source code reads as follows:
+
+.. literalinclude:: ex04.py
+ :linenos:
diff --git a/docs/examples/ex04_mesh.json b/docs/examples/ex04_mesh.json
new file mode 100644
index 000000000..8064125bd
--- /dev/null
+++ b/docs/examples/ex04_mesh.json
@@ -0,0 +1 @@
+{"p": [[0.0, 1.0], [0.0, -1.0], [0.09801714055591637, -0.9951847266499028], [0.19509032241883761, -0.9807852803231266], [0.2902846778931259, -0.9569403355384724], [0.38268343326549487, -0.9238795321383267], [0.47139673797238385, -0.8819212637355989], [0.5555702343496238, -0.8314696114138532], [0.6343932856511478, -0.7730104521419754], [0.7071067827656586, -0.7071067796074364], [0.773010454973133, -0.6343932822013723], [0.8314696138527624, -0.5555702306995381], [0.8819212658270285, -0.47139673405959426], [0.9238795338254099, -0.38268342919251586], [0.9569403368184075, -0.2902846736737454], [0.9807852811883714, -0.1950903180689581], [0.995184727094887, -0.09801713603791533], [1.0, 4.5837337757057685e-09], [0.9951847262476915, 0.09801714463963587], [0.980785279615548, 0.19509032597607545], [0.9569403346479524, 0.290284680828777], [0.9238795311905263, 0.38268343555368756], [0.8819212628015901, 0.4713967397197912], [0.8314696106349614, 0.5555702355153178], [0.7730104516368775, 0.6343932862666112], [0.7071067794805702, 0.7071067828925248], [0.634393282512301, 0.7730104547179604], [0.5555702315457777, 0.8314696132873234], [0.4713967355549627, 0.8819212650277374], [0.38268343125990784, 0.9238795329690681], [0.2902846763567968, 0.9569403360045128], [0.195090321428782, 0.9807852805200609], [0.09801714007694617, 0.9951847266970772], [0.0, 0.9000000000002774], [0.0, 0.8000000000005548], [0.0, 0.7000000000008322], [0.0, 0.6000000000011096], [0.0, 0.500000000001387], [0.0, 0.40000000000166436], [0.0, 0.3000000000019418], [0.0, 0.20000000000221918], [0.0, 0.10000000000249643], [0.0, 2.752242878045763e-12], [0.0, -0.09999999999750342], [0.0, -0.19999999999778084], [0.0, -0.29999999999805826], [0.0, -0.39999999999833546], [0.0, -0.4999999999986129], [0.0, -0.5999999999988903], [0.0, -0.699999999999168], [0.0, -0.7999999999994454], [0.0, -0.8999999999997226], [0.9159764153258458, 0.05383987386802238], [0.8314138091584239, -0.39785805095354165], [0.8264043258240841, 0.39254386345989056], [0.08657154386342972, -0.04993181939059333], [0.08825881533998217, 0.24854065189332641], [0.0866467339337712, -0.350125666409955], [0.6151960220089083, 0.6914850111403807], [0.611656851878296, -0.6767753615474275], [0.08792437843765298, 0.546553776392746], [0.40214879730821057, -0.8329018889819535], [0.0864975927752666, -0.5501598487639185], [0.4031239663030456, 0.8299929282327413], [0.08552443345288481, 0.7412874524475859], [0.08387010326216801, -0.7452678551627052], [0.9064888004049059, -0.13475830751070036], [0.883006687779775, 0.21591742664103838], [0.7318295262956729, 0.5462376817273685], [0.7310857439949313, -0.5534884214698972], [0.22203717119110486, 0.8864220786001916], [0.222037172228374, -0.88642207782408], [0.12735200151713993, 0.8953724394140379], [0.1264872216201089, -0.8973542598638573], [0.08603181563097448, -0.6492168486242395], [0.1728569295463326, -0.5999642076514183], [0.1729345952627196, -0.500480145586956], [0.2591331121034306, -0.5504140700880753], [0.2592714087128541, -0.4510637781789997], [0.3451538888162647, -0.5008525665656047], [0.3453660173120208, -0.4016689274236115], [0.43112082579940325, -0.45140415758047114], [0.4313132869549403, -0.35228060068609146], [0.3456417826557642, -0.30202707863747263], [0.4316914876509484, -0.25313887750341396], [0.5172439124243171, -0.30297215045699327], [0.5175060234485399, -0.20398089749625326], [0.4320078005036593, -0.1541618871667654], [0.5178249360980043, -0.10497328859987322], [0.4323517746522692, -0.05513829699300258], [0.5181810627769347, -0.005926661553542319], [0.4326873030792038, 0.0439068356495514], [0.5185129534634814, 0.0931212502063908], [0.43303358988383006, 0.14295033231019513], [0.5188407112939315, 0.19216073999298827], [0.4333907822760996, 0.2419806037163972], [0.5191570773658916, 0.29117889053487644], [0.4336167327285912, 0.34098097416379747], [0.3466297470498466, 0.09502752487561306], [0.6036466835467159, -0.05588497882022823], [0.6043483291878866, 0.14226361071690696], [0.6030083111892219, -0.253739769204638], [0.5194358558074463, 0.3901995369490897], [0.4339131655358476, 0.44001601280296204], [0.34799605056576544, 0.3907547089718833], [0.3483736150587413, 0.4897859002884928], [0.3461095808402675, -0.10314722872110887], [0.3451380056642011, -0.5999795967606039], [0.43423779016255415, 0.5389955255642274], [0.3487855783748346, 0.5887556881868931], [0.6027356100339158, -0.35253659554114325], [0.6048915795006872, 0.34040851731280114], [0.26231771840506496, 0.4405891472347052], [0.26200491468128767, 0.34136888541754196], [0.4341509042932531, 0.6378572118248705], [0.17613909944252154, 0.39128549747516345], [0.6867867143117476, -0.3005776176954207], [0.6883938442262401, -0.20422841407407338], [0.7887504482032948, -0.24363163143334435], [0.6039194807805792, 0.4371233138451532], [0.34717469685206775, 0.6842928036187731], [0.2630163234652642, 0.6380173006925967], [0.26231461886291657, 0.7347372834745751], [0.6821568960369653, -0.38932236004593684], [0.6024613510918597, -0.4500576975192292], [0.7718417301410816, -0.14476793206076038], [0.6840585160728114, 0.38136305112776203], [0.6872609976885783, 0.28791528538360167], [0.5168363646187115, 0.5833344105436389], [0.17550407052083822, 0.6900813696897046], [0.17619557798382957, 0.4914856785795687], [0.42016802198142106, 0.7346497809142203], [0.17302511054291092, -0.40062265779056283], [0.17237477181374233, -0.6988080454878537], [0.1714952098378757, -0.8012321435395398], [0.25599371682296734, -0.7410015264285403], [0.5190057914509976, 0.4879530306260381], [0.618290339602385, 0.536875949487155], [0.4308956703486012, -0.5504811711457109], [0.43130650380920504, -0.6485484160152994], [0.5175687390320705, -0.6021217093581201], [0.5092206012248724, -0.7035085970951118], [0.620737462322611, -0.555765633253074], [0.3429556442908701, -0.6929402428430802], [0.25848406726552225, -0.6479093927842566], [0.6885382570987393, -0.10392056646537914], [0.769286205655087, -0.05057534362618522], [0.6040111157907584, 0.043242214936212164], [0.6897389251565653, 0.09246051544826223], [0.6895451688254398, 0.1909341621919352], [0.7853297010802, 0.13990094156251878], [0.7756810503792961, 0.2413610566339871], [0.7701883492540708, 0.03853734593193141], [0.5190243392045759, -0.5017121515376488], [0.6040073106437357, 0.2408102010221849], [0.687568256083656, -0.006023468765897795], [0.08680748951240319, 0.14978745087550502], [0.17485631474456298, 0.19832056843973647], [0.1730517611561072, 0.10005310843209428], [0.26087327445847697, 0.14679457604766757], [0.08657886954653793, 0.04998478457738213], [0.17291704371860217, -0.00010105186068302731], [0.1729970659840663, -0.10002826078482492], [0.2593646476016979, -0.04985929515148815], [0.08660511306657999, -0.1498281482654842], [0.17305207548825177, -0.20001247954407836], [0.08658138073658056, -0.24999061371589165], [0.17306526794695942, -0.3003972076357125], [0.2594694617813657, -0.25019468739155504], [0.517316554251502, -0.4018272255535962], [0.08648158353870185, -0.4502095033908465], [0.6031530092679103, -0.1544546524434075], [0.26271089280937426, 0.539695638559989], [0.25930650815864587, -0.35099572284298564], [0.3481417000647074, 0.29225483630730303], [0.08775490432998115, 0.6445771191516094], [0.17312433591087684, 0.7949143361266013], [0.08917340860778944, 0.4445319633762133], [0.17610593477516273, 0.5911079955205527], [0.34591484982492143, -0.20259666030335644], [0.2596845571558199, 0.04809501117136811], [0.259197333421417, -0.15100131147531112], [0.34645554142539703, -0.004036758239607918], [0.3476513202248709, 0.19411470667553155], [0.5165083519158323, 0.6874008844465628], [0.26285614408485514, 0.24520073416291943], [0.8702436047190465, 0.12724803991148598], [0.1765715058218143, 0.29520431385172685], [0.08940011612821333, 0.34630648225343424], [0.4166298503101692, -0.7429585203863014], [0.8549693515455798, 0.30930121195949895], [0.8950418208143169, -0.2241963012461981], [0.860025803539879, -0.3124749658989302], [0.6724940254177364, 0.6208040077969013], [0.7854994399143052, 0.4754049055547044], [0.9215823721068173, -0.045274434635522595], [0.785420238160406, -0.4789576989502909], [0.31501693368973693, 0.863121801168385], [0.31408164224217594, -0.8662391964201028], [0.688719459186926, -0.6255058956158415], [0.844562532980393, -0.007366724481464558], [0.6969550069593997, 0.46476187067423613], [0.7641662461091518, 0.3352234361081791], [0.766062195642033, -0.3394573773342276], [0.6967061893284834, -0.4732097765718286], [0.48759115339420067, 0.7864769275402721], [0.33363382533391917, 0.7769765419565695], [0.4861251931921746, -0.7917065301079622], [0.5603470596005116, 0.7537033611358295], [0.25400471431432753, 0.8150646543307231], [0.3311957483405575, -0.7829906385661866], [0.8445851642438251, -0.08836500410921283], [0.5617849379877878, -0.7574792990349553], [0.07305726944797225, 0.8238422029156591], [0.07358418293812316, -0.8256623870972735], [0.5883410243610833, 0.6196215619927677], [0.2566752832473187, -0.8195973033555625], [0.841341592761485, -0.16714383527204318], [0.7514167069759504, 0.40985942538495446], [0.7553021899259311, -0.4137557280795954], [0.8372601206719111, 0.0704318953584988], [0.9360577467536026, 0.13890639314725461]], "t": [[32, 33, 72], [2, 51, 73], [67, 150, 151], [67, 150, 186], [0, 32, 33], [1, 2, 51], [124, 142, 153], [31, 70, 72], [3, 71, 73], [31, 32, 72], [2, 3, 73], [68, 137, 193], [70, 72, 176], [59, 140, 142], [118, 191, 192], [68, 137, 201], [140, 142, 153], [20, 21, 190], [59, 140, 141], [118, 191, 217], [124, 142, 204], [71, 73, 134], [72, 176, 213], [69, 142, 204], [131, 184, 205], [184, 205, 208], [37, 60, 177], [60, 130, 177], [69, 142, 199], [20, 67, 190], [21, 54, 190], [64, 129, 175], [64, 129, 176], [118, 192, 203], [129, 175, 178], [38, 39, 188], [60, 175, 178], [196, 204, 219], [37, 38, 177], [60, 130, 178], [98, 159, 180], [67, 151, 190], [167, 168, 173], [83, 168, 173], [106, 163, 181], [98, 159, 183], [98, 180, 182], [38, 177, 188], [163, 180, 182], [69, 196, 204], [57, 132, 167], [132, 167, 173], [162, 163, 181], [146, 152, 155], [158, 159, 180], [106, 163, 182], [127, 149, 151], [127, 149, 154], [135, 143, 144], [106, 179, 181], [75, 133, 144], [133, 135, 144], [57, 132, 170], [165, 167, 168], [99, 145, 155], [145, 146, 155], [45, 57, 166], [124, 153, 169], [57, 166, 167], [35, 64, 175], [165, 166, 167], [168, 179, 181], [65, 133, 134], [158, 160, 161], [83, 168, 179], [102, 119, 136], [128, 136, 137], [62, 74, 75], [99, 145, 171], [119, 136, 137], [133, 134, 135], [75, 77, 144], [65, 74, 133], [100, 147, 148], [96, 111, 154], [48, 62, 74], [107, 138, 139], [111, 127, 154], [148, 152, 155], [147, 148, 155], [74, 75, 133], [92, 100, 147], [44, 45, 166], [81, 153, 169], [102, 103, 136], [79, 107, 138], [88, 99, 171], [46, 57, 170], [156, 158, 160], [100, 149, 154], [138, 140, 153], [62, 76, 170], [41, 156, 160], [76, 132, 170], [138, 139, 140], [77, 79, 107], [139, 140, 141], [75, 76, 77], [62, 75, 76], [43, 44, 164], [81, 138, 153], [110, 124, 169], [49, 65, 74], [77, 78, 79], [114, 120, 131], [79, 81, 138], [100, 148, 149], [94, 96, 154], [108, 114, 128], [56, 156, 157], [85, 110, 169], [111, 119, 126], [161, 162, 163], [157, 158, 159], [148, 149, 150], [34, 35, 64], [101, 110, 116], [49, 50, 65], [55, 160, 161], [96, 102, 111], [102, 111, 119], [110, 116, 123], [108, 109, 114], [79, 80, 81], [99, 147, 155], [44, 164, 166], [85, 101, 110], [90, 92, 147], [92, 94, 100], [85, 86, 101], [158, 161, 180], [116, 117, 118], [47, 48, 62], [149, 150, 151], [88, 90, 99], [94, 100, 154], [156, 157, 158], [55, 161, 162], [103, 108, 136], [82, 84, 85], [76, 77, 78], [107, 139, 143], [90, 99, 147], [36, 37, 60], [40, 56, 156], [148, 150, 152], [42, 55, 160], [78, 79, 80], [109, 114, 120], [117, 118, 125], [101, 116, 117], [81, 82, 169], [48, 49, 74], [84, 85, 86], [39, 40, 56], [111, 126, 127], [96, 97, 102], [103, 105, 108], [45, 46, 57], [80, 81, 82], [42, 43, 55], [101, 117, 171], [94, 95, 96], [162, 164, 165], [117, 145, 171], [108, 128, 136], [41, 42, 160], [92, 93, 94], [82, 85, 169], [40, 41, 156], [90, 91, 92], [43, 55, 164], [105, 108, 109], [110, 123, 124], [97, 102, 103], [47, 62, 170], [88, 89, 90], [86, 87, 88], [125, 145, 146], [86, 88, 171], [77, 107, 144], [82, 83, 84], [46, 47, 170], [95, 96, 97], [162, 165, 181], [55, 162, 164], [86, 101, 171], [93, 94, 95], [164, 165, 166], [64, 176, 213], [120, 121, 122], [91, 92, 93], [117, 125, 145], [89, 90, 91], [84, 86, 87], [87, 88, 89], [15, 16, 66], [19, 20, 67], [103, 104, 105], [80, 82, 83], [23, 24, 68], [76, 78, 132], [21, 22, 54], [17, 18, 52], [10, 11, 69], [12, 13, 53], [25, 26, 58], [35, 36, 175], [97, 103, 104], [8, 9, 59], [109, 120, 121], [137, 193, 215], [107, 143, 144], [91, 93, 98], [112, 115, 130], [28, 29, 63], [87, 89, 106], [5, 6, 61], [104, 112, 113], [122, 129, 176], [112, 130, 172], [105, 109, 172], [78, 80, 173], [104, 105, 112], [78, 132, 173], [80, 83, 173], [95, 97, 174], [112, 113, 115], [36, 60, 175], [97, 104, 174], [121, 122, 129], [109, 121, 172], [30, 31, 70], [3, 4, 71], [104, 113, 174], [105, 112, 172], [83, 84, 179], [115, 130, 177], [39, 56, 188], [87, 106, 179], [91, 98, 182], [89, 106, 182], [56, 157, 187], [93, 98, 183], [121, 172, 178], [95, 174, 183], [165, 168, 181], [157, 159, 185], [84, 87, 179], [121, 129, 178], [130, 172, 178], [161, 163, 180], [33, 72, 213], [51, 73, 214], [89, 91, 182], [114, 128, 184], [93, 95, 183], [56, 187, 188], [157, 185, 187], [141, 189, 207], [73, 134, 214], [174, 183, 185], [114, 131, 184], [141, 207, 212], [139, 143, 189], [159, 183, 185], [68, 194, 201], [143, 189, 210], [194, 201, 218], [113, 185, 187], [113, 174, 185], [113, 115, 187], [53, 192, 203], [115, 187, 188], [115, 177, 188], [59, 142, 199], [197, 206, 209], [118, 125, 217], [119, 137, 201], [139, 141, 189], [128, 184, 215], [61, 189, 210], [63, 131, 206], [116, 123, 203], [58, 184, 215], [120, 131, 206], [65, 134, 214], [63, 197, 206], [33, 34, 213], [50, 51, 214], [125, 146, 211], [70, 176, 209], [14, 15, 191], [16, 17, 195], [22, 23, 194], [13, 14, 192], [24, 25, 193], [11, 12, 196], [15, 66, 191], [16, 66, 195], [23, 68, 194], [22, 54, 194], [24, 68, 193], [17, 52, 195], [11, 69, 196], [13, 53, 192], [12, 53, 196], [25, 58, 193], [29, 30, 197], [4, 5, 198], [150, 152, 220], [29, 63, 197], [5, 61, 198], [123, 204, 219], [30, 70, 197], [4, 71, 198], [127, 151, 202], [10, 69, 199], [9, 59, 199], [126, 127, 202], [150, 186, 220], [63, 131, 205], [59, 141, 212], [28, 63, 205], [146, 152, 200], [26, 58, 208], [6, 61, 207], [9, 10, 199], [19, 67, 221], [18, 52, 221], [8, 59, 212], [119, 126, 201], [14, 191, 192], [58, 184, 208], [125, 211, 217], [151, 190, 202], [67, 186, 221], [52, 186, 221], [61, 189, 207], [116, 118, 203], [123, 124, 204], [61, 198, 210], [71, 134, 216], [27, 28, 205], [120, 122, 206], [134, 135, 216], [52, 195, 200], [122, 176, 209], [135, 143, 210], [58, 193, 215], [126, 201, 218], [198, 210, 216], [123, 203, 219], [6, 7, 207], [26, 27, 208], [54, 190, 202], [122, 206, 209], [128, 137, 215], [7, 8, 212], [34, 64, 213], [50, 65, 214], [18, 19, 221], [66, 195, 211], [27, 205, 208], [7, 207, 212], [70, 197, 209], [52, 186, 220], [195, 200, 211], [126, 202, 218], [66, 211, 217], [135, 210, 216], [66, 191, 217], [71, 198, 216], [146, 200, 211], [54, 194, 218], [53, 196, 219], [54, 202, 218], [52, 200, 220], [152, 200, 220], [53, 203, 219]], "boundaries": {"contact": [0, 2, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31, 34, 37, 40, 43, 46, 49, 52, 55, 58, 61, 64, 67, 70, 73, 76, 79, 82, 85, 88, 91], "dirichlet": [1, 3, 96, 99, 102, 105, 108, 111, 114, 117, 120, 123, 126, 129, 132, 135, 138, 141, 144, 147]}, "subdomains": {"domain": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389]}}
\ No newline at end of file
diff --git a/docs/examples/ex06.rst b/docs/examples/ex06.rst
index f4188d899..6976cbef0 100644
--- a/docs/examples/ex06.rst
+++ b/docs/examples/ex06.rst
@@ -15,6 +15,14 @@ interpolates any solution vector. The resulting mesh is non-conforming,
i.e. the connectivity between neighboring elements is lost, and hence it can
be used only for visualisation purposes.
+.. note::
+
+ As of 0.4.0, this functionality is included in :func:`~skfem.visuals.matplotlib.plot`,
+ i.e. inputting :class:`~skfem.assembly.InteriorBasis` instead of :class:`~skfem.mesh.Mesh`
+ uses :meth:`~skfem.assembly.InteriorBasis.refinterp` automatically.
+ The steps in this example are still useful when, e.g., exporting to different
+ formats for visualization purposes.
+
As an example, we solve the Poisson equation in a unit square with zero boundary
conditions and biquadratic basis on quadrilateral elements. The quadrilateral
elements are defined using an isoparametric local-to-global mapping.
diff --git a/docs/index.rst b/docs/index.rst
index 67284fa28..ef75f5ce7 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -19,5 +19,5 @@ questions, do not hesitate to drop in and say hello at our `Gitter chat
gettingstarted
tutorial
- learning
+ examples
api
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index b0e1c1b87..79420ba13 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -6,4 +6,4 @@ Tutorials
.. toctree::
tutorials/walkthrough.rst
- examples/ex01.rst
+ tutorials/poisson.rst
diff --git a/docs/examples/ex01.rst b/docs/tutorials/poisson.rst
similarity index 89%
rename from docs/examples/ex01.rst
rename to docs/tutorials/poisson.rst
index d80f6cfa6..02ba586f9 100644
--- a/docs/examples/ex01.rst
+++ b/docs/tutorials/poisson.rst
@@ -39,7 +39,7 @@ times:
After creating the mesh, we evaluate the finite element basis at the global
quadrature points.
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:start-at: ElementTriP1
:end-at: InteriorBasis
@@ -54,13 +54,13 @@ The bilinear and linear forms are defined using the decorators
:func:`~skfem.assembly.linear_form`. It is important to have the order of the
form arguments correct.
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:start-at: bilinear_form
- :end-at: 1.0*v
+ :end-at: 1. * v
All assembly operations are performed using the function :func:`~skfem.assembly.asm`.
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:start-at: asm(laplace, basis)
:end-at: asm(load, basis)
@@ -71,21 +71,21 @@ We are left with solving the assembled linear system.
We eliminate the boundary degrees-of-freedom using
:func:`~skfem.utils.condense` and call :func:`~skfem.utils.solve`.
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:start-at: solve
:end-at: solve
The solution can now be visualised using :meth:`~skfem.visuals.matplotlib.plot`.
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:start-at: main
:end-at: solution.png
-.. figure:: ex01_solution.png
+.. figure:: ../examples/ex01_solution.png
The solution of Poisson equation.
The complete source code reads as follows:
-.. literalinclude:: ex01.py
+.. literalinclude:: ../examples/ex01.py
:linenos:
diff --git a/docs/tutorials/walkthrough.rst b/docs/tutorials/walkthrough.rst
index 38c9fdcf9..8b1ea5da9 100644
--- a/docs/tutorials/walkthrough.rst
+++ b/docs/tutorials/walkthrough.rst
@@ -1,7 +1,7 @@
.. _basic-features:
-Basic features
---------------
+Overview of basic features
+--------------------------
This is an overview of the basic features of scikit-fem that are needed in
most finite element computations.
@@ -136,7 +136,7 @@ See :ref:`learning` for more use cases and instructions.
Setting essential boundary conditions
#####################################
-.. note::
+.. warning::
Using the assembled matrices requires basic understanding of
the finite element method. In particular, to understand
diff --git a/skfem/__init__.py b/skfem/__init__.py
index b37000ec3..7cb6ce084 100644
--- a/skfem/__init__.py
+++ b/skfem/__init__.py
@@ -1,9 +1,12 @@
+"""Support for wildcard import."""
+
from skfem.mesh import *
from skfem.assembly import *
from skfem.mapping import *
from skfem.element import *
from skfem.utils import *
+
__all__ = ['Mesh',
'Mesh2D',
'Mesh3D',
@@ -17,8 +20,10 @@
'InteriorBasis',
'MortarBasis',
'asm',
+ 'Mapping',
'MappingAffine',
'MappingIsoparametric',
+ 'MortarPair',
'adaptive_theta',
'bilinear_form',
'build_pc_ilu',
@@ -44,6 +49,7 @@
'ElementQuad0',
'ElementQuad1',
'ElementQuad2',
+ 'ElementQuadDG',
'ElementTetN0',
'ElementTetP0',
'ElementTetP1',
@@ -56,4 +62,4 @@
'ElementTriRT0',
'ElementVectorH1',
'ElementLineP1',
- 'ElementLineP2',]
+ 'ElementLineP2']
diff --git a/skfem/assembly/basis/basis.py b/skfem/assembly/basis/basis.py
index 82f412993..629c8e4f5 100644
--- a/skfem/assembly/basis/basis.py
+++ b/skfem/assembly/basis/basis.py
@@ -1,3 +1,4 @@
+import warnings
from typing import List, Optional, NamedTuple, Any
import numpy as np
@@ -28,6 +29,19 @@ def __init__(self, mesh, elem, mapping, intorder):
self.mapping = mapping
self._build_dofnum(mesh, elem)
+
+ # human readable names
+ self.dofnames = elem.dofnames
+
+ # global degree-of-freedom location
+ if hasattr(elem, 'doflocs'):
+ doflocs = self.mapping.F(elem.doflocs.T)
+ self.doflocs = np.zeros((doflocs.shape[0], self.N))
+ for itr in range(doflocs.shape[0]):
+ for jtr in range(self.element_dofs.shape[0]):
+ self.doflocs[itr, self.element_dofs[jtr]] =\
+ doflocs[itr, :, jtr]
+
self.mesh = mesh
self.elem = elem
@@ -104,17 +118,8 @@ def _build_dofnum(self, mesh, element):
# interior dofs
self.element_dofs = np.vstack((self.element_dofs, self.interior_dofs))
- # number-of-dofs and human readable names
+ # total dofs
self.N = np.max(self.element_dofs) + 1
- self.dofnames = element.dofnames
-
- if hasattr(element, 'doflocs'):
- doflocs = self.mapping.F(element.doflocs.T)
- self.doflocs = np.zeros((doflocs.shape[0], self.N))
- for itr in range(doflocs.shape[0]):
- for jtr in range(self.element_dofs.shape[0]):
- self.doflocs[itr, self.element_dofs[jtr]] =\
- doflocs[itr, :, jtr]
def complement_dofs(self, *D):
if type(D[0]) is dict:
@@ -188,6 +193,12 @@ def get_dofs(self, facets: Optional[Any] = None):
Mesh.facets_satisfying for each entry to get facets and then find
dofs for those.
+ Returns
+ -------
+ Dofs
+ A subset of degrees-of-freedom as :class:`skfem.assembly.dofs.Dofs`.
+
+
"""
if facets is None:
facets = self.mesh.boundary_facets()
diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py
index 9164e1369..fdc600b31 100644
--- a/skfem/assembly/basis/mortar_basis.py
+++ b/skfem/assembly/basis/mortar_basis.py
@@ -1,65 +1,76 @@
-from typing import Optional
-
import numpy as np
from numpy import ndarray
-from skfem.quadrature import get_quadrature
+from ...quadrature import get_quadrature
from .basis import Basis
from .interior_basis import InteriorBasis
class MortarBasis(Basis):
- """Global basis functions evaluated at integration points on the mortar
- boundary. """
+ """Global basis functions at integration points on the mortar boundary."""
+
def __init__(self,
mesh,
elem,
mapping,
- intorder: Optional[int] = None,
- side: int = 0):
+ intorder):
+ """Initialize a basis for assembling mortar matrices.
+
+ Parameters
+ ----------
+ mesh
+ An object of type :class:`~skfem.mesh.Mesh`.
+ elem
+ An object of type :class:`~skfem.element.Element`.
+ mapping
+ Mapping to the relevant facets of the mesh, see
+ :class:`~skfem.mapping.MortarPair`.
+ intorder
+ Integration order, i.e. the degree of polynomials that are
+ integrated exactly by the used quadrature. Please use equivalent
+ integration orders on both sides of the mortar boundary.
+
+ """
super(MortarBasis, self).__init__(mesh, elem, mapping, intorder)
- self.ib1 = InteriorBasis(mesh.mesh1, elem)
- self.ib2 = InteriorBasis(mesh.mesh2, elem)
-
self.X, self.W = get_quadrature(self.brefdom, self.intorder)
- self.find = np.nonzero(self.mesh.f2t[1, :] != -1)[0]
- self.tind = self.mesh.f2t[side, self.find]
+ # global facets where basis is evaluated
+ self.find = self.mapping.find
+ self.tind = self.mesh.f2t[0, self.find]
# boundary refdom to global facet
- x = self.mapping.G(self.X, find=self.find)
+ x = self.mapping.G(self.X)
+
# global facet to refdom facet
Y = self.mapping.invF(x, tind=self.tind)
- self.normals = np.repeat(mesh.normals[:, :, None], len(self.W), axis=2)
-
- self.nelems = len(self.find)
+ # normals are defined in the mortar mapping
+ self.normals = np.repeat(self.mapping.normals[:, :, None],
+ len(self.W),
+ axis=2)
self.basis = [self.elem.gbasis(self.mapping, Y, j, self.tind)
for j in range(self.Nbfun)]
+ self.nelems = len(self.find)
self.dx = np.abs(self.mapping.detDG(self.X, find=self.find)) *\
np.tile(self.W, (self.nelems, 1))
- if side == 0:
- self.element_dofs = self.ib1.element_dofs[:, self.tind]
- elif side == 1:
- self.element_dofs = self.ib2.element_dofs[:, self.tind - mesh.mesh1.t.shape[1]] + self.ib1.N
-
- self.N = self.ib1.N + self.ib2.N
+ self.element_dofs = self.element_dofs[:, self.tind]
def default_parameters(self):
"""Return default parameters for `~skfem.assembly.asm`."""
- return {'x':self.global_coordinates(),
- 'h':self.mesh_parameters(),
- 'n':self.normals}
-
+ return {'x': self.global_coordinates(),
+ 'h': self.mesh_parameters(),
+ 'n': self.normals}
+
def global_coordinates(self) -> ndarray:
- return self.mapping.G(self.X, find=self.find)
+ return self.mapping.G(self.X)
def mesh_parameters(self) -> ndarray:
if self.mesh.dim() == 1:
return np.array([0.0])
else:
- return np.abs(self.mapping.detDG(self.X, self.find)) ** (1.0 / (self.mesh.dim() - 1))
+ return (np.abs(self.mapping.detDG(self.X, find=self.find))
+ ** (1.0 / (self.mesh.dim() - 1)))
diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py
index 71f76c25c..2c63b6dce 100644
--- a/skfem/element/__init__.py
+++ b/skfem/element/__init__.py
@@ -21,9 +21,9 @@
from .element_tri import ElementTriP1, ElementTriP2, ElementTriDG,\
ElementTriP0, ElementTriRT0, ElementTriMorley,\
ElementTriArgyris
-from .element_quad import ElementQuad0, ElementQuad1, ElementQuad2
+from .element_quad import ElementQuad0, ElementQuad1, ElementQuad2,\
+ ElementQuadDG
from .element_tet import ElementTetP0, ElementTetP1, ElementTetP2,\
ElementTetRT0, ElementTetN0
from .element_hex import ElementHex1
from .element_line import ElementLineP1, ElementLineP2
-
diff --git a/skfem/element/element_hex/element_hex1.py b/skfem/element/element_hex/element_hex1.py
index af2d0f79a..874ef14e6 100644
--- a/skfem/element/element_hex/element_hex1.py
+++ b/skfem/element/element_hex/element_hex1.py
@@ -7,6 +7,14 @@ class ElementHex1(ElementH1):
dim = 3
maxdeg = 3
dofnames = ['u']
+ doflocs = np.array([[-1., -1., -1.],
+ [-1., -1., 1.],
+ [-1., 1., -1.],
+ [ 1., -1., -1.],
+ [-1., 1., 1.],
+ [ 1., -1., 1.],
+ [ 1., 1., -1.],
+ [ 1., 1., 1.]])
def lbasis(self, X, i):
x, y, z = X
diff --git a/skfem/element/element_quad/__init__.py b/skfem/element/element_quad/__init__.py
index e3788daa7..58d9f2de4 100644
--- a/skfem/element/element_quad/__init__.py
+++ b/skfem/element/element_quad/__init__.py
@@ -1,3 +1,4 @@
from .element_quad0 import ElementQuad0
from .element_quad1 import ElementQuad1
from .element_quad2 import ElementQuad2
+from .element_quad_dg import ElementQuadDG
diff --git a/skfem/element/element_quad/element_quad_dg.py b/skfem/element/element_quad/element_quad_dg.py
new file mode 100644
index 000000000..01f1d023d
--- /dev/null
+++ b/skfem/element/element_quad/element_quad_dg.py
@@ -0,0 +1,16 @@
+from ..element_h1 import ElementH1
+
+
+class ElementQuadDG(ElementH1):
+ dim = 2
+
+ def __init__(self, elem):
+ # change all dofs to interior dofs
+ self.elem = elem
+ self.maxdeg = elem.maxdeg
+ self.interior_dofs = (4 * elem.nodal_dofs +
+ 4 * elem.facet_dofs +
+ elem.interior_dofs)
+
+ def lbasis(self, X, i):
+ return self.elem.lbasis(X, i)
diff --git a/skfem/importers/__init__.py b/skfem/importers/__init__.py
index d65017e61..6f8efa2b7 100644
--- a/skfem/importers/__init__.py
+++ b/skfem/importers/__init__.py
@@ -1 +1,3 @@
-from .meshio import from_meshio
+import warnings
+
+from ..io import *
diff --git a/skfem/importers/meshio.py b/skfem/importers/meshio.py
index 032d350d6..37d0bf473 100644
--- a/skfem/importers/meshio.py
+++ b/skfem/importers/meshio.py
@@ -1,120 +1,6 @@
-"""Import any formats supported by meshio."""
-
import warnings
-from collections import OrderedDict
-
-import meshio
-import numpy as np
-
-import skfem
-
-MESH_TYPE_MAPPING = OrderedDict([
- ('tetra', skfem.MeshTet),
- ('hexahedron', skfem.MeshHex),
- ('triangle', skfem.MeshTri),
- ('quad', skfem.MeshQuad),
- ('line', skfem.MeshLine),
-])
-
-
-def from_meshio(m, force_mesh_type=None):
- """Convert meshio mesh into :class:`skfem.mesh.Mesh`.
-
- Parameters
- ----------
- m
- The mesh object from meshio.
- force_mesh_type
- An optional string forcing the mesh type if automatic detection
- fails. See skfem.importers.meshio.MESH_TYPE_MAPPING for possible values.
-
- """
-
- if force_mesh_type is None:
- for k, v in MESH_TYPE_MAPPING.items():
- # find first if match
- if k in m.cells:
- meshio_type, mesh_type = (k, MESH_TYPE_MAPPING[k])
- break
- else:
- meshio_type, mesh_type = (force_mesh_type,
- MESH_TYPE_MAPPING[force_mesh_type])
-
- def strip_extra_coordinates(p):
- if meshio_type == "line":
- return p[:, :1]
- if meshio_type in ("quad", "triangle"):
- return p[:, :2]
- return p
-
- # create p and t
- p = np.ascontiguousarray(strip_extra_coordinates(m.points).T)
- t = np.ascontiguousarray(m.cells[meshio_type].T)
-
- mtmp = mesh_type(p,
- (t[[0, 4, 3, 1, 7, 5, 2, 6]]
- if meshio_type == 'hexahedron'
- # vtk requires a different ordering
- else t))
-
- try:
- # element to boundary element type mapping
- bnd_type = {
- 'line': 'vertex',
- 'triangle': 'line',
- 'quad': 'line',
- 'tetra': 'triangle',
- 'hexahedron': 'quad',
- }[meshio_type]
-
- def find_tagname(tag):
- for key in m.field_data:
- if m.field_data[key][0] == tag:
- return key
- return None
-
- # find subdomains
- if meshio_type in m.cell_data and 'gmsh:physical' in m.cell_data[meshio_type]:
- elements_tag = m.cell_data[meshio_type]['gmsh:physical']
-
- subdomains = {}
- tags = np.unique(elements_tag)
-
- for tag in tags:
- t_set = np.nonzero(tag == elements_tag)[0]
- subdomains[find_tagname(tag)] = t_set
-
- # find tagged boundaries
- if bnd_type in m.cell_data and 'gmsh:physical' in m.cell_data[bnd_type]:
- facets = m.cells[bnd_type]
- facets_tag = m.cell_data[bnd_type]['gmsh:physical']
- bndfacets = mtmp.boundary_facets()
-
- # put meshio facets to dict
- dic = {tuple(np.sort(facets[i])): facets_tag[i]
- for i in range(facets.shape[0])}
-
- # get index of corresponding Mesh.facets for each meshio
- # facet found in the dict
- index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i]
- for i in bndfacets
- if tuple(np.sort(mtmp.facets[:, i])) in dic])
-
- # read meshio tag numbers and names
- tags = index[:, 0]
- boundaries = {}
- for tag in np.unique(tags):
- tagindex = np.nonzero(tags == tag)[0]
- boundaries[find_tagname(tag)] = index[tagindex, 1]
-
- mtmp.boundaries = boundaries
- mtmp.subdomains = subdomains
-
- except Exception as e:
- warnings.warn("Unable to load tagged boundaries/subdomains.")
-
- return mtmp
+warnings.warn("skfem.importers is removed in the next release. "
+ "Use skfem.io instead.", DeprecationWarning)
-def from_file(filename):
- return from_meshio(meshio.read(filename))
+from ..io.meshio import *
diff --git a/skfem/io/__init__.py b/skfem/io/__init__.py
new file mode 100644
index 000000000..d65017e61
--- /dev/null
+++ b/skfem/io/__init__.py
@@ -0,0 +1 @@
+from .meshio import from_meshio
diff --git a/skfem/io/json.py b/skfem/io/json.py
new file mode 100644
index 000000000..0d7fb337e
--- /dev/null
+++ b/skfem/io/json.py
@@ -0,0 +1,34 @@
+"""Import mesh from JSON as defined by :class:`skfem.mesh.to_dict`."""
+
+import json
+
+from skfem.mesh import *
+
+
+def from_file(filename: str):
+ with open(filename, 'r') as handle:
+ d = json.load(handle)
+
+ # detect dimension and number of vertices
+ dim = len(d['p'][0])
+ nverts = len(d['t'][0])
+
+ if dim == 1:
+ mesh_type = MeshLine
+ elif dim == 2 and nverts == 3:
+ mesh_type = MeshTri
+ elif dim == 2 and nverts == 4:
+ mesh_type = MeshQuad
+ elif dim == 3 and nverts == 4:
+ mesh_type = MeshTet
+ elif dim == 3 and nverts == 8:
+ mesh_type = MeshHex
+ else:
+ raise NotImplementedError("The given mesh is not supported.")
+
+ return mesh_type.from_dict(d)
+
+
+def to_file(mesh: Mesh, filename: str):
+ with open(filename, 'w') as handle:
+ d = json.dump(mesh.to_dict(), handle)
diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
new file mode 100644
index 000000000..977a13e5d
--- /dev/null
+++ b/skfem/io/meshio.py
@@ -0,0 +1,121 @@
+"""Import any formats supported by meshio."""
+
+import warnings
+from collections import OrderedDict
+
+import meshio
+import numpy as np
+
+import skfem
+
+
+MESH_TYPE_MAPPING = OrderedDict([
+ ('tetra', skfem.MeshTet),
+ ('hexahedron', skfem.MeshHex),
+ ('triangle', skfem.MeshTri),
+ ('quad', skfem.MeshQuad),
+ ('line', skfem.MeshLine),
+])
+
+
+def from_meshio(m, force_mesh_type=None):
+ """Convert meshio mesh into :class:`skfem.mesh.Mesh`.
+
+ Parameters
+ ----------
+ m
+ The mesh object from meshio.
+ force_mesh_type
+ An optional string forcing the mesh type if automatic detection
+ fails. See skfem.importers.meshio.MESH_TYPE_MAPPING for possible values.
+
+ """
+
+ if force_mesh_type is None:
+ for k, v in MESH_TYPE_MAPPING.items():
+ # find first if match
+ if k in m.cells:
+ meshio_type, mesh_type = (k, MESH_TYPE_MAPPING[k])
+ break
+ else:
+ meshio_type, mesh_type = (force_mesh_type,
+ MESH_TYPE_MAPPING[force_mesh_type])
+
+ def strip_extra_coordinates(p):
+ if meshio_type == "line":
+ return p[:, :1]
+ if meshio_type in ("quad", "triangle"):
+ return p[:, :2]
+ return p
+
+ # create p and t
+ p = np.ascontiguousarray(strip_extra_coordinates(m.points).T)
+ t = np.ascontiguousarray(m.cells[meshio_type].T)
+
+ mtmp = mesh_type(p,
+ (t[[0, 4, 3, 1, 7, 5, 2, 6]]
+ if meshio_type == 'hexahedron'
+ # vtk requires a different ordering
+ else t))
+
+ try:
+ # element to boundary element type mapping
+ bnd_type = {
+ 'line': 'vertex',
+ 'triangle': 'line',
+ 'quad': 'line',
+ 'tetra': 'triangle',
+ 'hexahedron': 'quad',
+ }[meshio_type]
+
+ def find_tagname(tag):
+ for key in m.field_data:
+ if m.field_data[key][0] == tag:
+ return key
+ return None
+
+ # find subdomains
+ if meshio_type in m.cell_data and 'gmsh:physical' in m.cell_data[meshio_type]:
+ elements_tag = m.cell_data[meshio_type]['gmsh:physical']
+
+ subdomains = {}
+ tags = np.unique(elements_tag)
+
+ for tag in tags:
+ t_set = np.nonzero(tag == elements_tag)[0]
+ subdomains[find_tagname(tag)] = t_set
+
+ # find tagged boundaries
+ if bnd_type in m.cell_data and 'gmsh:physical' in m.cell_data[bnd_type]:
+ facets = m.cells[bnd_type]
+ facets_tag = m.cell_data[bnd_type]['gmsh:physical']
+ bndfacets = mtmp.boundary_facets()
+
+ # put meshio facets to dict
+ dic = {tuple(np.sort(facets[i])): facets_tag[i]
+ for i in range(facets.shape[0])}
+
+ # get index of corresponding Mesh.facets for each meshio
+ # facet found in the dict
+ index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i]
+ for i in bndfacets
+ if tuple(np.sort(mtmp.facets[:, i])) in dic])
+
+ # read meshio tag numbers and names
+ tags = index[:, 0]
+ boundaries = {}
+ for tag in np.unique(tags):
+ tagindex = np.nonzero(tags == tag)[0]
+ boundaries[find_tagname(tag)] = index[tagindex, 1]
+
+ mtmp.boundaries = boundaries
+ mtmp.subdomains = subdomains
+
+ except Exception as e:
+ warnings.warn("Unable to load tagged boundaries/subdomains.")
+
+ return mtmp
+
+
+def from_file(filename):
+ return from_meshio(meshio.read(filename))
diff --git a/skfem/mapping/__init__.py b/skfem/mapping/__init__.py
index f40c9e142..c3087c8ea 100644
--- a/skfem/mapping/__init__.py
+++ b/skfem/mapping/__init__.py
@@ -8,3 +8,4 @@
from .mapping import Mapping
from .mapping_affine import MappingAffine
from .mapping_isoparametric import MappingIsoparametric
+from .mortar_pair import MortarPair
diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py
index 542abe436..b3aa5d9f2 100644
--- a/skfem/mapping/mapping_affine.py
+++ b/skfem/mapping/mapping_affine.py
@@ -152,7 +152,12 @@ def G(self, X, find=None):
else:
B, c = self.B[:, :, find], self.c[:, find]
- return (np.einsum('ijk,jl', B, X).T + c.T).T
+ if len(X.shape) == 2:
+ return (np.einsum('ijk,jl', B, X).T + c.T).T
+ elif len(X.shape) == 3:
+ return (np.einsum('ijk,jkl->ikl', B, X).T + c.T).T
+ else:
+ raise Exception("Wrong dimension of input.")
def detDG(self, X, find=None):
if find is None:
@@ -160,7 +165,7 @@ def detDG(self, X, find=None):
else:
detDG = self.detB[find]
- return np.tile(detDG, (X.shape[1], 1)).T
+ return np.tile(detDG, (X.shape[-1], 1)).T
def normals(self, X, tind, find, t2f):
if self.dim == 1:
diff --git a/skfem/mapping/mapping_isoparametric.py b/skfem/mapping/mapping_isoparametric.py
index 4102eab08..690fc58d8 100644
--- a/skfem/mapping/mapping_isoparametric.py
+++ b/skfem/mapping/mapping_isoparametric.py
@@ -81,7 +81,7 @@ def bndmap(i, X, find=None):
out += p[i, facets[itr, :]][:, None] * phi
return out
else:
- out = np.zeros((len(find), X.shape[1]))
+ out = np.zeros((len(find), X.shape[-1]))
for itr in range(facets.shape[0]):
phi, _ = bndelem.lbasis(X, itr)
out += p[i, facets[itr, find]][:, None] * phi
@@ -95,7 +95,7 @@ def bndJ(i, j, X, find=None):
out += p[i, facets[itr, :]][:, None] * dphi[j]
return out
else:
- out = np.zeros((len(find), X.shape[1]))
+ out = np.zeros((len(find), X.shape[-1]))
for itr in range(facets.shape[0]):
_, dphi = bndelem.lbasis(X, itr)
out += p[i, facets[itr, find]][:, None] * dphi[j]
diff --git a/skfem/mapping/mortar_pair.py b/skfem/mapping/mortar_pair.py
new file mode 100644
index 000000000..1c6633fd0
--- /dev/null
+++ b/skfem/mapping/mortar_pair.py
@@ -0,0 +1,163 @@
+from typing import NamedTuple
+
+import numpy as np
+from numpy import ndarray
+
+from ..mesh.mesh2d import Mesh2D
+from .mapping import Mapping
+
+
+class MortarPair(NamedTuple):
+ """A pair of mappings for mortar methods.
+
+ In mortar methods, we are enforcing interface conditions on non-matching
+ meshes. To this end, we must generate matching quadrature points on both
+ sides of the interface.
+
+ :class:`~skfem.mapping.MortarPair` is a pair of
+ :class:`~skfem.mapping.Mapping` objects defined so that mapping any local
+ point in the reference element using both of the mappings, we get matching
+ global points on the original non-matching meshes.
+
+ Attributes
+ ----------
+ mapping1
+ Mapping to the facets of the first mesh.
+ mapping2
+ Mapping to the facets of the second mesh.
+
+ """
+
+ mapping1: Mapping = None
+ mapping2: Mapping = None
+
+ @classmethod
+ def init_2D(cls,
+ mesh1: Mesh2D,
+ mesh2: Mesh2D,
+ boundary1: ndarray,
+ boundary2: ndarray,
+ tangent: ndarray = None):
+ """Create mortar mappings for two 2D meshes via projection.
+
+ Parameters
+ ----------
+ mesh1
+ An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`.
+ mesh2
+ An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`.
+ boundary1
+ A subset of facets to use from mesh1.
+ boundary2
+ A subset of facets to use from mesh2.
+ tangent
+ A tangent vector defining the direction of the projection.
+ If not given, use first facet of boundary1 as the vector.
+
+ """
+ from ..mesh import MeshLine
+
+ if tangent is None:
+ tangent = (mesh1.p[:, mesh1.facets[0, boundary1[0]]] -
+ mesh1.p[:, mesh1.facets[1, boundary1[0]]])
+
+ tangent /= np.linalg.norm(tangent)
+
+ # find unique nodes on the two boundaries
+ p1_ix = np.unique(mesh1.facets[:, boundary1].flatten())
+ p2_ix = np.unique(mesh2.facets[:, boundary2].flatten())
+ p1 = mesh1.p[:, p1_ix]
+ p2 = mesh2.p[:, p2_ix]
+
+ def proj(p):
+ """Project onto the line defined by 'tangent'."""
+ return np.outer(tangent, tangent) @ p
+
+ def param(p):
+ """Calculate signed distances of projected points from origin."""
+ y = proj(p)
+ return np.linalg.norm(y, axis=0) * np.sign(np.dot(tangent, y))
+
+ # find unique supermesh facets by combining nodes from both sides
+ param_p1 = param(p1)
+ param_p2 = param(p2)
+ _, ix = np.unique(np.concatenate((param_p1, param_p2)), return_index=True)
+ ixorig = np.concatenate((p1_ix, p2_ix + mesh1.p.shape[1]))[ix]
+ p = np.array([np.hstack((param(mesh1.p), param(mesh2.p)))])
+ t = np.array([ixorig[:-1], ixorig[1:]])
+
+ # create 1-dimensional supermesh from the intersections of the projected
+ # facet elements
+ p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))]
+ range_max = np.min([np.max(param_p1), np.max(param_p2)])
+ range_min = np.max([np.min(param_p1), np.min(param_p2)])
+ p = np.array([p[0, (p[0] <= range_max) * (p[0] >= range_min)]])
+ t = np.array([np.arange(p.shape[1] - 1), np.arange(1, p.shape[1])])
+ m_super = MeshLine(p, t)
+
+ # helper meshes for creating the mappings
+ m1 = MeshLine(np.sort(param_p1), np.array([np.arange(p1.shape[1] - 1),
+ np.arange(1, p1.shape[1])]))
+ m2 = MeshLine(np.sort(param_p2), np.array([np.arange(p2.shape[1] - 1),
+ np.arange(1, p2.shape[1])]))
+
+ # construct normals by rotating 'tangent'
+ normal = np.array([tangent[1], -tangent[0]])
+ normals = normal[:, None].repeat(t.shape[1], axis=1)
+
+ # initialize mappings (for orienting)
+ map_super = m_super.mapping()
+ map_m1 = m1.mapping()
+ map_m2 = m2.mapping()
+ map_mesh1 = mesh1.mapping()
+ map_mesh2 = mesh2.mapping()
+
+ # matching of elements in the supermesh and the helper meshes
+ mps = map_super.F(np.array([[0.5]]))
+ ix1 = np.digitize(mps[0, :, 0], m1.p[0]) - 1
+ ix2 = np.digitize(mps[0, :, 0], m2.p[0]) - 1
+
+ # for each element, map two points to global coordinates, reparametrize
+ # the points, and flip corresponding helper mesh element indices if
+ # sorting is wrong
+ f1mps = .5 * (mesh1.p[:, mesh1.facets[0, boundary1]] +
+ mesh1.p[:, mesh1.facets[1, boundary1]])
+ sort_boundary1 = np.argsort(param(f1mps))
+ z1 = map_mesh1.G(map_m1.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix1),
+ find=boundary1[sort_boundary1][ix1])
+ ix1_flip = np.unique(ix1[param(z1[:, :, 1]) < param(z1[:, :, 0])])
+ m1.t[:, ix1_flip] = np.flipud(m1.t[:, ix1_flip])
+
+ f2mps = .5 * (mesh2.p[:, mesh2.facets[0, boundary2]] +
+ mesh2.p[:, mesh2.facets[1, boundary2]])
+ sort_boundary2 = np.argsort(param(f2mps))
+ z2 = map_mesh2.G(map_m2.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix2),
+ find=boundary2[sort_boundary2][ix2])
+ ix2_flip = np.unique(ix2[param(z2[:, :, 1]) < param(z2[:, :, 0])])
+ m2.t[:, ix2_flip] = np.flipud(m2.t[:, ix2_flip])
+
+ # create new mapping objects with G replaced by the supermesh mapping
+ new_map1 = mesh1.mapping()
+ new_map2 = mesh2.mapping()
+ map_m1 = m1.mapping()
+ map_m2 = m2.mapping()
+ new_map1.G = lambda X: map_mesh1.G(map_m1.invF(map_super.F(X), tind=ix1),
+ find=boundary1[sort_boundary1][ix1])
+ new_map2.G = lambda X: map_mesh2.G(map_m2.invF(map_super.F(X), tind=ix2),
+ find=boundary2[sort_boundary2][ix2])
+
+ # these are used by :class:`skfem.assembly.basis.MortarBasis`
+ new_map1.find = boundary1[sort_boundary1][ix1]
+ new_map2.find = boundary2[sort_boundary2][ix2]
+ new_map1.normals = normals
+ new_map2.normals = normals
+
+ # method for calculating the lengths of the segments ('detDG')
+ segments1 = new_map1.G(np.array([[0.0, 1.0]]))
+ segments2 = new_map2.G(np.array([[0.0, 1.0]]))
+ new_map1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 +
+ np.diff(segments1[1], axis=1)**2)
+ new_map2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 +
+ np.diff(segments2[1], axis=1)**2)
+
+ return cls(mapping1=new_map1, mapping2=new_map2)
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index 00a00d36a..3db862b4c 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -312,7 +312,9 @@ def element_finder(self) -> Callable[[ndarray], ndarray]:
raise NotImplementedError("element_finder not implemented "
"for the given Mesh type.")
- def nodes_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
+ def nodes_satisfying(self,
+ test: Callable[[ndarray], bool],
+ boundaries_only: bool = False) -> ndarray:
"""Return nodes that satisfy some condition.
Parameters
@@ -320,11 +322,18 @@ def nodes_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
test
A function which returns True for the set of nodes that are to be
included in the return set.
+ boundaries_only
+ If True, include only boundary facets.
"""
- return np.nonzero(test(self.p))[0]
-
- def facets_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
+ nodes = np.nonzero(test(self.p))[0]
+ if boundaries_only:
+ nodes = np.intersect1d(nodes, self.boundary_nodes())
+ return nodes
+
+ def facets_satisfying(self,
+ test: Callable[[ndarray], bool],
+ boundaries_only: bool = False) -> ndarray:
"""Return facets whose midpoints satisfy some condition.
Parameters
@@ -332,11 +341,16 @@ def facets_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
test
A function which returns True for the facet midpoints that are to
be included in the return set.
+ boundaries_only
+ If True, include only boundary facets.
"""
midp = [np.sum(self.p[itr, self.facets], axis=0) / self.facets.shape[0]
for itr in range(self.p.shape[0])]
- return np.nonzero(test(np.array(midp)))[0]
+ facets = np.nonzero(test(np.array(midp)))[0]
+ if boundaries_only:
+ facets = np.intersect1d(facets, self.boundary_facets())
+ return facets
def elements_satisfying(self,
test: Callable[[ndarray], bool]) -> ndarray:
@@ -357,7 +371,7 @@ def elements_satisfying(self,
def from_dict(cls: Type[MeshType], d) -> MeshType:
"""Initialize a mesh from a dictionary."""
if 'p' not in d or 't' not in d:
- raise ValueError("Dictionary must contain keys 'p' and 't'")
+ raise ValueError("Dictionary must contain keys 'p' and 't'.")
else:
d['p'] = np.array(d['p']).T
d['t'] = np.array(d['t']).T
diff --git a/tests/test_assembly.py b/tests/test_assembly.py
index 544d9b3fb..49b525522 100644
--- a/tests/test_assembly.py
+++ b/tests/test_assembly.py
@@ -1,4 +1,3 @@
-"""Unit tests for assembly operations."""
import unittest
import numpy as np
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index 39ff41683..75c716db1 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -1,6 +1,3 @@
-# -*- coding: utf-8 -*-
-"""Check the convergence rates of the elements."""
-
import unittest
import numpy as np
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 9f7f3c4fc..63c902059 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -1,84 +1,70 @@
-"""These tests run the examples and check that their output stays constant."""
import unittest
import numpy as np
class TestEx01(unittest.TestCase):
- """Run examples/ex01.py"""
-
def runTest(self):
import docs.examples.ex01 as ex01
self.assertAlmostEqual(np.max(ex01.x), 0.07344576657)
class TestEx02(unittest.TestCase):
- """Run examples/ex02.py"""
-
def runTest(self):
import docs.examples.ex02 as ex02
self.assertAlmostEqual(np.max(ex02.x), 0.001217973811129439)
class TestEx03(unittest.TestCase):
- """Run examples/ex03.py"""
-
def runTest(self):
import docs.examples.ex03 as ex03
self.assertAlmostEqual(ex03.L[0], 0.00418289)
-class TestEx05(unittest.TestCase):
- """Run examples/ex05.py"""
+class TestEx04(unittest.TestCase):
+ def runTest(self):
+ import docs.examples.ex04 as ex04
+ self.assertAlmostEqual(np.max(ex04.vonmises1), 64.57919892367978)
+ self.assertAlmostEqual(np.max(ex04.vonmises2), 67.91419753783893)
+
+class TestEx05(unittest.TestCase):
def runTest(self):
import docs.examples.ex05 as ex05
self.assertAlmostEqual(np.max(ex05.x), 0.93570751751091152)
class TestEx06(unittest.TestCase):
- """Run examples/ex06.py"""
-
def runTest(self):
import docs.examples.ex06 as ex06
self.assertAlmostEqual(np.max(ex06.x), 0.073651530833125131)
class TestEx07(unittest.TestCase):
- """Run examples/ex07.py"""
-
def runTest(self):
import docs.examples.ex07 as ex07
self.assertAlmostEqual(np.max(ex07.x), 0.07869083767545548)
class TestEx08(unittest.TestCase):
- """Run examples/ex08.py"""
-
def runTest(self):
- pass
+ import docs.examples.ex08 as ex08
# only run the initialization, nothing to test
class TestEx09(unittest.TestCase):
- """Run examples/ex09.py"""
-
def runTest(self):
import docs.examples.ex09 as ex09
self.assertAlmostEqual(np.max(ex09.x), 0.055596791644282988)
class TestEx10(unittest.TestCase):
- """Run examples/ex10.py"""
-
def runTest(self):
import docs.examples.ex10 as ex10
self.assertAlmostEqual(np.mean(ex10.x), 0.277931521728906)
class TestEx11(unittest.TestCase):
- """Run examples/ex11.py"""
-
def runTest(self):
import docs.examples.ex11 as ex11
u = ex11.u
@@ -105,8 +91,6 @@ def runTest(self):
class TestEx14(unittest.TestCase):
- """Run examples/ex14.py"""
-
def runTest(self):
import docs.examples.ex14
u = docs.examples.ex14.u
@@ -115,35 +99,35 @@ def runTest(self):
class TestEx15(unittest.TestCase):
- """Run examples/ex15.py"""
-
def runTest(self):
- import docs.examples.ex15
- self.assertTrue(np.max(docs.examples.ex15.x) - 0.1234567 < 1e-5)
+ import docs.examples.ex15 as ex15
+ self.assertTrue(np.max(ex15.x) - 0.1234567 < 1e-5)
class TestEx16(unittest.TestCase):
- """Run examples/ex16.py"""
-
def runTest(self):
- import docs.examples.ex16
- self.assertTrue(np.linalg.norm(np.array([0, 2, 6, 12, 20, 30]) - docs.examples.ex16.ks) < 0.4)
- self.assertTrue(docs.examples.ex16.ks[-1], 30.309720458315521)
+ import docs.examples.ex16 as ex16
+ self.assertTrue(np.linalg.norm(np.array([0, 2, 6, 12, 20, 30])
+ - ex16.ks) < 0.4)
+ self.assertTrue(ex16.ks[-1], 30.309720458315521)
class TestEx17(unittest.TestCase):
def runTest(self):
import docs.examples.ex17 as ex
+ # TODO improve
class TestEx18(unittest.TestCase):
def runTest(self):
import docs.examples.ex18 as ex
+ # TODO improve
class TestEx19(unittest.TestCase):
def runTest(self):
import docs.examples.ex19 as ex
+ # TODO improve
class TestEx20(unittest.TestCase):
@@ -177,22 +161,23 @@ def runTest(self):
self.assertAlmostEqual(max(ex.lmbda_list), ex.turning_point,
delta=5e-5)
+
class TestEx24(unittest.TestCase):
def runTest(self):
- pass
+ import docs.examples.ex24 as ex24
class TestEx25(unittest.TestCase):
def runTest(self):
- import docs.examples.ex25 as ex
- mu = np.mean(ex.t)
+ import docs.examples.ex25 as ex25
+ mu = np.mean(ex25.t)
self.assertAlmostEqual(mu, 0.4642600944590631, places=5)
- self.assertAlmostEqual(np.mean(ex.t0), mu, places=2)
+ self.assertAlmostEqual(np.mean(ex25.t0), mu, places=2)
class TestEx26(unittest.TestCase):
def runTest(self):
- pass
+ import docs.examples.ex26 as ex26
class TestEx27(unittest.TestCase):
diff --git a/tests/test_manufactured.py b/tests/test_manufactured.py
index d9c6e0cd5..6fe667b2b 100644
--- a/tests/test_manufactured.py
+++ b/tests/test_manufactured.py
@@ -1,4 +1,5 @@
"""Solve problems that have manufactured solutions."""
+
import unittest
import numpy as np
diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index 44f3d1035..6b2d807ce 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -73,5 +73,54 @@ def initialize_meshes(self):
return m
+class TestMortarPair(unittest.TestCase):
+ """Check that mapped points match."""
+
+ mesh1_type = MeshTri
+ mesh2_type = MeshTri
+ nrefs1 = 2
+ nrefs2 = 3
+ translate_y = 0.0
+
+ def init_meshes(self):
+ m1 = self.mesh1_type()
+ m1.refine(self.nrefs1)
+ m2 = self.mesh2_type()
+ m2.refine(self.nrefs2)
+ m2.translate((1.0, self.translate_y))
+ return m1, m2
+
+ def runTest(self):
+ m1, m2 = self.init_meshes()
+ mp = MortarPair.init_2D(m1, m2,
+ m1.facets_satisfying(lambda x: x[0] == 1.),
+ m2.facets_satisfying(lambda x: x[0] == 1.))
+ test_points = np.array([np.linspace(0., 1., 7)])
+ self.assertTrue((mp.mapping1.G(test_points) -
+ mp.mapping1.G(test_points) < 1e-10).all())
+
+
+class TestMortarPairTriQuad(TestMortarPair):
+ mesh1_type = MeshTri
+ mesh2_type = MeshQuad
+
+
+class TestMortarPairQuadQuad(TestMortarPair):
+ mesh1_type = MeshQuad
+ mesh2_type = MeshQuad
+
+
+class TestMortarPairNoMatch1(TestMortarPair):
+ mesh1_type = MeshQuad
+ mesh2_type = MeshTri
+ translate_y = 0.1
+
+
+class TestMortarPairNoMatch2(TestMortarPair):
+ mesh1_type = MeshQuad
+ mesh2_type = MeshTri
+ translate_y = -np.pi / 10.
+
+
if __name__ == '__main__':
unittest.main()
diff --git a/tests/test_utils.py b/tests/test_utils.py
index afab46d2f..d482bddc3 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,4 +1,3 @@
-"""Unit tests for utils package."""
import unittest
import numpy as np
diff --git a/tests/test_visuals.py b/tests/test_visuals.py
index e36868f64..cebd1873f 100644
--- a/tests/test_visuals.py
+++ b/tests/test_visuals.py
@@ -1,6 +1,5 @@
import unittest
-
from skfem.mesh import *
from skfem.visuals.matplotlib import *