No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Jan 24, 06:18

View Options

	diff --git a/docs/details/figures/MeshQuad4/FineLayer/element-numbers.py b/docs/details/figures/MeshQuad4/FineLayer/element-numbers.py
	index 2e93989..aa8a5ec 100644
	--- a/docs/details/figures/MeshQuad4/FineLayer/element-numbers.py
	+++ b/docs/details/figures/MeshQuad4/FineLayer/element-numbers.py
	@@ -1,33 +1,33 @@

	import GooseFEM as gf
	import matplotlib.pyplot as plt
	import GooseMPL as gplt
	import numpy as np

	plt.style.use(['goose'])

	# --------------------------------------------------------------------------------------------------

	mesh = gf.Mesh.Quad4.FineLayer(6, 18)
	coor = mesh.coor()
	conn = mesh.conn()

	cindex = np.arange(mesh.nelem())

	# --------------------------------------------------------------------------------------------------

	fig, ax = plt.subplots(figsize=(10,10))

	im = gplt.patch(coor=coor, conn=conn, cindex=cindex, cmap='jet')

	for e in range(conn.shape[0]):
	x = np.mean(coor[conn[e,:], 0])
	y = np.mean(coor[conn[e,:], 1])
	ax.text(x, y, str(e), ha='center', va='center')

	ax.set_aspect('equal')
	ax.get_xaxis().set_visible(False)
	ax.get_yaxis().set_visible(False)

	-plt.savefig('element_numbers.svg')
	+plt.savefig('element-numbers.svg')
	plt.close()
	diff --git a/docs/details/figures/MeshQuad4/FineLayer/element_numbers.svg b/docs/details/figures/MeshQuad4/FineLayer/element-numbers.svg
	similarity index 100%
	rename from docs/details/figures/MeshQuad4/FineLayer/element_numbers.svg
	rename to docs/details/figures/MeshQuad4/FineLayer/element-numbers.svg
	diff --git a/docs/details/figures/MeshQuad4/FineLayer/elementgrid.py b/docs/details/figures/MeshQuad4/FineLayer/elementgrid.py
	new file mode 100644
	index 0000000..97570df
	--- /dev/null
	+++ b/docs/details/figures/MeshQuad4/FineLayer/elementgrid.py
	@@ -0,0 +1,27 @@
	+
	+import GooseFEM as gf
	+import matplotlib.pyplot as plt
	+import GooseMPL as gplt
	+import numpy as np
	+
	+plt.style.use(['goose'])
	+
	+mesh = gf.Mesh.Quad4.FineLayer(12, 18)
	+coor = mesh.coor()
	+conn = mesh.conn()
	+layer = mesh.elementsMiddleLayer()
	+select = mesh.elementgrid_around_ravel(layer[7], 2)
	+
	+I = np.zeros((conn.shape[0]))
	+I[select] = 1
	+
	+fig, ax = plt.subplots()
	+
	+im = gplt.patch(coor=coor, conn=conn, cindex=I, cmap='Reds', clim=[0, 1])
	+
	+ax.set_aspect('equal')
	+ax.get_xaxis().set_visible(False)
	+ax.get_yaxis().set_visible(False)
	+
	+fig.savefig('elementgrid.svg')
	+plt.close(fig)
	diff --git a/docs/details/figures/MeshQuad4/FineLayer/elementgrid.svg b/docs/details/figures/MeshQuad4/FineLayer/elementgrid.svg
	new file mode 100644
	index 0000000..9770059
	--- /dev/null
	+++ b/docs/details/figures/MeshQuad4/FineLayer/elementgrid.svg
	@@ -0,0 +1,718 @@
	+<?xml version="1.0" encoding="utf-8" standalone="no"?>
	+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
	+ "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
	+<!-- Created with matplotlib (https://matplotlib.org/) -->
	+<svg height="432pt" version="1.1" viewBox="0 0 576 432" width="576pt" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
	+ <metadata>
	+ <rdf:RDF xmlns:cc="http://creativecommons.org/ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
	+ <cc:Work>
	+ <dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/>
	+ <dc:date>2020-12-18T18:27:41.433494</dc:date>
	+ <dc:format>image/svg+xml</dc:format>
	+ <dc:creator>
	+ <cc:Agent>
	+ <dc:title>Matplotlib v3.3.3, https://matplotlib.org/</dc:title>
	+ </cc:Agent>
	+ </dc:creator>
	+ </cc:Work>
	+ </rdf:RDF>
	+ </metadata>
	+ <defs>
	+ <style type="text/css">*{stroke-linecap:butt;stroke-linejoin:round;}</style>
	+ </defs>
	+ <g id="figure_1">
	+ <g id="patch_1">
	+ <path d="M 0 432
	+L 576 432
	+L 576 0
	+L 0 0
	+z
	+" style="fill:none;"/>
	+ </g>
	+ <g id="axes_1">
	+ <g id="patch_2">
	+ <path d="M 176.445714 411.22
	+L 399.554286 411.22
	+L 399.554286 20.78
	+L 176.445714 20.78
	+z
	+" style="fill:none;"/>
	+ </g>
	+ <g id="PatchCollection_1">
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 378.683333
	+L 241.519048 378.683333
	+L 241.519048 332.202381
	+L 195.038095 332.202381
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 378.683333
	+L 288 378.683333
	+L 288 332.202381
	+L 241.519048 332.202381
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 378.683333
	+L 334.480952 378.683333
	+L 334.480952 332.202381
	+L 288 332.202381
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 378.683333
	+L 380.961905 378.683333
	+L 380.961905 332.202381
	+L 334.480952 332.202381
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 332.202381
	+L 241.519048 332.202381
	+L 241.519048 285.721429
	+L 195.038095 285.721429
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 332.202381
	+L 288 332.202381
	+L 288 285.721429
	+L 241.519048 285.721429
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 332.202381
	+L 334.480952 332.202381
	+L 334.480952 285.721429
	+L 288 285.721429
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 332.202381
	+L 380.961905 332.202381
	+L 380.961905 285.721429
	+L 334.480952 285.721429
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 285.721429
	+L 241.519048 285.721429
	+L 226.025397 270.227778
	+L 210.531746 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 285.721429
	+L 241.519048 254.734127
	+L 226.025397 254.734127
	+L 226.025397 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 270.227778
	+L 226.025397 270.227778
	+L 226.025397 254.734127
	+L 210.531746 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 285.721429
	+L 210.531746 270.227778
	+L 210.531746 254.734127
	+L 195.038095 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 285.721429
	+L 288 285.721429
	+L 272.506349 270.227778
	+L 257.012698 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 285.721429
	+L 288 254.734127
	+L 272.506349 254.734127
	+L 272.506349 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 270.227778
	+L 272.506349 270.227778
	+L 272.506349 254.734127
	+L 257.012698 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 285.721429
	+L 257.012698 270.227778
	+L 257.012698 254.734127
	+L 241.519048 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 285.721429
	+L 334.480952 285.721429
	+L 318.987302 270.227778
	+L 303.493651 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 285.721429
	+L 334.480952 254.734127
	+L 318.987302 254.734127
	+L 318.987302 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 270.227778
	+L 318.987302 270.227778
	+L 318.987302 254.734127
	+L 303.493651 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 285.721429
	+L 303.493651 270.227778
	+L 303.493651 254.734127
	+L 288 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 285.721429
	+L 380.961905 285.721429
	+L 365.468254 270.227778
	+L 349.974603 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 380.961905 285.721429
	+L 380.961905 254.734127
	+L 365.468254 254.734127
	+L 365.468254 270.227778
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 270.227778
	+L 365.468254 270.227778
	+L 365.468254 254.734127
	+L 349.974603 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 285.721429
	+L 349.974603 270.227778
	+L 349.974603 254.734127
	+L 334.480952 254.734127
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 254.734127
	+L 210.531746 254.734127
	+L 210.531746 239.240476
	+L 195.038095 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 254.734127
	+L 226.025397 254.734127
	+L 226.025397 239.240476
	+L 210.531746 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 254.734127
	+L 241.519048 254.734127
	+L 241.519048 239.240476
	+L 226.025397 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 254.734127
	+L 257.012698 254.734127
	+L 257.012698 239.240476
	+L 241.519048 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 254.734127
	+L 272.506349 254.734127
	+L 272.506349 239.240476
	+L 257.012698 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 254.734127
	+L 288 254.734127
	+L 288 239.240476
	+L 272.506349 239.240476
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 254.734127
	+L 303.493651 254.734127
	+L 303.493651 239.240476
	+L 288 239.240476
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 254.734127
	+L 318.987302 254.734127
	+L 318.987302 239.240476
	+L 303.493651 239.240476
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 254.734127
	+L 334.480952 254.734127
	+L 334.480952 239.240476
	+L 318.987302 239.240476
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 254.734127
	+L 349.974603 254.734127
	+L 349.974603 239.240476
	+L 334.480952 239.240476
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 254.734127
	+L 365.468254 254.734127
	+L 365.468254 239.240476
	+L 349.974603 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 254.734127
	+L 380.961905 254.734127
	+L 380.961905 239.240476
	+L 365.468254 239.240476
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 239.240476
	+L 210.531746 239.240476
	+L 210.531746 223.746825
	+L 195.038095 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 239.240476
	+L 226.025397 239.240476
	+L 226.025397 223.746825
	+L 210.531746 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 239.240476
	+L 241.519048 239.240476
	+L 241.519048 223.746825
	+L 226.025397 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 239.240476
	+L 257.012698 239.240476
	+L 257.012698 223.746825
	+L 241.519048 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 239.240476
	+L 272.506349 239.240476
	+L 272.506349 223.746825
	+L 257.012698 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 239.240476
	+L 288 239.240476
	+L 288 223.746825
	+L 272.506349 223.746825
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 239.240476
	+L 303.493651 239.240476
	+L 303.493651 223.746825
	+L 288 223.746825
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 239.240476
	+L 318.987302 239.240476
	+L 318.987302 223.746825
	+L 303.493651 223.746825
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 239.240476
	+L 334.480952 239.240476
	+L 334.480952 223.746825
	+L 318.987302 223.746825
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 239.240476
	+L 349.974603 239.240476
	+L 349.974603 223.746825
	+L 334.480952 223.746825
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 239.240476
	+L 365.468254 239.240476
	+L 365.468254 223.746825
	+L 349.974603 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 239.240476
	+L 380.961905 239.240476
	+L 380.961905 223.746825
	+L 365.468254 223.746825
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 223.746825
	+L 210.531746 223.746825
	+L 210.531746 208.253175
	+L 195.038095 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 223.746825
	+L 226.025397 223.746825
	+L 226.025397 208.253175
	+L 210.531746 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 223.746825
	+L 241.519048 223.746825
	+L 241.519048 208.253175
	+L 226.025397 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 223.746825
	+L 257.012698 223.746825
	+L 257.012698 208.253175
	+L 241.519048 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 223.746825
	+L 272.506349 223.746825
	+L 272.506349 208.253175
	+L 257.012698 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 223.746825
	+L 288 223.746825
	+L 288 208.253175
	+L 272.506349 208.253175
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 223.746825
	+L 303.493651 223.746825
	+L 303.493651 208.253175
	+L 288 208.253175
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 223.746825
	+L 318.987302 223.746825
	+L 318.987302 208.253175
	+L 303.493651 208.253175
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 223.746825
	+L 334.480952 223.746825
	+L 334.480952 208.253175
	+L 318.987302 208.253175
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 223.746825
	+L 349.974603 223.746825
	+L 349.974603 208.253175
	+L 334.480952 208.253175
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 223.746825
	+L 365.468254 223.746825
	+L 365.468254 208.253175
	+L 349.974603 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 223.746825
	+L 380.961905 223.746825
	+L 380.961905 208.253175
	+L 365.468254 208.253175
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 208.253175
	+L 210.531746 208.253175
	+L 210.531746 192.759524
	+L 195.038095 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 208.253175
	+L 226.025397 208.253175
	+L 226.025397 192.759524
	+L 210.531746 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 208.253175
	+L 241.519048 208.253175
	+L 241.519048 192.759524
	+L 226.025397 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 208.253175
	+L 257.012698 208.253175
	+L 257.012698 192.759524
	+L 241.519048 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 208.253175
	+L 272.506349 208.253175
	+L 272.506349 192.759524
	+L 257.012698 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 208.253175
	+L 288 208.253175
	+L 288 192.759524
	+L 272.506349 192.759524
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 208.253175
	+L 303.493651 208.253175
	+L 303.493651 192.759524
	+L 288 192.759524
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 208.253175
	+L 318.987302 208.253175
	+L 318.987302 192.759524
	+L 303.493651 192.759524
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 208.253175
	+L 334.480952 208.253175
	+L 334.480952 192.759524
	+L 318.987302 192.759524
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 208.253175
	+L 349.974603 208.253175
	+L 349.974603 192.759524
	+L 334.480952 192.759524
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 208.253175
	+L 365.468254 208.253175
	+L 365.468254 192.759524
	+L 349.974603 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 208.253175
	+L 380.961905 208.253175
	+L 380.961905 192.759524
	+L 365.468254 192.759524
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 192.759524
	+L 210.531746 192.759524
	+L 210.531746 177.265873
	+L 195.038095 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 192.759524
	+L 226.025397 192.759524
	+L 226.025397 177.265873
	+L 210.531746 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 192.759524
	+L 241.519048 192.759524
	+L 241.519048 177.265873
	+L 226.025397 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 192.759524
	+L 257.012698 192.759524
	+L 257.012698 177.265873
	+L 241.519048 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 192.759524
	+L 272.506349 192.759524
	+L 272.506349 177.265873
	+L 257.012698 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 192.759524
	+L 288 192.759524
	+L 288 177.265873
	+L 272.506349 177.265873
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 192.759524
	+L 303.493651 192.759524
	+L 303.493651 177.265873
	+L 288 177.265873
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 192.759524
	+L 318.987302 192.759524
	+L 318.987302 177.265873
	+L 303.493651 177.265873
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 192.759524
	+L 334.480952 192.759524
	+L 334.480952 177.265873
	+L 318.987302 177.265873
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 192.759524
	+L 349.974603 192.759524
	+L 349.974603 177.265873
	+L 334.480952 177.265873
	+z
	+" style="fill:#67000d;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 192.759524
	+L 365.468254 192.759524
	+L 365.468254 177.265873
	+L 349.974603 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 192.759524
	+L 380.961905 192.759524
	+L 380.961905 177.265873
	+L 365.468254 177.265873
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 177.265873
	+L 210.531746 177.265873
	+L 210.531746 161.772222
	+L 195.038095 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 177.265873
	+L 226.025397 177.265873
	+L 226.025397 161.772222
	+L 210.531746 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 226.025397 177.265873
	+L 241.519048 177.265873
	+L 241.519048 146.278571
	+L 226.025397 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 210.531746 161.772222
	+L 226.025397 161.772222
	+L 241.519048 146.278571
	+L 195.038095 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 177.265873
	+L 257.012698 177.265873
	+L 257.012698 161.772222
	+L 241.519048 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 177.265873
	+L 272.506349 177.265873
	+L 272.506349 161.772222
	+L 257.012698 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 272.506349 177.265873
	+L 288 177.265873
	+L 288 146.278571
	+L 272.506349 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 257.012698 161.772222
	+L 272.506349 161.772222
	+L 288 146.278571
	+L 241.519048 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 177.265873
	+L 303.493651 177.265873
	+L 303.493651 161.772222
	+L 288 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 177.265873
	+L 318.987302 177.265873
	+L 318.987302 161.772222
	+L 303.493651 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 318.987302 177.265873
	+L 334.480952 177.265873
	+L 334.480952 146.278571
	+L 318.987302 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 303.493651 161.772222
	+L 318.987302 161.772222
	+L 334.480952 146.278571
	+L 288 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 177.265873
	+L 349.974603 177.265873
	+L 349.974603 161.772222
	+L 334.480952 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 177.265873
	+L 365.468254 177.265873
	+L 365.468254 161.772222
	+L 349.974603 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 365.468254 177.265873
	+L 380.961905 177.265873
	+L 380.961905 146.278571
	+L 365.468254 161.772222
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 349.974603 161.772222
	+L 365.468254 161.772222
	+L 380.961905 146.278571
	+L 334.480952 146.278571
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 146.278571
	+L 241.519048 146.278571
	+L 241.519048 99.797619
	+L 195.038095 99.797619
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 146.278571
	+L 288 146.278571
	+L 288 99.797619
	+L 241.519048 99.797619
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 146.278571
	+L 334.480952 146.278571
	+L 334.480952 99.797619
	+L 288 99.797619
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 146.278571
	+L 380.961905 146.278571
	+L 380.961905 99.797619
	+L 334.480952 99.797619
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 195.038095 99.797619
	+L 241.519048 99.797619
	+L 241.519048 53.316667
	+L 195.038095 53.316667
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 241.519048 99.797619
	+L 288 99.797619
	+L 288 53.316667
	+L 241.519048 53.316667
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 288 99.797619
	+L 334.480952 99.797619
	+L 334.480952 53.316667
	+L 288 53.316667
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ <path clip-path="url(#pd28355a76e)" d="M 334.480952 99.797619
	+L 380.961905 99.797619
	+L 380.961905 53.316667
	+L 334.480952 53.316667
	+z
	+" style="fill:#fff5f0;stroke:#000000;"/>
	+ </g>
	+ <g id="patch_3">
	+ <path d="M 176.445714 411.22
	+L 176.445714 20.78
	+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
	+ </g>
	+ <g id="patch_4">
	+ <path d="M 399.554286 411.22
	+L 399.554286 20.78
	+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
	+ </g>
	+ <g id="patch_5">
	+ <path d="M 176.445714 411.22
	+L 399.554286 411.22
	+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
	+ </g>
	+ <g id="patch_6">
	+ <path d="M 176.445714 20.78
	+L 399.554286 20.78
	+" style="fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;stroke-width:0.8;"/>
	+ </g>
	+ </g>
	+ </g>
	+ <defs>
	+ <clipPath id="pd28355a76e">
	+ <rect height="390.44" width="223.108571" x="176.445714" y="20.78"/>
	+ </clipPath>
	+ </defs>
	+</svg>
	diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h
	index ffb5d84..d42ab2c 100644
	--- a/include/GooseFEM/MeshQuad4.h
	+++ b/include/GooseFEM/MeshQuad4.h
	@@ -1,269 +1,269 @@
	/*

	(c - GPLv3) T.W.J. de Geus (Tom) \| tom@geus.me \| www.geus.me \| github.com/tdegeus/GooseFEM

	*/

	#ifndef GOOSEFEM_MESHQUAD4_H
	#define GOOSEFEM_MESHQUAD4_H

	#include "config.h"

	namespace GooseFEM {
	namespace Mesh {
	namespace Quad4 {

	namespace Map {
	class FineLayer2Regular;
	} // namespace Map

	class Regular {
	public:
	Regular() = default;
	Regular(size_t nelx, size_t nely, double h = 1.0);

	// size
	size_t nelem() const; // number of elements
	size_t nnode() const; // number of nodes
	size_t nne() const; // number of nodes-per-element
	size_t ndim() const; // number of dimensions
	size_t nelx() const; // number of elements in x-direction
	size_t nely() const; // number of elements in y-direction
	double h() const; // edge size

	// type
	ElementType getElementType() const;

	// mesh
	xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
	xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]

	// boundary nodes: edges
	xt::xtensor<size_t, 1> nodesBottomEdge() const;
	xt::xtensor<size_t, 1> nodesTopEdge() const;
	xt::xtensor<size_t, 1> nodesLeftEdge() const;
	xt::xtensor<size_t, 1> nodesRightEdge() const;

	// boundary nodes: edges, without corners
	xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
	xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
	xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
	xt::xtensor<size_t, 1> nodesRightOpenEdge() const;

	// boundary nodes: corners (including aliases)
	size_t nodesBottomLeftCorner() const;
	size_t nodesBottomRightCorner() const;
	size_t nodesTopLeftCorner() const;
	size_t nodesTopRightCorner() const;
	size_t nodesLeftBottomCorner() const;
	size_t nodesLeftTopCorner() const;
	size_t nodesRightBottomCorner() const;
	size_t nodesRightTopCorner() const;

	// DOF-numbers for each component of each node (sequential)
	xt::xtensor<size_t, 2> dofs() const;

	// DOF-numbers for the case that the periodicity if fully eliminated
	xt::xtensor<size_t, 2> dofsPeriodic() const;

	// periodic node pairs [:,2]: (independent, dependent)
	xt::xtensor<size_t, 2> nodesPeriodic() const;

	// front-bottom-left node, used as reference for periodicity
	size_t nodesOrigin() const;

	// element numbers as matrix
	xt::xtensor<size_t, 2> elementgrid() const;

	private:
	double m_h; // elementary element edge-size (in all directions)
	size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h")
	size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h")
	size_t m_nelem; // number of elements
	size_t m_nnode; // number of nodes
	static const size_t m_nne = 4; // number of nodes-per-element
	static const size_t m_ndim = 2; // number of dimensions
	};

	class FineLayer {
	public:
	FineLayer() = default;
	FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1);

	// Reconstruct class for given coordinates / connectivity
	FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);

	// size
	size_t nelem() const; // number of elements
	size_t nnode() const; // number of nodes
	size_t nne() const; // number of nodes-per-element
	size_t ndim() const; // number of dimensions
	size_t nelx() const; // number of elements in x-direction
	size_t nely() const; // number of elements in y-direction
	double h() const; // edge size

	// type
	ElementType getElementType() const;

	// mesh
	xt::xtensor<double, 2> coor() const; // nodal positions [nnode, ndim]
	xt::xtensor<size_t, 2> conn() const; // connectivity [nelem, nne]

	- // element sets
	- // - elements in the middle (fine) layer
	+ // elements in the middle (fine) layer
	xt::xtensor<size_t, 1> elementsMiddleLayer() const;
	- // - select region of elements from 'matrix' of element numbers
	+
	+ // select region of elements from 'matrix' of element numbers
	xt::xtensor<size_t, 1> elementgrid_ravel(
	- std::array<size_t, 2> start_stop_rows,
	- std::array<size_t, 2> start_stop_cols) const;
	- // - select region of elements from 'matrix' of element numbers around an element
	+ std::vector<size_t> rows_range,
	+ std::vector<size_t> cols_range) const;
	+
	+ // select region of elements from 'matrix' of element numbers around an element
	xt::xtensor<size_t, 1> elementgrid_around_ravel(
	size_t element,
	- std::array<int, 2> start_stop_rows,
	- std::array<int, 2> start_stop_cols) const;
	- // xt::xtensor<size_t, 1> elementsAround(size_t e, int left, int right) const; // region-of-interests
	+ size_t size,
	+ bool periodic = true);

	// boundary nodes: edges
	xt::xtensor<size_t, 1> nodesBottomEdge() const;
	xt::xtensor<size_t, 1> nodesTopEdge() const;
	xt::xtensor<size_t, 1> nodesLeftEdge() const;
	xt::xtensor<size_t, 1> nodesRightEdge() const;

	// boundary nodes: edges, without corners
	xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
	xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
	xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
	xt::xtensor<size_t, 1> nodesRightOpenEdge() const;

	// boundary nodes: corners (including aliases)
	size_t nodesBottomLeftCorner() const;
	size_t nodesBottomRightCorner() const;
	size_t nodesTopLeftCorner() const;
	size_t nodesTopRightCorner() const;
	size_t nodesLeftBottomCorner() const;
	size_t nodesLeftTopCorner() const;
	size_t nodesRightBottomCorner() const;
	size_t nodesRightTopCorner() const;

	// DOF-numbers for each component of each node (sequential)
	xt::xtensor<size_t, 2> dofs() const;

	// DOF-numbers for the case that the periodicity if fully eliminated
	xt::xtensor<size_t, 2> dofsPeriodic() const;

	// periodic node pairs [:,2]: (independent, dependent)
	xt::xtensor<size_t, 2> nodesPeriodic() const;

	// front-bottom-left node, used as reference for periodicity
	size_t nodesOrigin() const;

	// mapping to 'roll' periodically in the x-direction,
	// returns element mapping, such that: new_elemvar = elemvar[elem_map]
	xt::xtensor<size_t, 1> roll(size_t n);

	private:
	double m_h; // elementary element edge-size (in all directions)
	double m_Lx; // mesh size in "x"
	size_t m_nelem; // number of elements
	size_t m_nnode; // number of nodes
	static const size_t m_nne = 4; // number of nodes-per-element
	static const size_t m_ndim = 2; // number of dimensions
	xt::xtensor<size_t, 1> m_nelx; // number of elements in "x" (*)
	xt::xtensor<size_t, 1> m_nnd; // total number of nodes in the main node layer (**)
	xt::xtensor<size_t, 1> m_nhx; // element size in x-direction (*)
	xt::xtensor<size_t, 1> m_nhy; // element size in y-direction (*)
	xt::xtensor<int, 1> m_refine; // refine direction (-1:no refine, 0:"x" (*)
	xt::xtensor<size_t, 1> m_startElem; // start element (*)
	xt::xtensor<size_t, 1> m_startNode; // start node (**)
	// (*) per element layer in "y"
	// (**) per node layer in "y"

	void init(size_t nelx, size_t nely, double h, size_t nfine = 1);
	void map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);

	friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular;
	};

	namespace Map {

	// Return "FineLayer"-class responsible for generating a connectivity
	// Throws if conversion is not possible
	GooseFEM::Mesh::Quad4::FineLayer FineLayer(
	const xt::xtensor<double, 2>& coor,
	const xt::xtensor<size_t, 2>& conn);

	class RefineRegular {
	public:
	// constructor
	RefineRegular() = default;
	RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny);

	// return the one of the two meshes
	GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const;
	GooseFEM::Mesh::Quad4::Regular getFineMesh() const;

	// elements of the Fine mesh per element of the Coarse mesh
	xt::xtensor<size_t, 2> getMap() const;

	// map field
	xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 1>& data) const; // scalar per el
	xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
	xt::xtensor<double, 4> mapToCoarse(const xt::xtensor<double, 4>& data) const; // tensor per intpnt

	// map field
	xt::xtensor<double, 1> mapToFine(const xt::xtensor<double, 1>& data) const; // scalar per el
	xt::xtensor<double, 2> mapToFine(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
	xt::xtensor<double, 4> mapToFine(const xt::xtensor<double, 4>& data) const; // tensor per intpnt

	private:
	// the meshes
	GooseFEM::Mesh::Quad4::Regular m_coarse;
	GooseFEM::Mesh::Quad4::Regular m_fine;

	// mapping
	xt::xtensor<size_t, 1> m_fine2coarse;
	xt::xtensor<size_t, 1> m_fine2coarse_index;
	xt::xtensor<size_t, 2> m_coarse2fine;
	};

	class FineLayer2Regular {
	public:
	// constructor
	FineLayer2Regular() = default;
	FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh);

	// return either of the meshes
	GooseFEM::Mesh::Quad4::Regular getRegularMesh() const;
	GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const;

	// elements of the Regular mesh per element of the FineLayer mesh
	// and the fraction by which the overlap is
	std::vector<std::vector<size_t>> getMap() const;
	std::vector<std::vector<double>> getMapFraction() const;

	// map field
	xt::xtensor<double, 1> mapToRegular(const xt::xtensor<double, 1>& data) const; // scalar per el
	xt::xtensor<double, 2> mapToRegular(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
	xt::xtensor<double, 4> mapToRegular(const xt::xtensor<double, 4>& data) const; // tensor per intpnt

	private:
	// the "FineLayer" mesh to map
	GooseFEM::Mesh::Quad4::FineLayer m_finelayer;

	// the new "Regular" mesh to which to map
	GooseFEM::Mesh::Quad4::Regular m_regular;

	// mapping
	std::vector<std::vector<size_t>> m_elem_regular;
	std::vector<std::vector<double>> m_frac_regular;
	};

	} // namespace Map

	} // namespace Quad4
	} // namespace Mesh
	} // namespace GooseFEM

	#include "MeshQuad4.hpp"

	#endif
	diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp
	index 620dd35..4f21b26 100644
	--- a/include/GooseFEM/MeshQuad4.hpp
	+++ b/include/GooseFEM/MeshQuad4.hpp
	@@ -1,1558 +1,1560 @@
	/*

	(c - GPLv3) T.W.J. de Geus (Tom) \| tom@geus.me \| www.geus.me \| github.com/tdegeus/GooseFEM

	*/

	#ifndef GOOSEFEM_MESHQUAD4_HPP
	#define GOOSEFEM_MESHQUAD4_HPP

	#include "MeshQuad4.h"

	namespace GooseFEM {
	namespace Mesh {
	namespace Quad4 {

	inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely)
	{
	GOOSEFEM_ASSERT(m_nelx >= 1ul);
	GOOSEFEM_ASSERT(m_nely >= 1ul);

	m_nnode = (m_nelx + 1) * (m_nely + 1);
	m_nelem = m_nelx * m_nely;
	}

	inline size_t Regular::nelem() const
	{
	return m_nelem;
	}

	inline size_t Regular::nnode() const
	{
	return m_nnode;
	}

	inline size_t Regular::nne() const
	{
	return m_nne;
	}

	inline size_t Regular::ndim() const
	{
	return m_ndim;
	}

	inline size_t Regular::nelx() const
	{
	return m_nelx;
	}

	inline size_t Regular::nely() const
	{
	return m_nely;
	}

	inline double Regular::h() const
	{
	return m_h;
	}

	inline ElementType Regular::getElementType() const
	{
	return ElementType::Quad4;
	}

	inline xt::xtensor<double, 2> Regular::coor() const
	{
	xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});

	xt::xtensor<double, 1> x =
	xt::linspace<double>(0.0, m_h * static_cast<double>(m_nelx), m_nelx + 1);
	xt::xtensor<double, 1> y =
	xt::linspace<double>(0.0, m_h * static_cast<double>(m_nely), m_nely + 1);

	size_t inode = 0;

	for (size_t iy = 0; iy < m_nely + 1; ++iy) {
	for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
	ret(inode, 0) = x(ix);
	ret(inode, 1) = y(iy);
	++inode;
	}
	}

	return ret;
	}

	inline xt::xtensor<size_t, 2> Regular::conn() const
	{
	xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});

	size_t ielem = 0;

	for (size_t iy = 0; iy < m_nely; ++iy) {
	for (size_t ix = 0; ix < m_nelx; ++ix) {
	ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix);
	ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1);
	ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix);
	ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1);
	++ielem;
	}
	}

	return ret;
	}

	inline xt::xtensor<size_t, 1> Regular::nodesBottomEdge() const
	{
	return xt::arange<size_t>(m_nelx + 1);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesTopEdge() const
	{
	return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesLeftEdge() const
	{
	return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesRightEdge() const
	{
	return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
	}

	inline xt::xtensor<size_t, 1> Regular::nodesBottomOpenEdge() const
	{
	return xt::arange<size_t>(1, m_nelx);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesTopOpenEdge() const
	{
	return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesLeftOpenEdge() const
	{
	return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
	}

	inline xt::xtensor<size_t, 1> Regular::nodesRightOpenEdge() const
	{
	return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
	}

	inline size_t Regular::nodesBottomLeftCorner() const
	{
	return 0;
	}

	inline size_t Regular::nodesBottomRightCorner() const
	{
	return m_nelx;
	}

	inline size_t Regular::nodesTopLeftCorner() const
	{
	return m_nely * (m_nelx + 1);
	}

	inline size_t Regular::nodesTopRightCorner() const
	{
	return m_nely * (m_nelx + 1) + m_nelx;
	}

	inline size_t Regular::nodesLeftBottomCorner() const
	{
	return nodesBottomLeftCorner();
	}

	inline size_t Regular::nodesLeftTopCorner() const
	{
	return nodesTopLeftCorner();
	}

	inline size_t Regular::nodesRightBottomCorner() const
	{
	return nodesBottomRightCorner();
	}

	inline size_t Regular::nodesRightTopCorner() const
	{
	return nodesTopRightCorner();
	}

	inline xt::xtensor<size_t, 2> Regular::nodesPeriodic() const
	{
	xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
	xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
	xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
	xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
	std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
	xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);

	ret(0, 0) = nodesBottomLeftCorner();
	ret(0, 1) = nodesBottomRightCorner();

	ret(1, 0) = nodesBottomLeftCorner();
	ret(1, 1) = nodesTopRightCorner();

	ret(2, 0) = nodesBottomLeftCorner();
	ret(2, 1) = nodesTopLeftCorner();

	size_t i = 3;

	xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
	xt::view(ret, xt::range(i, i + bot.size()), 1) = top;

	i += bot.size();

	xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
	xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;

	return ret;
	}

	inline size_t Regular::nodesOrigin() const
	{
	return nodesBottomLeftCorner();
	}

	inline xt::xtensor<size_t, 2> Regular::dofs() const
	{
	return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
	}

	inline xt::xtensor<size_t, 2> Regular::dofsPeriodic() const
	{
	xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
	xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
	xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
	xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);

	for (size_t j = 0; j < m_ndim; ++j) {
	xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
	}

	return GooseFEM::Mesh::renumber(ret);
	}

	inline xt::xtensor<size_t, 2> Regular::elementgrid() const
	{
	return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
	}

	inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine)
	{
	this->init(nelx, nely, h, nfine);
	}

	inline FineLayer::FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
	{
	this->map(coor, conn);
	}

	inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine)
	{
	GOOSEFEM_ASSERT(nelx >= 1ul);
	GOOSEFEM_ASSERT(nely >= 1ul);

	m_h = h;
	m_Lx = m_h * static_cast<double>(nelx);

	// compute element size in y-direction (use symmetry, compute upper half)

	// temporary variables
	size_t nmin, ntot;
	xt::xtensor<size_t, 1> nhx = xt::ones<size_t>({nely});
	xt::xtensor<size_t, 1> nhy = xt::ones<size_t>({nely});
	xt::xtensor<int, 1> refine = -1 * xt::ones<int>({nely});

	// minimum height in y-direction (half of the height because of symmetry)
	if (nely % 2 == 0) {
	nmin = nely / 2;
	}
	else {
	nmin = (nely + 1) / 2;
	}

	// minimum number of fine layers in y-direction (minimum 1, middle layer part of this half)
	if (nfine % 2 == 0) {
	nfine = nfine / 2 + 1;
	}
	else {
	nfine = (nfine + 1) / 2;
	}

	if (nfine < 1) {
	nfine = 1;
	}

	if (nfine > nmin) {
	nfine = nmin;
	}

	// loop over element layers in y-direction, try to coarsen using these rules:
	// (1) element size in y-direction <= distance to origin in y-direction
	// (2) element size in x-direction should fit the total number of elements in x-direction
	// (3) a certain number of layers have the minimum size "1" (are fine)
	for (size_t iy = nfine;;) {
	// initialize current size in y-direction
	if (iy == nfine) {
	ntot = nfine;
	}
	// check to stop
	if (iy >= nely \|\| ntot >= nmin) {
	nely = iy;
	break;
	}

	// rules (1,2) satisfied: coarsen in x-direction
	if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
	refine(iy) = 0;
	nhy(iy) *= 2;
	auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
	auto vnhx = xt::view(nhx, xt::range(iy, _));
	vnhy *= 3;
	vnhx *= 3;
	}

	// update the number of elements in y-direction
	ntot += nhy(iy);
	// proceed to next element layer in y-direction
	++iy;
	// check to stop
	if (iy >= nely \|\| ntot >= nmin) {
	nely = iy;
	break;
	}
	}

	// symmetrize, compute full information

	// allocate mesh constructor parameters
	m_nhx = xt::empty<size_t>({nely * 2 - 1});
	m_nhy = xt::empty<size_t>({nely * 2 - 1});
	m_refine = xt::empty<int>({nely * 2 - 1});
	m_nelx = xt::empty<size_t>({nely * 2 - 1});
	m_nnd = xt::empty<size_t>({nely * 2});
	m_startElem = xt::empty<size_t>({nely * 2 - 1});
	m_startNode = xt::empty<size_t>({nely * 2});

	// fill
	// - lower half
	for (size_t iy = 0; iy < nely; ++iy) {
	m_nhx(iy) = nhx(nely - iy - 1);
	m_nhy(iy) = nhy(nely - iy - 1);
	m_refine(iy) = refine(nely - iy - 1);
	}
	// - upper half
	for (size_t iy = 0; iy < nely - 1; ++iy) {
	m_nhx(iy + nely) = nhx(iy + 1);
	m_nhy(iy + nely) = nhy(iy + 1);
	m_refine(iy + nely) = refine(iy + 1);
	}

	// update size
	nely = m_nhx.size();

	// compute the number of elements per element layer in y-direction
	for (size_t iy = 0; iy < nely; ++iy) {
	m_nelx(iy) = nelx / m_nhx(iy);
	}

	// compute the number of nodes per node layer in y-direction
	for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
	m_nnd(iy) = m_nelx(iy) + 1;
	}
	for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
	m_nnd(iy + 1) = m_nelx(iy) + 1;
	}

	// compute mesh dimensions

	// initialize
	m_nnode = 0;
	m_nelem = 0;
	m_startNode(0) = 0;

	// loop over element layers (bottom -> middle, elements become finer)
	for (size_t i = 0; i < (nely - 1) / 2; ++i) {
	// - store the first element of the layer
	m_startElem(i) = m_nelem;
	// - add the nodes of this layer
	if (m_refine(i) == 0) {
	m_nnode += (3 * m_nelx(i) + 1);
	}
	else {
	m_nnode += (m_nelx(i) + 1);
	}
	// - add the elements of this layer
	if (m_refine(i) == 0) {
	m_nelem += (4 * m_nelx(i));
	}
	else {
	m_nelem += (m_nelx(i));
	}
	// - store the starting node of the next layer
	m_startNode(i + 1) = m_nnode;
	}

	// loop over element layers (middle -> top, elements become coarser)
	for (size_t i = (nely - 1) / 2; i < nely; ++i) {
	// - store the first element of the layer
	m_startElem(i) = m_nelem;
	// - add the nodes of this layer
	if (m_refine(i) == 0) {
	m_nnode += (5 * m_nelx(i) + 1);
	}
	else {
	m_nnode += (m_nelx(i) + 1);
	}
	// - add the elements of this layer
	if (m_refine(i) == 0) {
	m_nelem += (4 * m_nelx(i));
	}
	else {
	m_nelem += (m_nelx(i));
	}
	// - store the starting node of the next layer
	m_startNode(i + 1) = m_nnode;
	}
	// - add the top row of nodes
	m_nnode += m_nelx(nely - 1) + 1;
	}

	inline size_t FineLayer::nelem() const
	{
	return m_nelem;
	}

	inline size_t FineLayer::nnode() const
	{
	return m_nnode;
	}

	inline size_t FineLayer::nne() const
	{
	return m_nne;
	}

	inline size_t FineLayer::ndim() const
	{
	return m_ndim;
	}

	inline size_t FineLayer::nelx() const
	{
	return xt::amax(m_nelx)();
	}

	inline size_t FineLayer::nely() const
	{
	return xt::sum(m_nhy)();
	}

	inline double FineLayer::h() const
	{
	return m_h;
	}

	inline ElementType FineLayer::getElementType() const
	{
	return ElementType::Quad4;
	}

	inline xt::xtensor<double, 2> FineLayer::coor() const
	{
	// allocate output
	xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});

	// current node, number of element layers
	size_t inode = 0;
	size_t nely = static_cast<size_t>(m_nhy.size());

	// y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
	// - allocate
	xt::xtensor<double, 1> y = xt::empty<double>({nely + 1});
	// - initialize
	y(0) = 0.0;
	// - compute
	for (size_t iy = 1; iy < nely + 1; ++iy) {
	y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h;
	}

	// loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes

	for (size_t iy = 0;; ++iy) {
	// get positions along the x- and z-axis
	xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);

	// add nodes of the bottom layer of this element
	for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
	ret(inode, 0) = x(ix);
	ret(inode, 1) = y(iy);
	++inode;
	}

	// stop at middle layer
	if (iy == (nely - 1) / 2) {
	break;
	}

	// add extra nodes of the intermediate layer, for refinement in x-direction
	if (m_refine(iy) == 0) {
	// - get position offset in x- and y-direction
	double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
	double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
	// - add nodes of the intermediate layer
	for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
	for (size_t j = 0; j < 2; ++j) {
	ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
	ret(inode, 1) = y(iy) + dy;
	++inode;
	}
	}
	}
	}

	// loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes

	for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
	// get positions along the x- and z-axis
	xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);

	// add extra nodes of the intermediate layer, for refinement in x-direction
	if (m_refine(iy) == 0) {
	// - get position offset in x- and y-direction
	double dx = m_h * static_cast<double>(m_nhx(iy) / 3);
	double dy = m_h * static_cast<double>(m_nhy(iy) / 2);
	// - add nodes of the intermediate layer
	for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
	for (size_t j = 0; j < 2; ++j) {
	ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
	ret(inode, 1) = y(iy) + dy;
	++inode;
	}
	}
	}

	// add nodes of the top layer of this element
	for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
	ret(inode, 0) = x(ix);
	ret(inode, 1) = y(iy + 1);
	++inode;
	}
	}

	return ret;
	}

	inline xt::xtensor<size_t, 2> FineLayer::conn() const
	{
	// allocate output
	xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});

	// current element, number of element layers, starting nodes of each node layer
	size_t ielem = 0;
	size_t nely = static_cast<size_t>(m_nhy.size());
	size_t bot, mid, top;

	// loop over all element layers
	for (size_t iy = 0; iy < nely; ++iy) {
	// - get: starting nodes of bottom(, middle) and top layer
	bot = m_startNode(iy);
	mid = m_startNode(iy) + m_nnd(iy);
	top = m_startNode(iy + 1);

	// - define connectivity: no coarsening/refinement
	if (m_refine(iy) == -1) {
	for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
	ret(ielem, 0) = bot + (ix);
	ret(ielem, 1) = bot + (ix + 1);
	ret(ielem, 2) = top + (ix + 1);
	ret(ielem, 3) = top + (ix);
	ielem++;
	}
	}

	// - define connectivity: refinement along the x-direction (below the middle layer)
	else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
	for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
	// -- bottom element
	ret(ielem, 0) = bot + (ix);
	ret(ielem, 1) = bot + (ix + 1);
	ret(ielem, 2) = mid + (2 * ix + 1);
	ret(ielem, 3) = mid + (2 * ix);
	ielem++;
	// -- top-right element
	ret(ielem, 0) = bot + (ix + 1);
	ret(ielem, 1) = top + (3 * ix + 3);
	ret(ielem, 2) = top + (3 * ix + 2);
	ret(ielem, 3) = mid + (2 * ix + 1);
	ielem++;
	// -- top-center element
	ret(ielem, 0) = mid + (2 * ix);
	ret(ielem, 1) = mid + (2 * ix + 1);
	ret(ielem, 2) = top + (3 * ix + 2);
	ret(ielem, 3) = top + (3 * ix + 1);
	ielem++;
	// -- top-left element
	ret(ielem, 0) = bot + (ix);
	ret(ielem, 1) = mid + (2 * ix);
	ret(ielem, 2) = top + (3 * ix + 1);
	ret(ielem, 3) = top + (3 * ix);
	ielem++;
	}
	}

	// - define connectivity: coarsening along the x-direction (above the middle layer)
	else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) {
	for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
	// -- lower-left element
	ret(ielem, 0) = bot + (3 * ix);
	ret(ielem, 1) = bot + (3 * ix + 1);
	ret(ielem, 2) = mid + (2 * ix);
	ret(ielem, 3) = top + (ix);
	ielem++;
	// -- lower-center element
	ret(ielem, 0) = bot + (3 * ix + 1);
	ret(ielem, 1) = bot + (3 * ix + 2);
	ret(ielem, 2) = mid + (2 * ix + 1);
	ret(ielem, 3) = mid + (2 * ix);
	ielem++;
	// -- lower-right element
	ret(ielem, 0) = bot + (3 * ix + 2);
	ret(ielem, 1) = bot + (3 * ix + 3);
	ret(ielem, 2) = top + (ix + 1);
	ret(ielem, 3) = mid + (2 * ix + 1);
	ielem++;
	// -- upper element
	ret(ielem, 0) = mid + (2 * ix);
	ret(ielem, 1) = mid + (2 * ix + 1);
	ret(ielem, 2) = top + (ix + 1);
	ret(ielem, 3) = top + (ix);
	ielem++;
	}
	}
	}

	return ret;
	}

	inline xt::xtensor<size_t, 1> FineLayer::elementsMiddleLayer() const
	{
	size_t nely = m_nhy.size();
	size_t iy = (nely - 1) / 2;
	return m_startElem(iy) + xt::arange<size_t>(m_nelx(iy));
	}

	inline xt::xtensor<size_t, 1> FineLayer::elementgrid_ravel(
	- std::array<size_t, 2> start_stop_rows,
	- std::array<size_t, 2> start_stop_cols) const
	+ std::vector<size_t> start_stop_rows,
	+ std::vector<size_t> start_stop_cols) const
	{
	- GOOSEFEM_ASSERT(start_stop_rows[0] <= this->nely());
	- GOOSEFEM_ASSERT(start_stop_rows[1] <= this->nely());
	- GOOSEFEM_ASSERT(start_stop_cols[0] <= this->nelx());
	- GOOSEFEM_ASSERT(start_stop_cols[1] <= this->nelx());
	+ GOOSEFEM_ASSERT(start_stop_rows.size() == 0 \|\| start_stop_rows.size() == 2);
	+ GOOSEFEM_ASSERT(start_stop_cols.size() == 0 \|\| start_stop_cols.size() == 2);

	- if (start_stop_rows[0] == start_stop_rows[1] \|\| start_stop_cols[0] == start_stop_cols[1]) {
	+ std::array<size_t, 2> rows;
	+ std::array<size_t, 2> cols;
	+
	+ if (start_stop_rows.size() == 2) {
	+ std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin());
	+ GOOSEFEM_ASSERT(rows[0] <= this->nely());
	+ GOOSEFEM_ASSERT(rows[1] <= this->nely());
	+ }
	+ else {
	+ rows[0] = 0;
	+ rows[1] = this->nely();
	+ }
	+
	+ if (start_stop_cols.size() == 2) {
	+ std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin());
	+ GOOSEFEM_ASSERT(cols[0] <= this->nelx());
	+ GOOSEFEM_ASSERT(cols[1] <= this->nelx());
	+ }
	+ else {
	+ cols[0] = 0;
	+ cols[1] = this->nelx();
	+ }
	+
	+ if (rows[0] == rows[1] \|\| cols[0] == cols[1]) {
	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({0});
	return ret;
	}

	// Compute dimensions

	auto H = xt::cumsum(m_nhy);
	size_t yl = 0;
	- if (start_stop_rows[0] > 0) {
	- yl = xt::argmax(H > start_stop_rows[0])();
	+ if (rows[0] > 0) {
	+ yl = xt::argmax(H > rows[0])();
	}
	- size_t yu = xt::argmax(H >= start_stop_rows[1])();
	+ size_t yu = xt::argmax(H >= rows[1])();
	size_t hx = std::max(m_nhx(yl), m_nhx(yu));
	- size_t xl = (start_stop_cols[0] - start_stop_cols[0] % hx) / hx;
	- size_t xu = (start_stop_cols[1] - start_stop_cols[1] % hx) / hx;
	+ size_t xl = (cols[0] - cols[0] % hx) / hx;
	+ size_t xu = (cols[1] - cols[1] % hx) / hx;

	// Allocate output

	size_t N = 0;

	for (size_t i = yl; i <= yu; ++i) {
	// no refinement
	size_t n = (xu - xl) * hx / m_nhx(i);
	// refinement
	if (m_refine(i) != -1) {
	n *= 4;
	}
	N += n;
	}

	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({N});

	// Write output

	N = 0;

	for (size_t i = yl; i <= yu; ++i) {
	// no refinement
	size_t n = (xu - xl) * hx / m_nhx(i);
	+ size_t h = hx;
	// refinement
	if (m_refine(i) != -1) {
	n *= 4;
	+ h *= 4;
	}
	- xt::xtensor<size_t, 1> e = m_startElem(i) + xl * hx / m_nhx(i) + xt::arange<size_t>(n);
	+ xt::xtensor<size_t, 1> e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange<size_t>(n);
	xt::view(ret, xt::range(N, N + n)) = e;
	N += n;
	}

	return ret;
	}

	+inline xt::xtensor<size_t, 1> FineLayer::elementgrid_around_ravel(
	+ size_t e,
	+ size_t size,
	+ bool periodic)
	+{
	+ GOOSEFEM_WIP_ASSERT(periodic == true);
	+
	+ size_t iy = xt::argmin(m_startElem <= e)() - 1;
	+ size_t nel = m_nelx(iy);
	+
	+ GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
	+
	+ auto H = xt::cumsum(m_nhy);

	-// inline xt::xtensor<size_t, 1> FineLayer::elementsAround(size_t e, int left, int right) const
	-// {
	-// GOOSEFEM_ASSERT(right - left <= static_cast<decltype(left)>(this->nelx()));
	-// GOOSEFEM_ASSERT(right - left <= static_cast<decltype(left)>(this->nely()));
	-// GOOSEFEM_ASSERT(e < this->nelem());
	-// GOOSEFEM_WIP_ASSERT(left <= 0);
	-// GOOSEFEM_WIP_ASSERT(right >= 0);
	-
	-// size_t l = static_cast<size_t>(std::abs(left));
	-// size_t r = static_cast<size_t>(std::abs(right));
	-// size_t w = l + r;
	-
	-// size_t nely = m_nhy.size();
	-// size_t iy = xt::argmin(xt::greater(e, m_startElem))() - 1;
	-// size_t nel = m_nelx(iy);
	-// size_t mid = static_cast<size_t>(static_cast<int>(m_startElem(iy)) - left);
	-// size_t nroll = nel - (nel - nel % 2) / 2 + e - m_startElem(iy);
	-// GOOSEFEM_WIP_ASSERT(iy == (nely - 1) / 2);
	-// GOOSEFEM_WIP_ASSERT(iy >= l);
	-// GOOSEFEM_WIP_ASSERT(iy < r);
	-
	-// // Select elements layers around the middle layer
	-
	-// size_t nyb = 0;
	-// size_t nyt = 0;
	-// size_t nhb = 0;
	-// size_t nht = m_nhy(iy);
	-// if (left < 0) {
	-// while (nhb < l) {
	-// nyb += 1;
	-// nhb += m_nhy(iy - nyb);
	-// }
	-// }
	-// if (right > 0) {
	-// while (nht < r) {
	-// nyt += 1;
	-// nht += m_nhy(iy + nyt);
	-// }
	-// }
	-// nyt++;
	-
	-// // Allocate output
	-
	-// size_t N = 0;
	-
	-// for (size_t i = iy - nyb; i < iy + nyt; ++i) {
	-// if (m_refine(iy) == -1) {
	-// N += w / m_nhx(i); // no refinement
	-// }
	-// else {
	-// N += 4 * w / m_nhx(i); // refinement
	-// }
	-// }
	-
	-// xt::xtensor<size_t, 1> ret = xt::empty<size_t>({N});
	-
	-// // Fill output
	-
	-// N = 0;
	-// size_t n;
	-
	-// std::cout << w << std::endl;
	-// std::cout << m_nhx << std::endl;
	-
	-// for (size_t i = iy - nyb; i < iy + nyt; ++i) {
	-// if (m_refine(iy) == -1) {
	-// n = w / m_nhx(i); // no refinement
	-// }
	-// else {
	-// n = 4 * w / m_nhx(i); // refinement
	-// }
	-// xt::view(ret, xt::range(N, N + n)) = m_startElem(i) + xt::arange<size_t>(n);
	-// N += n;
	-// }
	-
	-// return ret;
	-// }
	+ if (2 * size >= H(H.size() - 1)) {
	+ return xt::arange<size_t>(this->nelem());
	+ }
	+
	+ size_t hy = H(iy);
	+ size_t l = xt::argmax(H > (hy - size - 1))();
	+ size_t u = xt::argmax(H >= (hy + size))();
	+ size_t lh = 0;
	+ if (l > 0) {
	+ lh = H(l - 1);
	+ }
	+ size_t uh = H(u);
	+
	+ size_t step = xt::amax(m_nhx)();
	+ size_t relx = (e - m_startElem(iy)) % step;
	+ size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
	+ size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
	+ size_t dx = m_nhx(u);
	+ size_t xl = mid - size;
	+ size_t xu = mid + size + 1;
	+ xl = xl - xl % dx;
	+ xu = xu - xu % dx;
	+ if (mid - xl < size) {
	+ if (xl < dx) {
	+ xl = 0;
	+ }
	+ else {
	+ xl -= dx;
	+ }
	+ }
	+ if (xu - mid < size) {
	+ if (xu > nel - dx) {
	+ xu = nel;
	+ }
	+ else {
	+ xu += dx;
	+ }
	+ }
	+
	+ auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu});
	+ auto map = this->roll(nroll);
	+ return xt::view(map, xt::keep(ret));
	+}

	inline xt::xtensor<size_t, 1> FineLayer::nodesBottomEdge() const
	{
	return m_startNode(0) + xt::arange<size_t>(m_nelx(0) + 1);
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesTopEdge() const
	{
	size_t nely = m_nhy.size();
	return m_startNode(nely) + xt::arange<size_t>(m_nelx(nely - 1) + 1);
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesLeftEdge() const
	{
	size_t nely = m_nhy.size();

	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});

	size_t i = 0;
	size_t j = (nely + 1) / 2;
	size_t k = (nely - 1) / 2;
	size_t l = nely;

	xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j));
	xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1));

	return ret;
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesRightEdge() const
	{
	size_t nely = m_nhy.size();

	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});

	size_t i = 0;
	size_t j = (nely + 1) / 2;
	size_t k = (nely - 1) / 2;
	size_t l = nely;

	xt::view(ret, xt::range(i, j)) =
	xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j));

	xt::view(ret, xt::range(k + 1, l + 1)) =
	xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l));

	return ret;
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesBottomOpenEdge() const
	{
	return m_startNode(0) + xt::arange<size_t>(1, m_nelx(0));
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesTopOpenEdge() const
	{
	size_t nely = m_nhy.size();

	return m_startNode(nely) + xt::arange<size_t>(1, m_nelx(nely - 1));
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesLeftOpenEdge() const
	{
	size_t nely = m_nhy.size();

	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});

	size_t i = 0;
	size_t j = (nely + 1) / 2;
	size_t k = (nely - 1) / 2;
	size_t l = nely;

	xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j));
	xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l));

	return ret;
	}

	inline xt::xtensor<size_t, 1> FineLayer::nodesRightOpenEdge() const
	{
	size_t nely = m_nhy.size();

	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});

	size_t i = 0;
	size_t j = (nely + 1) / 2;
	size_t k = (nely - 1) / 2;
	size_t l = nely;

	xt::view(ret, xt::range(i, j - 1)) =
	xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j));

	xt::view(ret, xt::range(k, l - 1)) =
	xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1));

	return ret;
	}

	inline size_t FineLayer::nodesBottomLeftCorner() const
	{
	return m_startNode(0);
	}

	inline size_t FineLayer::nodesBottomRightCorner() const
	{
	return m_startNode(0) + m_nelx(0);
	}

	inline size_t FineLayer::nodesTopLeftCorner() const
	{
	size_t nely = m_nhy.size();

	return m_startNode(nely);
	}

	inline size_t FineLayer::nodesTopRightCorner() const
	{
	size_t nely = m_nhy.size();

	return m_startNode(nely) + m_nelx(nely - 1);
	}

	inline size_t FineLayer::nodesLeftBottomCorner() const
	{
	return nodesBottomLeftCorner();
	}

	inline size_t FineLayer::nodesRightBottomCorner() const
	{
	return nodesBottomRightCorner();
	}

	inline size_t FineLayer::nodesLeftTopCorner() const
	{
	return nodesTopLeftCorner();
	}

	inline size_t FineLayer::nodesRightTopCorner() const
	{
	return nodesTopRightCorner();
	}

	inline xt::xtensor<size_t, 2> FineLayer::nodesPeriodic() const
	{
	xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
	xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
	xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
	xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
	std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
	xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);

	ret(0, 0) = nodesBottomLeftCorner();
	ret(0, 1) = nodesBottomRightCorner();

	ret(1, 0) = nodesBottomLeftCorner();
	ret(1, 1) = nodesTopRightCorner();

	ret(2, 0) = nodesBottomLeftCorner();
	ret(2, 1) = nodesTopLeftCorner();

	size_t i = 3;

	xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
	xt::view(ret, xt::range(i, i + bot.size()), 1) = top;

	i += bot.size();

	xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
	xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;

	return ret;
	}

	inline size_t FineLayer::nodesOrigin() const
	{
	return nodesBottomLeftCorner();
	}

	inline xt::xtensor<size_t, 2> FineLayer::dofs() const
	{
	return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
	}

	inline xt::xtensor<size_t, 2> FineLayer::dofsPeriodic() const
	{
	xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
	xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
	xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
	xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);

	for (size_t j = 0; j < m_ndim; ++j) {
	xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
	}

	return GooseFEM::Mesh::renumber(ret);
	}

	inline xt::xtensor<size_t, 1> FineLayer::roll(size_t n)
	{
	auto conn = this->conn();
	size_t nely = static_cast<size_t>(m_nhy.size());
	xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelem});

	// loop over all element layers
	for (size_t iy = 0; iy < nely; ++iy) {

	// no refinement
	size_t shift = n * (m_nelx(iy) / m_nelx(0));
	size_t nel = m_nelx(iy);

	// refinement
	if (m_refine(iy) != -1) {
	shift = n * (m_nelx(iy) / m_nelx(0)) * 4;
	nel = m_nelx(iy) * 4;
	}

	// element numbers of the layer, and roll them
	auto e = m_startElem(iy) + xt::arange<size_t>(nel);
	xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift);
	}

	return ret;
	}

	inline void FineLayer::map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
	{
	GOOSEFEM_ASSERT(coor.shape(1) == 2);
	GOOSEFEM_ASSERT(conn.shape(1) == 4);
	GOOSEFEM_ASSERT(conn.shape(0) > 0);
	GOOSEFEM_ASSERT(coor.shape(0) >= 4);
	GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));

	if (conn.shape(0) == 1) {
	this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0));
	GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
	GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
	return;
	}

	// Identify the middle layer

	size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2;
	size_t eleft = emid;
	size_t eright = emid;

	while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) {
	eleft--;
	}

	while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) {
	eright++;
	}

	GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0));

	// Get element sizes along the middle layer

	auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0);
	auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1);
	auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2);
	auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0);
	auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1);
	auto hx = xt::amin(dx)();
	auto hy = xt::amin(dy)();

	GOOSEFEM_CHECK(xt::allclose(hx, hy));
	GOOSEFEM_CHECK(xt::allclose(dx, hx));
	GOOSEFEM_CHECK(xt::allclose(dy, hy));

	// Extract shape and initialise

	size_t nelx = eright - eleft + 1;
	size_t nely = static_cast<size_t>((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx);
	this->init(nelx, nely, hx);
	GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
	GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
	GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange<size_t>(nelx))));
	}

	namespace Map {

	inline RefineRegular::RefineRegular(
	const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny)
	: m_coarse(mesh)
	{
	m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h());

	xt::xtensor<size_t, 2> elmat_coarse = m_coarse.elementgrid();
	xt::xtensor<size_t, 2> elmat_fine = m_fine.elementgrid();

	m_coarse2fine = xt::empty<size_t>({m_coarse.nelem(), nx * ny});

	for (size_t i = 0; i < elmat_coarse.shape(0); ++i) {
	for (size_t j = 0; j < elmat_coarse.shape(1); ++j) {
	xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view(
	elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx)));
	}
	}
	}

	inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const
	{
	return m_coarse;
	}

	inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const
	{
	return m_fine;
	}

	inline xt::xtensor<size_t, 2> RefineRegular::getMap() const
	{
	return m_coarse2fine;
	}

	inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 1>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());

	size_t m = m_coarse2fine.shape(0);
	size_t n = m_coarse2fine.shape(1);

	xt::xtensor<double, 2> ret = xt::empty<double>({m, n});

	for (size_t i = 0; i < m; ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	xt::view(ret, i) = xt::view(data, xt::keep(e));
	}

	return ret;
	}

	inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 2>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());

	size_t m = m_coarse2fine.shape(0);
	size_t n = m_coarse2fine.shape(1);
	size_t N = data.shape(1);

	xt::xtensor<double, 2> ret = xt::empty<double>({m, n * N});

	for (size_t i = 0; i < m; ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	for (size_t q = 0; q < data.shape(1); ++q) {
	xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
	}
	}

	return ret;
	}

	inline xt::xtensor<double, 4> RefineRegular::mapToCoarse(const xt::xtensor<double, 4>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());

	size_t m = m_coarse2fine.shape(0);
	size_t n = m_coarse2fine.shape(1);
	size_t N = data.shape(1);

	xt::xtensor<double, 4> ret = xt::empty<double>({m, n * N, data.shape(2), data.shape(3)});

	for (size_t i = 0; i < m; ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	for (size_t q = 0; q < data.shape(1); ++q) {
	xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
	}
	}

	return ret;
	}

	inline xt::xtensor<double, 1> RefineRegular::mapToFine(const xt::xtensor<double, 1>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));

	xt::xtensor<double, 1> ret = xt::empty<double>({m_coarse2fine.size()});

	for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	xt::view(ret, xt::keep(e)) = data(i);
	}

	return ret;
	}

	inline xt::xtensor<double, 2> RefineRegular::mapToFine(const xt::xtensor<double, 2>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));

	xt::xtensor<double, 2> ret = xt::empty<double>({m_coarse2fine.size(), data.shape(1)});

	for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	xt::view(ret, xt::keep(e)) = xt::view(data, i);
	}

	return ret;
	}

	inline xt::xtensor<double, 4> RefineRegular::mapToFine(const xt::xtensor<double, 4>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));

	xt::xtensor<double, 4> ret =
	xt::empty<double>({m_coarse2fine.size(), data.shape(1), data.shape(2), data.shape(3)});

	for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
	auto e = xt::view(m_coarse2fine, i, xt::all());
	xt::view(ret, xt::keep(e)) = xt::view(data, i);
	}

	return ret;
	}

	inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh)
	: m_finelayer(mesh)
	{
	// ------------
	// Regular-mesh
	// ------------

	m_regular = GooseFEM::Mesh::Quad4::Regular(
	xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);

	// -------
	// mapping
	// -------

	// allocate mapping
	m_elem_regular.resize(m_finelayer.m_nelem);
	m_frac_regular.resize(m_finelayer.m_nelem);

	// alias
	xt::xtensor<size_t, 1> nhx = m_finelayer.m_nhx;
	xt::xtensor<size_t, 1> nhy = m_finelayer.m_nhy;
	xt::xtensor<size_t, 1> nelx = m_finelayer.m_nelx;
	xt::xtensor<size_t, 1> start = m_finelayer.m_startElem;

	// 'matrix' of element numbers of the Regular-mesh
	xt::xtensor<size_t, 2> elementgrid = m_regular.elementgrid();

	// cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh
	xt::xtensor<size_t, 1> cum_nhy =
	xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));

	// number of element layers in y-direction of the FineLayer-mesh
	size_t nely = nhy.size();

	// loop over layers of the FineLayer-mesh
	for (size_t iy = 0; iy < nely; ++iy) {
	// element numbers of the Regular-mesh along this layer of the FineLayer-mesh
	auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());

	// no coarsening/refinement
	// ------------------------

	if (m_finelayer.m_refine(iy) == -1) {
	// element numbers of the FineLayer-mesh along this layer
	xt::xtensor<size_t, 1> el_old = start(iy) + xt::arange<size_t>(nelx(iy));

	// loop along this layer of the FineLayer-mesh
	for (size_t ix = 0; ix < nelx(iy); ++ix) {
	// get the element numbers of the Regular-mesh for this element of the
	// FineLayer-mesh
	auto block =
	xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));

	// write to mapping
	for (auto& i : block) {
	m_elem_regular[el_old(ix)].push_back(i);
	m_frac_regular[el_old(ix)].push_back(1.0);
	}
	}
	}

	// refinement along the x-direction (below the middle layer)
	// ---------------------------------------------------------

	else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
	// element numbers of the FineLayer-mesh along this layer
	// rows: coarse block, columns element numbers per block
	xt::xtensor<size_t, 2> el_old =
	start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});

	// loop along this layer of the FineLayer-mesh
	for (size_t ix = 0; ix < nelx(iy); ++ix) {
	// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
	auto block =
	xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));

	// bottom: wide-to-narrow
	{
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));

	m_elem_regular[el_old(ix, 0)].push_back(e(0));
	m_frac_regular[el_old(ix, 0)].push_back(0.5);

	for (size_t k = 1; k < e.size() - 1; ++k) {
	m_elem_regular[el_old(ix, 0)].push_back(e(k));
	m_frac_regular[el_old(ix, 0)].push_back(1.0);
	}

	m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
	m_frac_regular[el_old(ix, 0)].push_back(0.5);
	}
	}

	// top: regular small element
	{
	auto e = xt::view(
	block,
	xt::range(nhy(iy) / 2, nhy(iy)),
	xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 2)].push_back(i);
	m_frac_regular[el_old(ix, 2)].push_back(1.0);
	}
	}

	// left
	{
	// left-bottom: narrow-to-wide
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e = xt::view(block, j, xt::range(0, j + 1));

	for (size_t k = 0; k < e.size() - 1; ++k) {
	m_elem_regular[el_old(ix, 3)].push_back(e(k));
	m_frac_regular[el_old(ix, 3)].push_back(1.0);
	}

	m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
	m_frac_regular[el_old(ix, 3)].push_back(0.5);
	}

	// left-top: regular
	{
	auto e = xt::view(
	block,
	xt::range(nhy(iy) / 2, nhy(iy)),
	xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 3)].push_back(i);
	m_frac_regular[el_old(ix, 3)].push_back(1.0);
	}
	}
	}

	// right
	{
	// right-bottom: narrow-to-wide
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy)));

	m_elem_regular[el_old(ix, 1)].push_back(e(0));
	m_frac_regular[el_old(ix, 1)].push_back(0.5);

	for (size_t k = 1; k < e.size(); ++k) {
	m_elem_regular[el_old(ix, 1)].push_back(e(k));
	m_frac_regular[el_old(ix, 1)].push_back(1.0);
	}
	}

	// right-top: regular
	{
	auto e = xt::view(
	block,
	xt::range(nhy(iy) / 2, nhy(iy)),
	xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 1)].push_back(i);
	m_frac_regular[el_old(ix, 1)].push_back(1.0);
	}
	}
	}
	}
	}

	// coarsening along the x-direction (above the middle layer)
	else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
	// element numbers of the FineLayer-mesh along this layer
	// rows: coarse block, columns element numbers per block
	xt::xtensor<size_t, 2> el_old =
	start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});

	// loop along this layer of the FineLayer-mesh
	for (size_t ix = 0; ix < nelx(iy); ++ix) {
	// get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
	auto block =
	xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));

	// top: narrow-to-wide
	{
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e = xt::view(
	block,
	nhy(iy) / 2 + j,
	xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));

	m_elem_regular[el_old(ix, 3)].push_back(e(0));
	m_frac_regular[el_old(ix, 3)].push_back(0.5);

	for (size_t k = 1; k < e.size() - 1; ++k) {
	m_elem_regular[el_old(ix, 3)].push_back(e(k));
	m_frac_regular[el_old(ix, 3)].push_back(1.0);
	}

	m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
	m_frac_regular[el_old(ix, 3)].push_back(0.5);
	}
	}

	// bottom: regular small element
	{
	auto e = xt::view(
	block,
	xt::range(0, nhy(iy) / 2),
	xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 1)].push_back(i);
	m_frac_regular[el_old(ix, 1)].push_back(1.0);
	}
	}

	// left
	{
	// left-bottom: regular
	{
	auto e = xt::view(
	block,
	xt::range(0, nhy(iy) / 2),
	xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 0)].push_back(i);
	m_frac_regular[el_old(ix, 0)].push_back(1.0);
	}
	}

	// left-top: narrow-to-wide
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e =
	xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));

	for (size_t k = 0; k < e.size() - 1; ++k) {
	m_elem_regular[el_old(ix, 0)].push_back(e(k));
	m_frac_regular[el_old(ix, 0)].push_back(1.0);
	}

	m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
	m_frac_regular[el_old(ix, 0)].push_back(0.5);
	}
	}

	// right
	{
	// right-bottom: regular
	{
	auto e = xt::view(
	block,
	xt::range(0, nhy(iy) / 2),
	xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));

	for (auto& i : e) {
	m_elem_regular[el_old(ix, 2)].push_back(i);
	m_frac_regular[el_old(ix, 2)].push_back(1.0);
	}
	}

	// right-top: narrow-to-wide
	for (size_t j = 0; j < nhy(iy) / 2; ++j) {
	auto e = xt::view(
	block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));

	m_elem_regular[el_old(ix, 2)].push_back(e(0));
	m_frac_regular[el_old(ix, 2)].push_back(0.5);

	for (size_t k = 1; k < e.size(); ++k) {
	m_elem_regular[el_old(ix, 2)].push_back(e(k));
	m_frac_regular[el_old(ix, 2)].push_back(1.0);
	}
	}
	}
	}
	}
	}
	}

	inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const
	{
	return m_regular;
	}

	inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const
	{
	return m_finelayer;
	}

	inline std::vector<std::vector<size_t>> FineLayer2Regular::getMap() const
	{
	return m_elem_regular;
	}

	inline std::vector<std::vector<double>> FineLayer2Regular::getMapFraction() const
	{
	return m_frac_regular;
	}

	inline xt::xtensor<double, 1>
	FineLayer2Regular::mapToRegular(const xt::xtensor<double, 1>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());

	xt::xtensor<double, 1> ret = xt::zeros<double>({m_regular.nelem()});

	for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
	for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
	ret(m_elem_regular[e][i]) += m_frac_regular[e][i] * data(e);
	}
	}

	return ret;
	}

	inline xt::xtensor<double, 2>
	FineLayer2Regular::mapToRegular(const xt::xtensor<double, 2>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());

	xt::xtensor<double, 2> ret = xt::zeros<double>({m_regular.nelem(), data.shape(1)});

	for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
	for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
	xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
	}
	}

	return ret;
	}

	inline xt::xtensor<double, 4>
	FineLayer2Regular::mapToRegular(const xt::xtensor<double, 4>& data) const
	{
	GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());

	xt::xtensor<double, 4> ret =
	xt::zeros<double>({m_regular.nelem(), data.shape(1), data.shape(2), data.shape(3)});

	for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
	for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
	xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
	}
	}

	return ret;
	}

	} // namespace Map

	} // namespace Quad4
	} // namespace Mesh
	} // namespace GooseFEM

	#endif
	diff --git a/python/MeshQuad4.hpp b/python/MeshQuad4.hpp
	index 9bebf32..b88643c 100644
	--- a/python/MeshQuad4.hpp
	+++ b/python/MeshQuad4.hpp
	@@ -1,263 +1,247 @@
	/* =================================================================================================

	(c - GPLv3) T.W.J. de Geus (Tom) \| tom@geus.me \| www.geus.me \| github.com/tdegeus/GooseFEM

	================================================================================================= */

	#include <GooseFEM/GooseFEM.h>
	#include <pybind11/pybind11.h>
	+#include <pybind11/stl.h>
	#include <pyxtensor/pyxtensor.hpp>

	namespace py = pybind11;

	void init_MeshQuad4(py::module& m)
	{

	py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")

	.def(
	py::init<size_t, size_t, double>(),
	"Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
	"size 'h'",
	py::arg("nx"),
	py::arg("ny"),
	py::arg("h") = 1.)

	.def("coor", &GooseFEM::Mesh::Quad4::Regular::coor)

	.def("conn", &GooseFEM::Mesh::Quad4::Regular::conn)

	.def("nelem", &GooseFEM::Mesh::Quad4::Regular::nelem)

	.def("nnode", &GooseFEM::Mesh::Quad4::Regular::nnode)

	.def("nne", &GooseFEM::Mesh::Quad4::Regular::nne)

	.def("ndim", &GooseFEM::Mesh::Quad4::Regular::ndim)

	.def("nelx", &GooseFEM::Mesh::Quad4::Regular::nelx)

	.def("nely", &GooseFEM::Mesh::Quad4::Regular::nely)

	.def("h", &GooseFEM::Mesh::Quad4::Regular::h)

	.def("getElementType", &GooseFEM::Mesh::Quad4::Regular::getElementType)

	.def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge)

	.def("nodesTopEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge)

	.def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge)

	.def("nodesRightEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge)

	.def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge)

	.def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge)

	.def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge)

	.def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge)

	.def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner)

	.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)

	.def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner)

	.def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner)

	.def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner)

	.def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner)

	.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)

	.def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner)

	.def("dofs", &GooseFEM::Mesh::Quad4::Regular::dofs)

	.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)

	.def("nodesOrigin", &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)

	.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)

	.def("elementgrid", &GooseFEM::Mesh::Quad4::Regular::elementgrid)

	.def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
	return "<GooseFEM.Mesh.Quad4.Regular>";
	});

	py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")

	.def(
	py::init<size_t, size_t, double, size_t>(),
	"FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical "
	"direction",
	py::arg("nx"),
	py::arg("ny"),
	py::arg("h") = 1.,
	py::arg("nfine") = 1)

	.def(
	py::init<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(),
	"Map connectivity to generating FineLayer-object.",
	py::arg("coor"),
	py::arg("conn"))

	.def("coor", &GooseFEM::Mesh::Quad4::FineLayer::coor)
	-
	.def("conn", &GooseFEM::Mesh::Quad4::FineLayer::conn)
	-
	.def("nelem", &GooseFEM::Mesh::Quad4::FineLayer::nelem)
	-
	.def("nnode", &GooseFEM::Mesh::Quad4::FineLayer::nnode)
	-
	.def("nne", &GooseFEM::Mesh::Quad4::FineLayer::nne)
	-
	.def("ndim", &GooseFEM::Mesh::Quad4::FineLayer::ndim)
	-
	.def("nelx", &GooseFEM::Mesh::Quad4::FineLayer::nelx)
	-
	.def("nely", &GooseFEM::Mesh::Quad4::FineLayer::nely)
	-
	.def("h", &GooseFEM::Mesh::Quad4::FineLayer::h)
	-
	.def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
	-
	.def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer)

	- .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
	+ .def(
	+ "elementgrid_ravel",
	+ &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_ravel,
	+ py::arg("rows_range"),
	+ py::arg("cols_range"))

	- .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
	+ .def(
	+ "elementgrid_around_ravel",
	+ &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_around_ravel,
	+ py::arg("element"),
	+ py::arg("size"),
	+ py::arg("periodic") = true)

	+ .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
	+ .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
	.def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge)
	-
	.def("nodesRightEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge)
	-
	.def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge)
	-
	.def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge)
	-
	.def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge)
	-
	.def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge)
	-
	.def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner)
	-
	.def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
	-
	.def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner)
	-
	.def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner)
	-
	.def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner)
	-
	.def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner)
	-
	.def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
	-
	.def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner)
	-
	.def("dofs", &GooseFEM::Mesh::Quad4::FineLayer::dofs)
	-
	.def("nodesPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
	-
	.def("nodesOrigin", &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
	-
	.def("dofsPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
	-
	.def("roll", &GooseFEM::Mesh::Quad4::FineLayer::roll)

	.def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
	return "<GooseFEM.Mesh.Quad4.FineLayer>";
	});
	}

	void init_MeshQuad4Map(py::module& m)
	{

	py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")

	.def(
	py::init<const GooseFEM::Mesh::Quad4::Regular&, size_t, size_t>(),
	"Refine a regular mesh",
	py::arg("mesh"),
	py::arg("nx"),
	py::arg("ny"))

	.def("getCoarseMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)

	.def("getFineMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)

	.def("getMap", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)

	.def(
	"mapToCoarse",
	py::overload_cast<const xt::xtensor<double, 1>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))

	.def(
	"mapToCoarse",
	py::overload_cast<const xt::xtensor<double, 2>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))

	.def(
	"mapToCoarse",
	py::overload_cast<const xt::xtensor<double, 4>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToCoarse, py::const_))

	.def(
	"mapToFine",
	py::overload_cast<const xt::xtensor<double, 1>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))

	.def(
	"mapToFine",
	py::overload_cast<const xt::xtensor<double, 2>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))

	.def(
	"mapToFine",
	py::overload_cast<const xt::xtensor<double, 4>&>(
	&GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine, py::const_))

	.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) {
	return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>";
	});

	py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")

	.def(
	py::init<const GooseFEM::Mesh::Quad4::FineLayer&>(),
	"Map a FineLayer mesh to a Regular mesh",
	py::arg("mesh"))

	.def("getRegularMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)

	.def("getFineLayerMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)

	.def("getMap", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)

	.def("getMapFraction", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)

	.def(
	"mapToRegular",
	py::overload_cast<const xt::xtensor<double, 1>&>(
	&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))

	.def(
	"mapToRegular",
	py::overload_cast<const xt::xtensor<double, 2>&>(
	&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))

	.def(
	"mapToRegular",
	py::overload_cast<const xt::xtensor<double, 4>&>(
	&GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular, py::const_))

	.def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) {
	return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>";
	});
	}
	diff --git a/test/basic/MeshQuad4.cpp b/test/basic/MeshQuad4.cpp
	index ae81cfc..aa9571f 100644
	--- a/test/basic/MeshQuad4.cpp
	+++ b/test/basic/MeshQuad4.cpp
	@@ -1,574 +1,657 @@

	#include <catch2/catch.hpp>
	#include <xtensor/xrandom.hpp>
	#include <xtensor/xmath.hpp>
	#include <GooseFEM/GooseFEM.h>

	#define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12));

	TEST_CASE("GooseFEM::MeshQuad4", "MeshQuad4.h")
	{

	SECTION("Regular")
	{
	xt::xtensor<double, 2> coor_ = {{0., 0.}, {1., 0.}, {2., 0.}, {3., 0.}, {4., 0.}, {5., 0.},
	{0., 1.}, {1., 1.}, {2., 1.}, {3., 1.}, {4., 1.}, {5., 1.},
	{0., 2.}, {1., 2.}, {2., 2.}, {3., 2.}, {4., 2.}, {5., 2.},
	{0., 3.}, {1., 3.}, {2., 3.}, {3., 3.}, {4., 3.}, {5., 3.}};

	xt::xtensor<double, 2> conn_ = {{0, 1, 7, 6},
	{1, 2, 8, 7},
	{2, 3, 9, 8},
	{3, 4, 10, 9},
	{4, 5, 11, 10},
	{6, 7, 13, 12},
	{7, 8, 14, 13},
	{8, 9, 15, 14},
	{9, 10, 16, 15},
	{10, 11, 17, 16},
	{12, 13, 19, 18},
	{13, 14, 20, 19},
	{14, 15, 21, 20},
	{15, 16, 22, 21},
	{16, 17, 23, 22}};

	size_t nelem_ = 15;
	size_t nnode_ = 24;
	size_t nne_ = 4;
	size_t ndim_ = 2;
	size_t nnodePeriodic_ = 15;

	xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3, 4, 5};
	xt::xtensor<size_t, 1> nodesTopEdge_ = {18, 19, 20, 21, 22, 23};
	xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 6, 12, 18};
	xt::xtensor<size_t, 1> nodesRightEdge_ = {5, 11, 17, 23};
	xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2, 3, 4};
	xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {19, 20, 21, 22};
	xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {6, 12};
	xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {11, 17};

	size_t nodesBottomLeftCorner_ = 0;
	size_t nodesBottomRightCorner_ = 5;
	size_t nodesTopLeftCorner_ = 18;
	size_t nodesTopRightCorner_ = 23;
	size_t nodesLeftBottomCorner_ = 0;
	size_t nodesLeftTopCorner_ = 18;
	size_t nodesRightBottomCorner_ = 5;
	size_t nodesRightTopCorner_ = 23;

	xt::xtensor<size_t, 2> dofs_ = {{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}};

	xt::xtensor<size_t, 2> nodesPeriodic_ = {
	{0, 5}, {0, 23}, {0, 18}, {1, 19}, {2, 20}, {3, 21}, {4, 22}, {6, 11}, {12, 17}};

	size_t nodesOrigin_ = 0;

	xt::xtensor<size_t, 2> dofsPeriodic_ = {
	{0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}, {10, 11}, {12, 13},
	{14, 15}, {16, 17}, {18, 19}, {10, 11}, {20, 21}, {22, 23}, {24, 25}, {26, 27},
	{28, 29}, {20, 21}, {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {0, 1}};

	GooseFEM::Mesh::Quad4::Regular mesh(5, 3);

	xt::xtensor<double, 2> coor = mesh.coor();
	xt::xtensor<size_t, 2> conn = mesh.conn();

	size_t nelem = mesh.nelem();
	size_t nnode = mesh.nnode();
	size_t nne = mesh.nne();
	size_t ndim = mesh.ndim();
	size_t nnodePeriodic = mesh.nnode() - mesh.nodesPeriodic().shape(0);

	xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
	xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
	xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
	xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
	xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
	xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
	xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
	xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();

	size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
	size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
	size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
	size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
	size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
	size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
	size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
	size_t nodesRightTopCorner = mesh.nodesRightTopCorner();

	xt::xtensor<size_t, 2> dofs = mesh.dofs();

	xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
	size_t nodesOrigin = mesh.nodesOrigin();
	xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();

	REQUIRE(nelem == nelem_);
	REQUIRE(nnode == nnode_);
	REQUIRE(nne == nne_);
	REQUIRE(ndim == ndim_);
	REQUIRE(nnodePeriodic == nnodePeriodic_);
	REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
	REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
	REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
	REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
	REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
	REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
	REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
	REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
	REQUIRE(nodesOrigin == nodesOrigin_);

	REQUIRE(xt::allclose(coor, coor_));

	REQUIRE(xt::all(xt::equal(conn, conn_)));
	REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
	REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
	REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
	REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
	REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
	REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
	}

	SECTION("FineLayer")
	{
	xt::xtensor<double, 2> coor_ = {
	{0., 0.}, {3., 0.}, {6., 0.}, {9., 0.}, {0., 3.}, {3., 3.}, {6., 3.}, {9., 3.},
	{0., 6.}, {3., 6.}, {6., 6.}, {9., 6.}, {1., 7.}, {2., 7.}, {4., 7.}, {5., 7.},
	{7., 7.}, {8., 7.}, {0., 8.}, {1., 8.}, {2., 8.}, {3., 8.}, {4., 8.}, {5., 8.},
	{6., 8.}, {7., 8.}, {8., 8.}, {9., 8.}, {0., 9.}, {1., 9.}, {2., 9.}, {3., 9.},
	{4., 9.}, {5., 9.}, {6., 9.}, {7., 9.}, {8., 9.}, {9., 9.}, {0., 10.}, {1., 10.},
	{2., 10.}, {3., 10.}, {4., 10.}, {5., 10.}, {6., 10.}, {7., 10.}, {8., 10.}, {9., 10.},
	{0., 11.}, {1., 11.}, {2., 11.}, {3., 11.}, {4., 11.}, {5., 11.}, {6., 11.}, {7., 11.},
	{8., 11.}, {9., 11.}, {0., 12.}, {1., 12.}, {2., 12.}, {3., 12.}, {4., 12.}, {5., 12.},
	{6., 12.}, {7., 12.}, {8., 12.}, {9., 12.}, {0., 13.}, {1., 13.}, {2., 13.}, {3., 13.},
	{4., 13.}, {5., 13.}, {6., 13.}, {7., 13.}, {8., 13.}, {9., 13.}, {1., 14.}, {2., 14.},
	{4., 14.}, {5., 14.}, {7., 14.}, {8., 14.}, {0., 15.}, {3., 15.}, {6., 15.}, {9., 15.},
	{0., 18.}, {3., 18.}, {6., 18.}, {9., 18.}, {0., 21.}, {3., 21.}, {6., 21.}, {9., 21.}};

	xt::xtensor<double, 2> conn_ = {
	{0, 1, 5, 4}, {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 9, 8},
	{5, 6, 10, 9}, {6, 7, 11, 10}, {8, 9, 13, 12}, {9, 21, 20, 13},
	{12, 13, 20, 19}, {8, 12, 19, 18}, {9, 10, 15, 14}, {10, 24, 23, 15},
	{14, 15, 23, 22}, {9, 14, 22, 21}, {10, 11, 17, 16}, {11, 27, 26, 17},
	{16, 17, 26, 25}, {10, 16, 25, 24}, {18, 19, 29, 28}, {19, 20, 30, 29},
	{20, 21, 31, 30}, {21, 22, 32, 31}, {22, 23, 33, 32}, {23, 24, 34, 33},
	{24, 25, 35, 34}, {25, 26, 36, 35}, {26, 27, 37, 36}, {28, 29, 39, 38},
	{29, 30, 40, 39}, {30, 31, 41, 40}, {31, 32, 42, 41}, {32, 33, 43, 42},
	{33, 34, 44, 43}, {34, 35, 45, 44}, {35, 36, 46, 45}, {36, 37, 47, 46},
	{38, 39, 49, 48}, {39, 40, 50, 49}, {40, 41, 51, 50}, {41, 42, 52, 51},
	{42, 43, 53, 52}, {43, 44, 54, 53}, {44, 45, 55, 54}, {45, 46, 56, 55},
	{46, 47, 57, 56}, {48, 49, 59, 58}, {49, 50, 60, 59}, {50, 51, 61, 60},
	{51, 52, 62, 61}, {52, 53, 63, 62}, {53, 54, 64, 63}, {54, 55, 65, 64},
	{55, 56, 66, 65}, {56, 57, 67, 66}, {58, 59, 69, 68}, {59, 60, 70, 69},
	{60, 61, 71, 70}, {61, 62, 72, 71}, {62, 63, 73, 72}, {63, 64, 74, 73},
	{64, 65, 75, 74}, {65, 66, 76, 75}, {66, 67, 77, 76}, {68, 69, 78, 84},
	{69, 70, 79, 78}, {70, 71, 85, 79}, {78, 79, 85, 84}, {71, 72, 80, 85},
	{72, 73, 81, 80}, {73, 74, 86, 81}, {80, 81, 86, 85}, {74, 75, 82, 86},
	{75, 76, 83, 82}, {76, 77, 87, 83}, {82, 83, 87, 86}, {84, 85, 89, 88},
	{85, 86, 90, 89}, {86, 87, 91, 90}, {88, 89, 93, 92}, {89, 90, 94, 93},
	{90, 91, 95, 94}};

	size_t nelem_ = 81;
	size_t nnode_ = 96;
	size_t nne_ = 4;
	size_t ndim_ = 2;
	size_t shape_x_ = 9;
	size_t shape_y_ = 21;

	xt::xtensor<size_t, 1> elementsMiddleLayer_ = {36, 37, 38, 39, 40, 41, 42, 43, 44};

	xt::xtensor<size_t, 1> nodesBottomEdge_ = {0, 1, 2, 3};
	xt::xtensor<size_t, 1> nodesTopEdge_ = {92, 93, 94, 95};
	xt::xtensor<size_t, 1> nodesLeftEdge_ = {0, 4, 8, 18, 28, 38, 48, 58, 68, 84, 88, 92};
	xt::xtensor<size_t, 1> nodesRightEdge_ = {3, 7, 11, 27, 37, 47, 57, 67, 77, 87, 91, 95};
	xt::xtensor<size_t, 1> nodesBottomOpenEdge_ = {1, 2};
	xt::xtensor<size_t, 1> nodesTopOpenEdge_ = {93, 94};
	xt::xtensor<size_t, 1> nodesLeftOpenEdge_ = {4, 8, 18, 28, 38, 48, 58, 68, 84, 88};
	xt::xtensor<size_t, 1> nodesRightOpenEdge_ = {7, 11, 27, 37, 47, 57, 67, 77, 87, 91};

	size_t nodesBottomLeftCorner_ = 0;
	size_t nodesBottomRightCorner_ = 3;
	size_t nodesTopLeftCorner_ = 92;
	size_t nodesTopRightCorner_ = 95;
	size_t nodesLeftBottomCorner_ = 0;
	size_t nodesLeftTopCorner_ = 92;
	size_t nodesRightBottomCorner_ = 3;
	size_t nodesRightTopCorner_ = 95;

	xt::xtensor<size_t, 2> dofs_ = {
	{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}};

	xt::xtensor<size_t, 2> nodesPeriodic_ = {{0, 3},
	{0, 95},
	{0, 92},
	{1, 93},
	{2, 94},
	{4, 7},
	{8, 11},
	{18, 27},
	{28, 37},
	{38, 47},
	{48, 57},
	{58, 67},
	{68, 77},
	{84, 87},
	{88, 91}};

	size_t nodesOrigin_ = 0;

	xt::xtensor<size_t, 2> dofsPeriodic_ = {
	{0, 1}, {2, 3}, {4, 5}, {0, 1}, {6, 7}, {8, 9}, {10, 11},
	{6, 7}, {12, 13}, {14, 15}, {16, 17}, {12, 13}, {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}, {30, 31},
	{48, 49}, {50, 51}, {52, 53}, {54, 55}, {56, 57}, {58, 59}, {60, 61},
	{62, 63}, {64, 65}, {48, 49}, {66, 67}, {68, 69}, {70, 71}, {72, 73},
	{74, 75}, {76, 77}, {78, 79}, {80, 81}, {82, 83}, {66, 67}, {84, 85},
	{86, 87}, {88, 89}, {90, 91}, {92, 93}, {94, 95}, {96, 97}, {98, 99},
	{100, 101}, {84, 85}, {102, 103}, {104, 105}, {106, 107}, {108, 109}, {110, 111},
	{112, 113}, {114, 115}, {116, 117}, {118, 119}, {102, 103}, {120, 121}, {122, 123},
	{124, 125}, {126, 127}, {128, 129}, {130, 131}, {132, 133}, {134, 135}, {136, 137},
	{120, 121}, {138, 139}, {140, 141}, {142, 143}, {144, 145}, {146, 147}, {148, 149},
	{150, 151}, {152, 153}, {154, 155}, {150, 151}, {156, 157}, {158, 159}, {160, 161},
	{156, 157}, {0, 1}, {2, 3}, {4, 5}, {0, 1}};

	GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18);

	xt::xtensor<double, 2> coor = mesh.coor();
	xt::xtensor<size_t, 2> conn = mesh.conn();

	size_t nelem = mesh.nelem();
	size_t nnode = mesh.nnode();
	size_t nne = mesh.nne();
	size_t ndim = mesh.ndim();
	size_t shape_x = mesh.nelx();
	size_t shape_y = mesh.nely();

	xt::xtensor<size_t, 1> elementsMiddleLayer = mesh.elementsMiddleLayer();

	xt::xtensor<size_t, 1> nodesBottomEdge = mesh.nodesBottomEdge();
	xt::xtensor<size_t, 1> nodesTopEdge = mesh.nodesTopEdge();
	xt::xtensor<size_t, 1> nodesLeftEdge = mesh.nodesLeftEdge();
	xt::xtensor<size_t, 1> nodesRightEdge = mesh.nodesRightEdge();
	xt::xtensor<size_t, 1> nodesBottomOpenEdge = mesh.nodesBottomOpenEdge();
	xt::xtensor<size_t, 1> nodesTopOpenEdge = mesh.nodesTopOpenEdge();
	xt::xtensor<size_t, 1> nodesLeftOpenEdge = mesh.nodesLeftOpenEdge();
	xt::xtensor<size_t, 1> nodesRightOpenEdge = mesh.nodesRightOpenEdge();

	size_t nodesBottomLeftCorner = mesh.nodesBottomLeftCorner();
	size_t nodesBottomRightCorner = mesh.nodesBottomRightCorner();
	size_t nodesTopLeftCorner = mesh.nodesTopLeftCorner();
	size_t nodesTopRightCorner = mesh.nodesTopRightCorner();
	size_t nodesLeftBottomCorner = mesh.nodesLeftBottomCorner();
	size_t nodesLeftTopCorner = mesh.nodesLeftTopCorner();
	size_t nodesRightBottomCorner = mesh.nodesRightBottomCorner();
	size_t nodesRightTopCorner = mesh.nodesRightTopCorner();

	xt::xtensor<size_t, 2> dofs = mesh.dofs();

	xt::xtensor<size_t, 2> nodesPeriodic = mesh.nodesPeriodic();
	size_t nodesOrigin = mesh.nodesOrigin();
	xt::xtensor<size_t, 2> dofsPeriodic = mesh.dofsPeriodic();

	REQUIRE(nelem == nelem_);
	REQUIRE(nnode == nnode_);
	REQUIRE(nne == nne_);
	REQUIRE(ndim == ndim_);
	REQUIRE(shape_x == shape_x_);
	REQUIRE(shape_y == shape_y_);
	REQUIRE(nodesBottomLeftCorner == nodesBottomLeftCorner_);
	REQUIRE(nodesBottomRightCorner == nodesBottomRightCorner_);
	REQUIRE(nodesTopLeftCorner == nodesTopLeftCorner_);
	REQUIRE(nodesTopRightCorner == nodesTopRightCorner_);
	REQUIRE(nodesLeftBottomCorner == nodesLeftBottomCorner_);
	REQUIRE(nodesLeftTopCorner == nodesLeftTopCorner_);
	REQUIRE(nodesRightBottomCorner == nodesRightBottomCorner_);
	REQUIRE(nodesRightTopCorner == nodesRightTopCorner_);
	REQUIRE(nodesOrigin == nodesOrigin_);

	REQUIRE(xt::allclose(coor, coor_));

	REQUIRE(xt::all(xt::equal(elementsMiddleLayer, elementsMiddleLayer_)));
	REQUIRE(xt::all(xt::equal(conn, conn_)));
	REQUIRE(xt::all(xt::equal(nodesBottomEdge, nodesBottomEdge_)));
	REQUIRE(xt::all(xt::equal(nodesTopEdge, nodesTopEdge_)));
	REQUIRE(xt::all(xt::equal(nodesLeftEdge, nodesLeftEdge_)));
	REQUIRE(xt::all(xt::equal(nodesRightEdge, nodesRightEdge_)));
	REQUIRE(xt::all(xt::equal(nodesBottomOpenEdge, nodesBottomOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesTopOpenEdge, nodesTopOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesLeftOpenEdge, nodesLeftOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesRightOpenEdge, nodesRightOpenEdge_)));
	REQUIRE(xt::all(xt::equal(nodesPeriodic, nodesPeriodic_)));
	REQUIRE(xt::all(xt::equal(dofsPeriodic, dofsPeriodic_)));
	}

	SECTION("FineLayer::elementgrid_ravel - uniform")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
	xt::xtensor<size_t, 1> a = {
	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};
	REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({0, 5}, {0, 5}))));
	+ REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({0, 5}, {}))));
	+ REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({}, {0, 5}))));
	+ REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({}, {}))));

	xt::xtensor<size_t, 1> b = {
	5, 6, 7, 8, 9,
	10, 11, 12, 13, 14,
	15, 16, 17, 18, 19};
	REQUIRE(xt::all(xt::equal(b, mesh.elementgrid_ravel({1, 4}, {0, 5}))));

	xt::xtensor<size_t, 1> c = {
	6, 7, 8,
	11, 12, 13,
	16, 17, 18};
	REQUIRE(xt::all(xt::equal(c, mesh.elementgrid_ravel({1, 4}, {1, 4}))));
	}

	SECTION("FineLayer::elementgrid_ravel - refined")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(6, 18);
	xt::xtensor<size_t, 1> a = {
	19, 20, 21, 22,
	25, 26, 27, 28,
	31, 32, 33, 34};
	REQUIRE(xt::all(xt::equal(a, mesh.elementgrid_ravel({9, 12}, {1, 5}))));

	xt::xtensor<size_t, 1> b = {
	0, 1,
	2, 3};
	REQUIRE(xt::all(xt::equal(b, mesh.elementgrid_ravel({0, 6}, {0, 6}))));

	xt::xtensor<size_t, 1> c = {
	50, 51,
	52, 53};
	REQUIRE(xt::all(xt::equal(c, mesh.elementgrid_ravel({15, 21}, {0, 6}))));

	xt::xtensor<size_t, 1> d = {
	4, 5, 6, 7, 8, 9, 10, 11};
	REQUIRE(xt::all(xt::equal(d, mesh.elementgrid_ravel({6, 8}, {0, 6}))));
	}

	+ SECTION("FineLayer::elementgrid_ravel - uniform")
	+ {
	+ GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);
	+ xt::xtensor<size_t, 1> a = {
	+ 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};
	+
	+ REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(10, 2)))));
	+ REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(11, 2)))));
	+ REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(12, 2)))));
	+ REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(13, 2)))));
	+ REQUIRE(xt::all(xt::equal(a, xt::sort(mesh.elementgrid_around_ravel(14, 2)))));
	+
	+ xt::xtensor<size_t, 1> b10 = {
	+ 5, 6, 9,
	+ 10, 11, 14,
	+ 15, 16, 19};
	+
	+ xt::xtensor<size_t, 1> b11 = {
	+ 5, 6, 7,
	+ 10, 11, 12,
	+ 15, 16, 17};
	+
	+ xt::xtensor<size_t, 1> b12 = {
	+ 6, 7, 8,
	+ 11, 12, 13,
	+ 16, 17, 18};
	+
	+ xt::xtensor<size_t, 1> b13 = {
	+ 7, 8, 9,
	+ 12, 13, 14,
	+ 17, 18, 19};
	+
	+ xt::xtensor<size_t, 1> b14 = {
	+ 5, 8, 9,
	+ 10, 13, 14,
	+ 15, 18, 19};
	+
	+ REQUIRE(xt::all(xt::equal(xt::sort(b10), xt::sort(mesh.elementgrid_around_ravel(10, 1)))));
	+ REQUIRE(xt::all(xt::equal(xt::sort(b11), xt::sort(mesh.elementgrid_around_ravel(11, 1)))));
	+ REQUIRE(xt::all(xt::equal(xt::sort(b12), xt::sort(mesh.elementgrid_around_ravel(12, 1)))));
	+ REQUIRE(xt::all(xt::equal(xt::sort(b13), xt::sort(mesh.elementgrid_around_ravel(13, 1)))));
	+ REQUIRE(xt::all(xt::equal(xt::sort(b14), xt::sort(mesh.elementgrid_around_ravel(14, 1)))));
	+
	+ REQUIRE(xt::all(xt::equal(10, xt::sort(mesh.elementgrid_around_ravel(10, 0)))));
	+ REQUIRE(xt::all(xt::equal(11, xt::sort(mesh.elementgrid_around_ravel(11, 0)))));
	+ REQUIRE(xt::all(xt::equal(12, xt::sort(mesh.elementgrid_around_ravel(12, 0)))));
	+ REQUIRE(xt::all(xt::equal(13, xt::sort(mesh.elementgrid_around_ravel(13, 0)))));
	+ REQUIRE(xt::all(xt::equal(14, xt::sort(mesh.elementgrid_around_ravel(14, 0)))));
	+ }
	+
	+ SECTION("FineLayer::elementgrid_around_ravel")
	+ {
	+ GooseFEM::Mesh::Quad4::FineLayer mesh(12, 18);
	+ xt::xtensor<size_t, 1> r0 = {
	+ 36, 37, 71,
	+ 48, 49, 59,
	+ 60, 61, 47,
	+ };
	+ xt::xtensor<size_t, 1> r1 = {
	+ 36, 37, 38,
	+ 48, 49, 50,
	+ 60, 61, 62,
	+ };
	+ xt::xtensor<size_t, 1> r12 = {
	+ 46, 47, 36,
	+ 58, 59, 48,
	+ 70, 71, 60,
	+ };
	+
	+ REQUIRE(xt::all(xt::equal(xt::sort(r0), xt::sort(mesh.elementgrid_around_ravel(48, 1)))));
	+ for (size_t n = 0; n < 10; ++n) {
	+ REQUIRE(xt::all(xt::equal(xt::sort(r1) + n, xt::sort(mesh.elementgrid_around_ravel(49 + n, 1)))));
	+ }
	+ REQUIRE(xt::all(xt::equal(xt::sort(r12), xt::sort(mesh.elementgrid_around_ravel(59, 1)))));
	+ }
	+
	SECTION("FineLayer - replica - trivial")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(1, 1);
	GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
	REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
	REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
	}

	SECTION("FineLayer - replica - equi-sized")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(4, 4);
	GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
	REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
	REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
	}

	SECTION("FineLayer - replica")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(27, 27);
	GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
	REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
	REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
	}

	SECTION("FineLayer - replica - one")
	{
	double h = xt::numeric_constants<double>::PI;
	GooseFEM::Mesh::Quad4::FineLayer mesh(1, 1, h);
	GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
	REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
	REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
	}

	SECTION("FineLayer - replica - large")
	{
	size_t N = 2 * std::pow(3, 6);
	double h = xt::numeric_constants<double>::PI;
	GooseFEM::Mesh::Quad4::FineLayer mesh(N, N, h);
	GooseFEM::Mesh::Quad4::FineLayer replica(mesh.coor(), mesh.conn());
	REQUIRE(xt::all(xt::equal(mesh.conn(), replica.conn())));
	REQUIRE(xt::allclose(mesh.coor(), replica.coor()));
	}

	SECTION("FineLayer - roll - trivial")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
	xt::xtensor<size_t, 1> m0 = {
	0, 1, 2,
	3, 4, 5,
	6, 7, 8
	};
	xt::xtensor<size_t, 1> m1 = {
	2, 0, 1,
	5, 3, 4,
	8, 6, 7
	};
	xt::xtensor<size_t, 1> m2 = {
	1, 2, 0,
	4, 5, 3,
	7, 8, 6,
	};
	REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
	REQUIRE(xt::all(xt::equal(mesh.roll(2), m2)));
	REQUIRE(xt::all(xt::equal(mesh.roll(3), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(4), m1)));
	REQUIRE(xt::all(xt::equal(mesh.roll(5), m2)));
	}

	SECTION("FineLayer - roll - a")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(6, 18);
	xt::xtensor<size_t, 1> m0 = {
	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
	};
	xt::xtensor<size_t, 1> m1 = {
	1, 0,
	3, 2,
	8, 9, 10, 11, 4, 5, 6, 7,
	15, 16, 17, 12, 13, 14,
	21, 22, 23, 18, 19, 20,
	27, 28, 29, 24, 25, 26,
	33, 34, 35, 30, 31, 32,
	39, 40, 41, 36, 37, 38,
	46, 47, 48, 49, 42, 43, 44, 45,
	51, 50,
	53, 52
	};
	REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
	REQUIRE(xt::all(xt::equal(mesh.roll(2), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(3), m1)));
	}

	SECTION("FineLayer - roll - b")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(9, 18);
	xt::xtensor<size_t, 1> m0 = {
	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
	};
	xt::xtensor<size_t, 1> m1 = {
	2, 0, 1,
	5, 3, 4,
	14, 15, 16, 17, 6, 7, 8, 9, 10, 11, 12, 13,
	24, 25, 26, 18, 19, 20, 21, 22, 23,
	33, 34, 35, 27, 28, 29, 30, 31, 32,
	42, 43, 44, 36, 37, 38, 39, 40, 41,
	51, 52, 53, 45, 46, 47, 48, 49, 50,
	60, 61, 62, 54, 55, 56, 57, 58, 59,
	71, 72, 73, 74, 63, 64, 65, 66, 67, 68, 69, 70,
	77, 75, 76,
	80, 78, 79
	};
	xt::xtensor<size_t, 1> m2 = {
	1, 2, 0,
	4, 5, 3,
	10, 11, 12, 13, 14, 15, 16, 17, 6, 7, 8, 9,
	21, 22, 23, 24, 25, 26, 18, 19, 20,
	30, 31, 32, 33, 34, 35, 27, 28, 29,
	39, 40, 41, 42, 43, 44, 36, 37, 38,
	48, 49, 50, 51, 52, 53, 45, 46, 47,
	57, 58, 59, 60, 61, 62, 54, 55, 56,
	67, 68, 69, 70, 71, 72, 73, 74, 63, 64, 65, 66,
	76, 77, 75,
	79, 80, 78
	};
	REQUIRE(xt::all(xt::equal(mesh.roll(0), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(1), m1)));
	REQUIRE(xt::all(xt::equal(mesh.roll(2), m2)));
	REQUIRE(xt::all(xt::equal(mesh.roll(3), m0)));
	REQUIRE(xt::all(xt::equal(mesh.roll(4), m1)));
	REQUIRE(xt::all(xt::equal(mesh.roll(5), m2)));
	}

	SECTION("Map::RefineRegular")
	{
	GooseFEM::Mesh::Quad4::Regular mesh(5, 4);

	GooseFEM::Mesh::Quad4::Map::RefineRegular refine(mesh, 5, 3);

	xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
	auto a_ = refine.mapToCoarse(refine.mapToFine(a));

	REQUIRE(xt::allclose(a, xt::mean(a_, {1})));

	xt::xtensor<double, 2> b =
	xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
	auto b_ = refine.mapToCoarse(refine.mapToFine(b));

	REQUIRE(xt::allclose(xt::mean(b, {1}), xt::mean(b_, {1})));

	xt::xtensor<double, 4> c =
	xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
	auto c_ = refine.mapToCoarse(refine.mapToFine(c));

	REQUIRE(xt::allclose(xt::mean(c, {1}), xt::mean(c_, {1})));
	}

	SECTION("Map::FineLayer2Regular")
	{
	GooseFEM::Mesh::Quad4::FineLayer mesh(5, 5);

	GooseFEM::Mesh::Quad4::Map::FineLayer2Regular map(mesh);

	xt::xtensor<double, 1> a = xt::random::rand<double>({mesh.nelem()});
	auto a_ = map.mapToRegular(a);

	REQUIRE(xt::allclose(a, a_));

	xt::xtensor<double, 2> b =
	xt::random::rand<double>(std::array<size_t, 2>{mesh.nelem(), 4ul});
	auto b_ = map.mapToRegular(b);

	REQUIRE(xt::allclose(b, b_));

	xt::xtensor<double, 4> c =
	xt::random::rand<double>(std::array<size_t, 4>{mesh.nelem(), 4ul, 3ul, 3ul});
	auto c_ = map.mapToRegular(c);

	REQUIRE(xt::allclose(c, c_));
	}
	}

No OneTemporaryActions

File Metadata

View Options

Event Timeline

No OneTemporary
Actions