Page MenuHomec4science

No OneTemporary

File Metadata

Created
Fri, Mar 29, 07:40
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_));
}
}

Event Timeline