Page MenuHomec4science
Files F95269353

No OneTemporary
Actions

File Metadata

Created
Sat, Dec 14, 05:47

View Options

diff --git a/docs/Vector.rst b/docs/Vector.rst
index a1d5ac7..db6a66e 100644
--- a/docs/Vector.rst
+++ b/docs/Vector.rst
@@ -1,7 +1,12 @@
.. _vector:
****************
GooseFEM::Vector
****************
+The vector class can be used to switch between three representations:
+
+* "dofval" ``[ndof]``: degrees-of-freedom.
+* "nodevec" ``[nnode, ndim]``: vectors per node.
+* "elemvec" ``[nelem, nne, ndim]``: vector per node for each element.
diff --git a/docs/figures/Vector/mesh.idraw b/docs/figures/Vector/mesh.idraw
new file mode 100644
index 0000000..c4e7e73
Binary files /dev/null and b/docs/figures/Vector/mesh.idraw differ
diff --git a/docs/figures/Vector/mesh.svg b/docs/figures/Vector/mesh.svg
new file mode 100644
index 0000000..380948f
--- /dev/null
+++ b/docs/figures/Vector/mesh.svg
@@ -0,0 +1,68 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0" y="0" width="207.448" height="217.375" viewBox="0, 0, 207.448, 217.375">
+ <g id="mesh">
+ <path d="M16.103,192.346 L204.734,192.346" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M16.103,98.03 L204.734,98.03" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M16.103,3.714 L204.734,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M204.734,192.346 L204.734,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M110.419,192.346 L110.419,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M16.103,192.346 L16.103,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M15.603,195.059 C13.828,195.059 12.389,193.621 12.389,191.846 C12.389,190.071 13.828,188.632 15.603,188.632 C17.378,188.632 18.816,190.071 18.816,191.846 C18.816,193.621 17.378,195.059 15.603,195.059 z" fill="#000000"/>
+ <path d="M15.603,100.744 C13.828,100.744 12.389,99.305 12.389,97.53 C12.389,95.755 13.828,94.316 15.603,94.316 C17.378,94.316 18.816,95.755 18.816,97.53 C18.816,99.305 17.378,100.744 15.603,100.744 z" fill="#000000"/>
+ <path d="M15.603,6.428 C13.828,6.428 12.389,4.989 12.389,3.214 C12.389,1.439 13.828,0 15.603,0 C17.378,0 18.816,1.439 18.816,3.214 C18.816,4.989 17.378,6.428 15.603,6.428 z" fill="#000000"/>
+ <path d="M109.919,195.059 C108.144,195.059 106.705,193.621 106.705,191.846 C106.705,190.071 108.144,188.632 109.919,188.632 C111.693,188.632 113.132,190.071 113.132,191.846 C113.132,193.621 111.693,195.059 109.919,195.059 z" fill="#000000"/>
+ <path d="M109.919,100.744 C108.144,100.744 106.705,99.305 106.705,97.53 C106.705,95.755 108.144,94.316 109.919,94.316 C111.693,94.316 113.132,95.755 113.132,97.53 C113.132,99.305 111.693,100.744 109.919,100.744 z" fill="#000000"/>
+ <path d="M109.919,6.428 C108.144,6.428 106.705,4.989 106.705,3.214 C106.705,1.439 108.144,0 109.919,0 C111.693,0 113.132,1.439 113.132,3.214 C113.132,4.989 111.693,6.428 109.919,6.428 z" fill="#000000"/>
+ <path d="M204.234,195.059 C202.459,195.059 201.021,193.621 201.021,191.846 C201.021,190.071 202.459,188.632 204.234,188.632 C206.009,188.632 207.448,190.071 207.448,191.846 C207.448,193.621 206.009,195.059 204.234,195.059 z" fill="#000000"/>
+ <path d="M204.234,100.744 C202.459,100.744 201.021,99.305 201.021,97.53 C201.021,95.755 202.459,94.316 204.234,94.316 C206.009,94.316 207.448,95.755 207.448,97.53 C207.448,99.305 206.009,100.744 204.234,100.744 z" fill="#000000"/>
+ <path d="M204.234,6.428 C202.459,6.428 201.021,4.989 201.021,3.214 C201.021,1.439 202.459,0 204.234,0 C206.009,0 207.448,1.439 207.448,3.214 C207.448,4.989 206.009,6.428 204.234,6.428 z" fill="#000000"/>
+ </g>
+ <g id="node_numbers">
+ <text transform="matrix(1, 0, 0, 1, 4.448, 205.875)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">0</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 96.911, 205.875)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">1</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 195.338, 205.875)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">2</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 4.448, 111.334)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">3</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 96.911, 111.334)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">4</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 195.338, 111.334)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">5</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 4.448, 19.635)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">6</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 96.911, 19.635)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">7</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 195.338, 19.635)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">8</tspan>
+ </text>
+ </g>
+ <g id="element_numbers">
+ <text transform="matrix(1, 0, 0, 1, 62.761, 144.688)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">0</tspan>
+ </text>
+ <path d="M63.261,155.69 C57.46,155.69 52.758,150.988 52.758,145.188 C52.758,139.387 57.46,134.686 63.261,134.686 C69.061,134.686 73.763,139.387 73.763,145.188 C73.763,150.988 69.061,155.69 63.261,155.69 z" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <text transform="matrix(1, 0, 0, 1, 157.076, 144.688)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">1</tspan>
+ </text>
+ <path d="M157.576,155.69 C151.776,155.69 147.074,150.988 147.074,145.188 C147.074,139.387 151.776,134.686 157.576,134.686 C163.377,134.686 168.079,139.387 168.079,145.188 C168.079,150.988 163.377,155.69 157.576,155.69 z" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <text transform="matrix(1, 0, 0, 1, 62.761, 50.372)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">2</tspan>
+ </text>
+ <path d="M63.261,61.374 C57.46,61.374 52.758,56.672 52.758,50.872 C52.758,45.072 57.46,40.37 63.261,40.37 C69.061,40.37 73.763,45.072 73.763,50.872 C73.763,56.672 69.061,61.374 63.261,61.374 z" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <text transform="matrix(1, 0, 0, 1, 157.076, 50.372)">
+ <tspan x="-4.448" y="4.5" font-family="Avenir-Book" font-size="16" fill="#000000">3</tspan>
+ </text>
+ <path d="M157.576,61.374 C151.776,61.374 147.074,56.672 147.074,50.872 C147.074,45.072 151.776,40.37 157.576,40.37 C163.377,40.37 168.079,45.072 168.079,50.872 C168.079,56.672 163.377,61.374 157.576,61.374 z" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ </g>
+</svg>
diff --git a/docs/figures/Vector/mesh_dofs-a.idraw b/docs/figures/Vector/mesh_dofs-a.idraw
new file mode 100644
index 0000000..b63db3a
Binary files /dev/null and b/docs/figures/Vector/mesh_dofs-a.idraw differ
diff --git a/docs/figures/Vector/mesh_dofs-a.svg b/docs/figures/Vector/mesh_dofs-a.svg
new file mode 100644
index 0000000..68799af
--- /dev/null
+++ b/docs/figures/Vector/mesh_dofs-a.svg
@@ -0,0 +1,105 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0" y="0" width="423.479" height="218.059" viewBox="0, 0, 423.479, 218.059">
+ <g id="mesh">
+ <path d="M232.133,192.346 L420.765,192.346" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,98.03 L420.765,98.03" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,3.714 L420.765,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M420.765,192.346 L420.765,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M326.449,192.346 L326.449,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,192.346 L232.133,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M231.634,195.059 C229.859,195.059 228.42,193.621 228.42,191.846 C228.42,190.071 229.859,188.632 231.634,188.632 C233.408,188.632 234.847,190.071 234.847,191.846 C234.847,193.621 233.408,195.059 231.634,195.059 z" fill="#000000"/>
+ <path d="M231.634,100.744 C229.859,100.744 228.42,99.305 228.42,97.53 C228.42,95.755 229.859,94.316 231.634,94.316 C233.408,94.316 234.847,95.755 234.847,97.53 C234.847,99.305 233.408,100.744 231.634,100.744 z" fill="#000000"/>
+ <path d="M231.634,6.428 C229.859,6.428 228.42,4.989 228.42,3.214 C228.42,1.439 229.859,0 231.634,0 C233.408,0 234.847,1.439 234.847,3.214 C234.847,4.989 233.408,6.428 231.634,6.428 z" fill="#000000"/>
+ <path d="M325.949,195.059 C324.174,195.059 322.736,193.621 322.736,191.846 C322.736,190.071 324.174,188.632 325.949,188.632 C327.724,188.632 329.163,190.071 329.163,191.846 C329.163,193.621 327.724,195.059 325.949,195.059 z" fill="#000000"/>
+ <path d="M325.949,100.744 C324.174,100.744 322.736,99.305 322.736,97.53 C322.736,95.755 324.174,94.316 325.949,94.316 C327.724,94.316 329.163,95.755 329.163,97.53 C329.163,99.305 327.724,100.744 325.949,100.744 z" fill="#000000"/>
+ <path d="M325.949,6.428 C324.174,6.428 322.736,4.989 322.736,3.214 C322.736,1.439 324.174,0 325.949,0 C327.724,0 329.163,1.439 329.163,3.214 C329.163,4.989 327.724,6.428 325.949,6.428 z" fill="#000000"/>
+ <path d="M420.265,195.059 C418.49,195.059 417.051,193.621 417.051,191.846 C417.051,190.071 418.49,188.632 420.265,188.632 C422.04,188.632 423.479,190.071 423.479,191.846 C423.479,193.621 422.04,195.059 420.265,195.059 z" fill="#000000"/>
+ <path d="M420.265,100.744 C418.49,100.744 417.051,99.305 417.051,97.53 C417.051,95.755 418.49,94.316 420.265,94.316 C422.04,94.316 423.479,95.755 423.479,97.53 C423.479,99.305 422.04,100.744 420.265,100.744 z" fill="#000000"/>
+ <path d="M420.265,6.428 C418.49,6.428 417.051,4.989 417.051,3.214 C417.051,1.439 418.49,0 420.265,0 C422.04,0 423.479,1.439 423.479,3.214 C423.479,4.989 422.04,6.428 420.265,6.428 z" fill="#000000"/>
+ </g>
+ <g id="dof_numers">
+ <text transform="matrix(1, 0, 0, 1, 74.305, 29.861)">
+ <tspan x="-10.08" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">(A)</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 217.3, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">0</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">1</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 311.513, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">2,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">3</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 406.305, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">4</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">5</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 217.3, 112.244)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">6</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,7</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 311.513, 112.244)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">8,9</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 397.035, 112.244)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">10</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,11</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 208.404, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">12</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,13</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 302.617, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">14,15</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 397.035, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">16</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,17</tspan>
+ </text>
+ </g>
+ <g id="node_numbers"/>
+ <g id="element_numbers">
+ <g>
+ <path d="M20.344,49.265 L129.266,49.265 L129.266,158.186 L20.344,158.186 L20.344,49.265 z" fill="#E0E0E0"/>
+ <path d="M20.344,49.265 L129.266,49.265 L129.266,158.186 L20.344,158.186 L20.344,49.265 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ </g>
+ <g>
+ <path d="M129.266,49.265 L144.204,49.265" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M144.204,52.265 L152.204,49.265 L144.204,46.265 M144.204,49.265 L152.204,49.265" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-opacity="1"/>
+ </g>
+ <g>
+ <path d="M129.266,158.186 L144.204,158.186" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M144.204,161.186 L152.204,158.186 L144.204,155.186 M144.204,158.186 L152.204,158.186" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-opacity="1"/>
+ </g>
+ <g>
+ <path d="M129.266,103.726 L144.204,103.726" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M144.204,106.726 L152.204,103.726 L144.204,100.726 M144.204,103.726 L152.204,103.726" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-opacity="1"/>
+ </g>
+ <path d="M12.947,171.034 L20.388,158.186 L27.829,171.034 L12.947,171.034 z" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M16.667,176.724 C15.096,176.724 13.822,175.451 13.822,173.879 C13.822,172.308 15.096,171.034 16.667,171.034 C18.239,171.034 19.513,172.308 19.513,173.879 C19.513,175.451 18.239,176.724 16.667,176.724 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M24.108,176.724 C22.537,176.724 21.263,175.451 21.263,173.879 C21.263,172.308 22.537,171.034 24.108,171.034 C25.68,171.034 26.954,172.308 26.954,173.879 C26.954,175.451 25.68,176.724 24.108,176.724 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M12.947,177.199 L27.829,177.199" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <path d="M12.907,170.988 L12.907,170.988" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M27.8,170.988 L27.8,170.988" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M121.835,171.034 L129.276,158.186 L136.717,171.034 L121.835,171.034 z" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M125.556,176.724 C123.984,176.724 122.71,175.451 122.71,173.879 C122.71,172.308 123.984,171.034 125.556,171.034 C127.127,171.034 128.401,172.308 128.401,173.879 C128.401,175.451 127.127,176.724 125.556,176.724 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M132.997,176.724 C131.425,176.724 130.151,175.451 130.151,173.879 C130.151,172.308 131.425,171.034 132.997,171.034 C134.568,171.034 135.842,172.308 135.842,173.879 C135.842,175.451 134.568,176.724 132.997,176.724 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M121.835,177.199 L136.717,177.199" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <path d="M121.795,170.988 L121.795,170.988" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M136.689,170.988 L136.689,170.988" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M6.665,41.861 L19.513,49.302 L6.665,56.743 L6.665,41.861 z" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M0.974,45.581 C0.974,44.01 2.248,42.736 3.82,42.736 C5.391,42.736 6.665,44.01 6.665,45.581 C6.665,47.152 5.391,48.426 3.82,48.426 C2.248,48.426 0.974,47.152 0.974,45.581 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M0.974,53.022 C0.974,51.451 2.248,50.177 3.82,50.177 C5.391,50.177 6.665,51.451 6.665,53.022 C6.665,54.593 5.391,55.867 3.82,55.867 C2.248,55.867 0.974,54.593 0.974,53.022 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M0.5,41.861 L0.5,56.743" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <path d="M6.711,41.82 L6.711,41.82" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M6.711,56.714 L6.711,56.714" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M6.665,150.749 L19.513,158.19 L6.665,165.631 L6.665,150.749 z" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M0.974,154.469 C0.974,152.898 2.248,151.624 3.82,151.624 C5.391,151.624 6.665,152.898 6.665,154.469 C6.665,156.041 5.391,157.315 3.82,157.315 C2.248,157.315 0.974,156.041 0.974,154.469 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M0.974,161.91 C0.974,160.339 2.248,159.065 3.82,159.065 C5.391,159.065 6.665,160.339 6.665,161.91 C6.665,163.482 5.391,164.756 3.82,164.756 C2.248,164.756 0.974,163.482 0.974,161.91 z" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M0.5,150.749 L0.5,165.631" fill-opacity="0" stroke="#000000" stroke-width="0.5"/>
+ <path d="M6.711,150.709 L6.711,150.709" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ <path d="M6.711,165.602 L6.711,165.602" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-linecap="round" stroke-linejoin="round"/>
+ </g>
+</svg>
diff --git a/docs/figures/Vector/mesh_dofs-b.idraw b/docs/figures/Vector/mesh_dofs-b.idraw
new file mode 100644
index 0000000..047bcee
Binary files /dev/null and b/docs/figures/Vector/mesh_dofs-b.idraw differ
diff --git a/docs/figures/Vector/mesh_dofs-b.svg b/docs/figures/Vector/mesh_dofs-b.svg
new file mode 100644
index 0000000..57aa403
--- /dev/null
+++ b/docs/figures/Vector/mesh_dofs-b.svg
@@ -0,0 +1,92 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
+<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0" y="0" width="423.479" height="218.059" viewBox="0, 0, 423.479, 218.059">
+ <defs>
+ <pattern id="Pattern_1" patternUnits="userSpaceOnUse" x="21.844" y="157.686" width="48" height="48" patternTransform="matrix(0.177, 0.177, -0.177, 0.177, 0, 0)">
+ <g transform="translate(0, -0)">
+ <path d="M12,0 L16,0 L16,48 L12,48 L12,0 z" fill="#000000"/>
+ <path d="M28,0 L32,0 L32,48 L28,48 L28,0 z" fill="#000000"/>
+ <path d="M44,0 L48,0 L48,48 L44,48 L44,0 z" fill="#000000"/>
+ </g>
+ </pattern>
+ </defs>
+ <g id="mesh">
+ <path d="M232.133,192.346 L420.765,192.346" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,98.03 L420.765,98.03" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,3.714 L420.765,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M420.765,192.346 L420.765,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M326.449,192.346 L326.449,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M232.133,192.346 L232.133,3.714" fill-opacity="0" stroke="#000000" stroke-width="1"/>
+ <path d="M231.634,195.059 C229.859,195.059 228.42,193.621 228.42,191.846 C228.42,190.071 229.859,188.632 231.634,188.632 C233.408,188.632 234.847,190.071 234.847,191.846 C234.847,193.621 233.408,195.059 231.634,195.059 z" fill="#000000"/>
+ <path d="M231.634,100.744 C229.859,100.744 228.42,99.305 228.42,97.53 C228.42,95.755 229.859,94.316 231.634,94.316 C233.408,94.316 234.847,95.755 234.847,97.53 C234.847,99.305 233.408,100.744 231.634,100.744 z" fill="#000000"/>
+ <path d="M231.634,6.428 C229.859,6.428 228.42,4.989 228.42,3.214 C228.42,1.439 229.859,0 231.634,0 C233.408,0 234.847,1.439 234.847,3.214 C234.847,4.989 233.408,6.428 231.634,6.428 z" fill="#000000"/>
+ <path d="M325.949,195.059 C324.174,195.059 322.736,193.621 322.736,191.846 C322.736,190.071 324.174,188.632 325.949,188.632 C327.724,188.632 329.163,190.071 329.163,191.846 C329.163,193.621 327.724,195.059 325.949,195.059 z" fill="#000000"/>
+ <path d="M325.949,100.744 C324.174,100.744 322.736,99.305 322.736,97.53 C322.736,95.755 324.174,94.316 325.949,94.316 C327.724,94.316 329.163,95.755 329.163,97.53 C329.163,99.305 327.724,100.744 325.949,100.744 z" fill="#000000"/>
+ <path d="M325.949,6.428 C324.174,6.428 322.736,4.989 322.736,3.214 C322.736,1.439 324.174,0 325.949,0 C327.724,0 329.163,1.439 329.163,3.214 C329.163,4.989 327.724,6.428 325.949,6.428 z" fill="#000000"/>
+ <path d="M420.265,195.059 C418.49,195.059 417.051,193.621 417.051,191.846 C417.051,190.071 418.49,188.632 420.265,188.632 C422.04,188.632 423.479,190.071 423.479,191.846 C423.479,193.621 422.04,195.059 420.265,195.059 z" fill="#000000"/>
+ <path d="M420.265,100.744 C418.49,100.744 417.051,99.305 417.051,97.53 C417.051,95.755 418.49,94.316 420.265,94.316 C422.04,94.316 423.479,95.755 423.479,97.53 C423.479,99.305 422.04,100.744 420.265,100.744 z" fill="#000000"/>
+ <path d="M420.265,6.428 C418.49,6.428 417.051,4.989 417.051,3.214 C417.051,1.439 418.49,0 420.265,0 C422.04,0 423.479,1.439 423.479,3.214 C423.479,4.989 422.04,6.428 420.265,6.428 z" fill="#000000"/>
+ </g>
+ <g id="dof_numers">
+ <text transform="matrix(1, 0, 0, 1, 217.3, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">0</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">1</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 311.513, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#690808">2</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">3</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 406.305, 206.559)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">4</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">5</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 217.3, 112.244)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">6,7</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 311.513, 112.244)">
+ <tspan x="-11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">8,9</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 397.035, 112.244)">
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#4118B8">6</tspan>
+ <tspan x="6.672" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="11.12" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#4118B8">7</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 208.404, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">10</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#690808">11</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 302.617, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#690808">12</tspan>
+ <tspan x="-2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">,</tspan>
+ <tspan x="2.224" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#690808">13</tspan>
+ </text>
+ <text transform="matrix(1, 0, 0, 1, 397.035, 17.928)">
+ <tspan x="-20.016" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#6A0808">14,15</tspan>
+ </text>
+ </g>
+ <g id="example">
+ <path d="M21.844,48.765 L130.766,48.765 L130.766,157.686 L21.844,157.686 L21.844,48.765 z" fill="#E0E0E0"/>
+ <text transform="matrix(1, 0, 0, 1, 76.305, 29.861)">
+ <tspan x="-10.08" y="4.5" font-family="Avenir-Oblique" font-size="16" fill="#000000">(B)</tspan>
+ </text>
+ <g>
+ <path d="M21.844,157.686 L130.766,157.686 L130.766,167.848 L21.844,167.848 L21.844,157.686 z" fill="url(#Pattern_1)"/>
+ <path d="M21.844,157.686 L130.766,157.686 L130.766,167.848 L21.844,167.848 L21.844,157.686 z" fill-opacity="0" stroke="#000000" stroke-width="0.96"/>
+ </g>
+ <path d="M22.344,39.103 L131.266,39.103 L131.266,49.265 L22.344,49.265 L22.344,39.103 z" fill-opacity="0" stroke="#000000" stroke-width="0.96"/>
+ <g>
+ <path d="M57.182,44.184 L87.788,44.184" fill-opacity="0" stroke="#000000" stroke-width="0.96"/>
+ <path d="M87.788,47.064 L95.468,44.184 L87.788,41.304 M87.788,44.184 L95.468,44.184" fill-opacity="0" stroke="#000000" stroke-width="0.96" stroke-opacity="1"/>
+ </g>
+ <path d="M131.266,49.265 L131.266,158.186" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ <path d="M131.266,49.265 L152.212,49.265" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ <path d="M131.266,158.186 L152.212,158.186" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ <path d="M22.344,49.265 L22.344,158.186" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ <path d="M22.344,49.265 L1.399,49.265" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ <path d="M22.344,158.186 L1.399,158.186" fill-opacity="0" stroke="#000000" stroke-width="1" stroke-dasharray="6,5"/>
+ </g>
+</svg>
diff --git a/docs/theory/examples/cartesian_linear-elasticity.py b/docs/theory/examples/cartesian_linear-elasticity.py
index f0e11e1..bb3532b 100644
--- a/docs/theory/examples/cartesian_linear-elasticity.py
+++ b/docs/theory/examples/cartesian_linear-elasticity.py
@@ -1,370 +1,370 @@
# ==================================================================================================
#
# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
#
# ==================================================================================================
import numpy as np
np.set_printoptions(linewidth=200)
# ==================================================================================================
# tensor products
# ==================================================================================================
ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
nx = 20
ny = 20
# mesh dimensions
nelem = nx * ny # number of elements
nnode = (nx+1) * (ny+1) # number of nodes
nne = 4 # number of nodes per element
ndim = 2 # number of dimensions
ndof = nnode * ndim # total number of degrees of freedom
# out-of-plane thickness
thick = 1.
# zero-initialise coordinates, displacements, and connectivity
coor = np.zeros((nnode,ndim), dtype='float')
disp = np.zeros((nnode,ndim), dtype='float')
conn = np.zeros((nelem,nne ), dtype='int' )
# coordinates
# - grid of points
x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
# - store from grid of points
coor[:,0] = x.ravel()
coor[:,1] = y.ravel()
# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(ny+1, nx+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
# - node sets
nodesLeft = inode[ :, 0]
nodesRight = inode[ :,-1]
nodesBottom = inode[ 0, :]
nodesTop = inode[-1, :]
# DOF-numbers per node
dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
# number of integration points
nip = 4
# ==================================================================================================
# gradient of the shape functions at each integration point
# ==================================================================================================
# allocate
dNx = np.empty((nelem,nip,nne,ndim))
vol = np.empty((nelem,nip))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal coordinates
xe = coor[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
# Jacobian
Je = np.einsum('mi,mj->ij', dNdxi, xe)
invJe = np.linalg.inv(Je)
detJe = np.linalg.det(Je)
# shape function gradients (w.r.t. the global coordinates)
dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
# store for later use
dNx[e,q,:,:] = dNdxe
vol[e,q] = w * detJe * thick
# ==================================================================================================
# stiffness tensor at each integration point (provides constitutive response and 'tangent')
# ==================================================================================================
# identity tensors (per integration point)
i = np.eye(3)
I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
I4s = (I4+I4rt)/2.
II = dyad22(I,I)
I4d = I4s-II/3.
# bulk and shear modulus (homogeneous)
K = 0.8333333333333333 * np.ones((nelem, nip))
G = 0.3846153846153846 * np.ones((nelem, nip))
# elasticity tensor (per integration point)
C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
# ==================================================================================================
# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
# zero-initialise strain tensor per integration point
# (plain strain -> all z-components are not written and should remain zero at all times)
Eps = np.zeros((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal displacements
ue = disp[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
dNdxe = dNx[e,q,:,:]
# displacement gradient
gradu = np.einsum('mi,mj->ij', dNdxe, ue)
# compute strain tensor, and store per integration point
# (use plain strain to convert 2-d to 3-d tensor)
Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
# constitutive response: strain tensor -> stress tensor (per integration point)
Sig = ddot42(C4, Eps)
# ==================================================================================================
# internal force from stress
# ==================================================================================================
# allocate internal force
fint = np.zeros((ndof))
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element internal force
fe = np.zeros((nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
# (plane strain: stress in z-direction irrelevant)
sig = Sig[e,q,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element internal force
# fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
- # assemble to global stiffness matrix
+ # assemble to global internal force column
iie = dofs[nodes,:].ravel()
fint[iie] += fe
# ==================================================================================================
# stiffness matrix from 'tangent'
# ==================================================================================================
# allocate stiffness matrix
K = np.zeros((ndof, ndof))
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element stiffness matrix
Ke = np.zeros((nne*ndim, nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
c4 = C4 [e,q,:2,:2,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element stiffness matrix
# Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * dV
Ke += ( np.einsum('mi,ijkl,nl->mjnk', dNdxe, c4, dNdxe) * dV ).reshape(nne*ndim, nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
# prescribed external force: zero on all free DOFs
# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
fext = np.zeros((ndof))
# DOF-partitioning: ['u'nknown, 'p'rescribed]
# - prescribed
iip = np.hstack((
dofs[nodesBottom,1],
dofs[nodesLeft ,0],
dofs[nodesRight ,0],
))
# - unknown
iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
up = np.hstack((
0.0 * np.ones((len(nodesBottom))),
0.0 * np.ones((len(nodesLeft ))),
0.1 * np.ones((len(nodesRight ))),
))
# residual force
r = fext - fint
# partition
# - stiffness matrix
Kuu = K[np.ix_(iiu, iiu)]
Kup = K[np.ix_(iiu, iip)]
Kpu = K[np.ix_(iip, iiu)]
Kpp = K[np.ix_(iip, iip)]
# - residual force
ru = r[iiu]
# solve for unknown displacement DOFs
uu = np.linalg.solve(Kuu, ru - Kup.dot(up))
# assemble
u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
# convert to nodal displacements
disp = u[dofs]
# ==================================================================================================
# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
# zero-initialise strain tensor per integration point
# (plain strain -> all z-components are not written and should remain zero at all times)
Eps = np.zeros((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal displacements
ue = disp[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
dNdxe = dNx[e,q,:,:]
# displacement gradient
gradu = np.einsum('mi,mj->ij', dNdxe, ue)
# compute strain tensor, and store per integration point
# (use plain strain to convert 2-d to 3-d tensor)
Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
# constitutive response: strain tensor -> stress tensor (per integration point)
Sig = ddot42(C4, Eps)
# ==================================================================================================
# internal force from stress, reaction (external) force from internal force on fixed nodes
# ==================================================================================================
# allocate internal force
fint = np.zeros((ndof))
# loop over nodes
for e, nodes in enumerate(conn):
# allocate element internal force
fe = np.zeros((nne*ndim))
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
# (plane strain: stress in z-direction irrelevant)
sig = Sig[e,q,:2,:2]
dNdxe = dNx[e,q,:,:]
dV = vol[e,q]
# assemble to element internal force
# fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
fint[iie] += fe
# reaction force
fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=2, figsize=(16,8))
# reconstruct external force as nodal vectors
fext = fext[dofs]
# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()
diff --git a/include/GooseFEM/ElementQuad4.h b/include/GooseFEM/ElementQuad4.h
index 2342622..76f20b3 100644
--- a/include/GooseFEM/ElementQuad4.h
+++ b/include/GooseFEM/ElementQuad4.h
@@ -1,170 +1,174 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_ELEMENTQUAD4_H
#define GOOSEFEM_ELEMENTQUAD4_H
// -------------------------------------------------------------------------------------------------
#include "config.h"
// =================================================================================================
namespace GooseFEM {
namespace Element {
namespace Quad4 {
// -------------------------------------------------------------------------------------------------
inline double inv(const xt::xtensor_fixed<double, xt::xshape<2,2>> &A,
xt::xtensor_fixed<double, xt::xshape<2,2>> &Ainv);
// -------------------------------------------------------------------------------------------------
namespace Gauss {
inline size_t nip(); // number of integration points
inline xt::xtensor<double,2> xi(); // integration point coordinates (local coordinates)
inline xt::xtensor<double,1> w(); // integration point weights
}
// -------------------------------------------------------------------------------------------------
namespace Nodal {
inline size_t nip(); // number of integration points
inline xt::xtensor<double,2> xi(); // integration point coordinates (local coordinates)
inline xt::xtensor<double,1> w(); // integration point weights
}
// -------------------------------------------------------------------------------------------------
namespace MidPoint {
inline size_t nip(); // number of integration points
inline xt::xtensor<double,2> xi(); // integration point coordinates (local coordinates)
inline xt::xtensor<double,1> w(); // integration point weights
}
// -------------------------------------------------------------------------------------------------
class Quadrature
{
public:
// fixed dimensions:
// ndim = 2 - number of dimensions
// nne = 4 - number of nodes per element
//
// naming convention:
// "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim]
// "elemvec" - nodal vectors stored per element - [nelem, nne, ndim]
// "qtensor" - integration point tensor - [nelem, nip, ndim, ndim]
// "qscalar" - integration point scalar - [nelem, nip]
// constructor: integration point coordinates and weights are optional (default: Gauss)
Quadrature() = default;
Quadrature(const xt::xtensor<double,3> &x);
Quadrature(const xt::xtensor<double,3> &x, const xt::xtensor<double,2> &xi, const xt::xtensor<double,1> &w);
// update the nodal positions (shape of "x" should match the earlier definition)
void update_x(const xt::xtensor<double,3> &x);
// return dimensions
size_t nelem() const; // number of elements
size_t nne() const; // number of nodes per element
size_t ndim() const; // number of dimension
size_t nip() const; // number of integration points
+ // return shape function gradients
+
+ xt::xtensor<double,4> gradN() const;
+
// return integration volume
void dV(xt::xtensor<double,2> &qscalar) const;
void dV(xt::xtensor<double,4> &qtensor) const; // same volume for all tensor components
// dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part
void gradN_vector(const xt::xtensor<double,3> &elemvec,
xt::xtensor<double,4> &qtensor) const;
void gradN_vector_T(const xt::xtensor<double,3> &elemvec,
xt::xtensor<double,4> &qtensor) const;
void symGradN_vector(const xt::xtensor<double,3> &elemvec,
xt::xtensor<double,4> &qtensor) const;
// integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV"
void int_N_scalar_NT_dV(const xt::xtensor<double,2> &qscalar,
xt::xtensor<double,3> &elemmat) const;
// integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV"
void int_gradN_dot_tensor2_dV(const xt::xtensor<double,4> &qtensor,
xt::xtensor<double,3> &elemvec) const;
// integral of the dot product "elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV"
void int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor<double,6> &qtensor,
xt::xtensor<double,3> &elemmat) const;
// auto allocation of the functions above
xt::xtensor<double,2> dV() const;
xt::xtensor<double,4> dVtensor() const;
xt::xtensor<double,4> gradN_vector(const xt::xtensor<double,3> &elemvec) const;
xt::xtensor<double,4> gradN_vector_T(const xt::xtensor<double,3> &elemvec) const;
xt::xtensor<double,4> symGradN_vector(const xt::xtensor<double,3> &elemvec) const;
xt::xtensor<double,3> int_N_scalar_NT_dV(const xt::xtensor<double,2> &qscalar) const;
xt::xtensor<double,3> int_gradN_dot_tensor2_dV(const xt::xtensor<double,4> &qtensor) const;
xt::xtensor<double,3> int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor<double,6> &qtensor) const;
private:
// compute "vol" and "dNdx" based on current "x"
void compute_dN();
private:
// dimensions (flexible)
size_t m_nelem; // number of elements
size_t m_nip; // number of integration points
// dimensions (fixed for this element type)
static const size_t m_nne=4; // number of nodes per element
static const size_t m_ndim=2; // number of dimensions
// data arrays
xt::xtensor<double,3> m_x; // nodal positions stored per element [nelem, nne, ndim]
xt::xtensor<double,1> m_w; // weight of each integration point [nip]
xt::xtensor<double,2> m_xi; // local coordinate of each integration point [nip, ndim]
xt::xtensor<double,2> m_N; // shape functions [nip, nne]
xt::xtensor<double,3> m_dNxi; // shape function gradients w.r.t. local coordinate [nip, nne, ndim]
xt::xtensor<double,4> m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim]
xt::xtensor<double,2> m_vol; // integration point volume [nelem, nip]
};
// -------------------------------------------------------------------------------------------------
}}} // namespace ...
// =================================================================================================
#include "ElementQuad4.hpp"
// =================================================================================================
#endif
diff --git a/include/GooseFEM/ElementQuad4.hpp b/include/GooseFEM/ElementQuad4.hpp
index 1078238..81e3aa7 100644
--- a/include/GooseFEM/ElementQuad4.hpp
+++ b/include/GooseFEM/ElementQuad4.hpp
@@ -1,674 +1,676 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_ELEMENTQUAD4_HPP
#define GOOSEFEM_ELEMENTQUAD4_HPP
// -------------------------------------------------------------------------------------------------
#include "ElementQuad4.h"
// =================================================================================================
namespace GooseFEM {
namespace Element {
namespace Quad4 {
// =================================================================================================
inline double inv(const xt::xtensor_fixed<double, xt::xshape<2,2>> &A,
xt::xtensor_fixed<double, xt::xshape<2,2>> &Ainv)
{
// compute determinant
double det = A(0,0) * A(1,1) - A(0,1) * A(1,0);
// compute inverse
Ainv(0,0) = A(1,1) / det;
Ainv(0,1) = -1. * A(0,1) / det;
Ainv(1,0) = -1. * A(1,0) / det;
Ainv(1,1) = A(0,0) / det;
return det;
}
// =================================================================================================
namespace Gauss {
// -------------------------------------------------------------------------------------------------
inline size_t nip()
{
return 4;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,2> xi()
{
size_t nip = 4;
size_t ndim = 2;
xt::xtensor<double,2> xi = xt::empty<double>({nip,ndim});
xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.);
xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.);
xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.);
xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.);
return xi;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,1> w()
{
size_t nip = 4;
xt::xtensor<double,1> w = xt::empty<double>({nip});
w(0) = 1.;
w(1) = 1.;
w(2) = 1.;
w(3) = 1.;
return w;
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
namespace Nodal {
// -------------------------------------------------------------------------------------------------
inline size_t nip()
{
return 4;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,2> xi()
{
size_t nip = 4;
size_t ndim = 2;
xt::xtensor<double,2> xi = xt::empty<double>({nip,ndim});
xi(0,0) = -1.; xi(0,1) = -1.;
xi(1,0) = +1.; xi(1,1) = -1.;
xi(2,0) = +1.; xi(2,1) = +1.;
xi(3,0) = -1.; xi(3,1) = +1.;
return xi;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,1> w()
{
size_t nip = 4;
xt::xtensor<double,1> w = xt::empty<double>({nip});
w(0) = 1.;
w(1) = 1.;
w(2) = 1.;
w(3) = 1.;
return w;
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
namespace MidPoint {
// -------------------------------------------------------------------------------------------------
inline size_t nip()
{
return 1;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,2> xi()
{
size_t nip = 1;
size_t ndim = 2;
xt::xtensor<double,2> xi = xt::empty<double>({nip,ndim});
xi(0,0) = 0.; xi(0,1) = 0.;
return xi;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,1> w()
{
size_t nip = 1;
xt::xtensor<double,1> w = xt::empty<double>({nip});
w(0) = 1.;
return w;
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
inline Quadrature::Quadrature(const xt::xtensor<double,3> &x) :
Quadrature(x, Gauss::xi(), Gauss::w()) {};
// -------------------------------------------------------------------------------------------------
inline Quadrature::Quadrature(const xt::xtensor<double,3> &x, const xt::xtensor<double,2> &xi,
const xt::xtensor<double,1> &w) : m_x(x), m_w(w), m_xi(xi)
{
// check input
assert( m_x.shape()[1] == m_nne );
assert( m_x.shape()[2] == m_ndim );
// extract number of elements and number of integration points
m_nelem = m_x.shape()[0];
m_nip = m_w.size();
// check input
assert( m_xi.shape()[0] == m_nip );
assert( m_xi.shape()[1] == m_ndim );
assert( m_w .size() == m_nip );
// allocate arrays
m_N = xt::empty<double>({ m_nip, m_nne });
m_dNxi = xt::empty<double>({ m_nip, m_nne, m_ndim});
m_dNx = xt::empty<double>({m_nelem, m_nip, m_nne, m_ndim});
m_vol = xt::empty<double>({m_nelem, m_nip });
// shape functions
for ( size_t q = 0 ; q < m_nip ; ++q )
{
m_N(q,0) = .25 * (1.-m_xi(q,0)) * (1.-m_xi(q,1));
m_N(q,1) = .25 * (1.+m_xi(q,0)) * (1.-m_xi(q,1));
m_N(q,2) = .25 * (1.+m_xi(q,0)) * (1.+m_xi(q,1));
m_N(q,3) = .25 * (1.-m_xi(q,0)) * (1.+m_xi(q,1));
}
// shape function gradients in local coordinates
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - dN / dxi_0
m_dNxi(q,0,0) = -.25*(1.-m_xi(q,1));
m_dNxi(q,1,0) = +.25*(1.-m_xi(q,1));
m_dNxi(q,2,0) = +.25*(1.+m_xi(q,1));
m_dNxi(q,3,0) = -.25*(1.+m_xi(q,1));
// - dN / dxi_1
m_dNxi(q,0,1) = -.25*(1.-m_xi(q,0));
m_dNxi(q,1,1) = -.25*(1.+m_xi(q,0));
m_dNxi(q,2,1) = +.25*(1.+m_xi(q,0));
m_dNxi(q,3,1) = +.25*(1.-m_xi(q,0));
}
// compute the shape function gradients, based on "x"
compute_dN();
}
// -------------------------------------------------------------------------------------------------
inline size_t Quadrature::nelem() const { return m_nelem; };
inline size_t Quadrature::nne() const { return m_nne; };
inline size_t Quadrature::ndim() const { return m_ndim; };
inline size_t Quadrature::nip() const { return m_nip; };
+inline xt::xtensor<double,4> Quadrature::gradN() const { return m_dNx; };
+
// -------------------------------------------------------------------------------------------------
inline void Quadrature::dV(xt::xtensor<double,2> &qscalar) const
{
assert( qscalar.shape()[0] == m_nelem );
assert( qscalar.shape()[1] == m_nip );
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
for ( size_t q = 0 ; q < m_nip ; ++q )
qscalar(e,q) = m_vol(e,q);
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::dV(xt::xtensor<double,4> &qtensor) const
{
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
for ( size_t q = 0 ; q < m_nip ; ++q )
for ( size_t i = 0 ; i < m_ndim ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
qtensor(e,q,i,j) = m_vol(e,q);
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::update_x(const xt::xtensor<double,3> &x)
{
assert( x.shape()[0] == m_nelem );
assert( x.shape()[1] == m_nne );
assert( x.shape()[2] == m_ndim );
assert( x.size() == m_x.size() );
// update positions
xt::noalias(m_x) = x;
// update the shape function gradients for the new "x"
compute_dN();
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::compute_dN()
{
#pragma omp parallel
{
// allocate local variables
xt::xtensor_fixed<double, xt::xshape<2,2>> J, Jinv;
// loop over all elements (in parallel)
#pragma omp for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias nodal positions
auto x = xt::adapt(&m_x(e,0,0), xt::xshape<m_nne,m_ndim>());
// loop over integration points
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNxi = xt::adapt(&m_dNxi( q,0,0), xt::xshape<m_nne,m_ndim>());
auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape<m_nne,m_ndim>());
// - Jacobian (loops unrolled for efficiency)
// J(i,j) += dNxi(m,i) * x(m,j);
J(0,0) = dNxi(0,0)*x(0,0) + dNxi(1,0)*x(1,0) + dNxi(2,0)*x(2,0) + dNxi(3,0)*x(3,0);
J(0,1) = dNxi(0,0)*x(0,1) + dNxi(1,0)*x(1,1) + dNxi(2,0)*x(2,1) + dNxi(3,0)*x(3,1);
J(1,0) = dNxi(0,1)*x(0,0) + dNxi(1,1)*x(1,0) + dNxi(2,1)*x(2,0) + dNxi(3,1)*x(3,0);
J(1,1) = dNxi(0,1)*x(0,1) + dNxi(1,1)*x(1,1) + dNxi(2,1)*x(2,1) + dNxi(3,1)*x(3,1);
// - determinant and inverse of the Jacobian
double Jdet = inv(J, Jinv);
// - shape function gradients wrt global coordinates (loops partly unrolled for efficiency)
// dNx(m,i) += Jinv(i,j) * dNxi(m,j);
for ( size_t m = 0 ; m < m_nne ; ++m )
{
dNx(m,0) = Jinv(0,0) * dNxi(m,0) + Jinv(0,1) * dNxi(m,1);
dNx(m,1) = Jinv(1,0) * dNxi(m,0) + Jinv(1,1) * dNxi(m,1);
}
// - integration point volume
m_vol(e,q) = m_w(q) * Jdet;
}
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::gradN_vector(
const xt::xtensor<double,3> &elemvec, xt::xtensor<double,4> &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape<m_nne,m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape<m_nne ,m_ndim>());
auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape<m_ndim,m_ndim>());
// - evaluate dyadic product (loops unrolled for efficiency)
// gradu(i,j) += dNx(m,i) * u(m,j)
gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
gradu(0,1) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1);
gradu(1,0) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0);
gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::gradN_vector_T(
const xt::xtensor<double,3> &elemvec, xt::xtensor<double,4> &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape<m_nne,m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape<m_nne ,m_ndim>());
auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape<m_ndim,m_ndim>());
// - evaluate transpose of dyadic product (loops unrolled for efficiency)
// gradu(j,i) += dNx(m,i) * u(m,j)
gradu(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
gradu(1,0) = dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1);
gradu(0,1) = dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0);
gradu(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::symGradN_vector(
const xt::xtensor<double,3> &elemvec, xt::xtensor<double,4> &qtensor) const
{
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nne );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias element vector (e.g. nodal displacements)
auto u = xt::adapt(&elemvec(e,0,0), xt::xshape<m_nne,m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape<m_nne ,m_ndim>());
auto eps = xt::adapt(&qtensor(e,q,0,0), xt::xshape<m_ndim,m_ndim>());
// - evaluate symmetrized dyadic product (loops unrolled for efficiency)
// gradu(i,j) += dNx(m,i) * u(m,j)
// eps (j,i) = 0.5 * ( gradu(i,j) + gradu(j,i) )
eps(0,0) = dNx(0,0)*u(0,0) + dNx(1,0)*u(1,0) + dNx(2,0)*u(2,0) + dNx(3,0)*u(3,0);
eps(1,1) = dNx(0,1)*u(0,1) + dNx(1,1)*u(1,1) + dNx(2,1)*u(2,1) + dNx(3,1)*u(3,1);
eps(0,1) = ( dNx(0,0)*u(0,1) + dNx(1,0)*u(1,1) + dNx(2,0)*u(2,1) + dNx(3,0)*u(3,1) +
dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0) ) * 0.5;
eps(1,0) = eps(0,1);
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::int_N_scalar_NT_dV(
const xt::xtensor<double,2> &qscalar, xt::xtensor<double,3> &elemmat) const
{
assert( qscalar.shape()[0] == m_nelem );
assert( qscalar.shape()[1] == m_nip );
assert( elemmat.shape()[0] == m_nelem );
assert( elemmat.shape()[1] == m_nne*m_ndim );
assert( elemmat.shape()[2] == m_nne*m_ndim );
// zero-initialize: matrix of matrices
elemmat.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias (e.g. mass matrix)
auto M = xt::adapt(&elemmat(e,0,0), xt::xshape<m_nne*m_ndim,m_nne*m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto N = xt::adapt(&m_N(q,0), xt::xshape<m_nne>());
auto& vol = m_vol (e,q);
auto& rho = qscalar(e,q);
// - evaluate scalar product, for all dimensions, and assemble
// M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV
for ( size_t m = 0 ; m < m_nne ; ++m ) {
for ( size_t n = 0 ; n < m_nne ; ++n ) {
M(m*m_ndim+0, n*m_ndim+0) += N(m) * rho * N(n) * vol;
M(m*m_ndim+1, n*m_ndim+1) += N(m) * rho * N(n) * vol;
}
}
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor<double,4> &qtensor,
xt::xtensor<double,3> &elemvec) const
{
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nip );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
assert( elemvec.shape()[0] == m_nelem );
assert( elemvec.shape()[1] == m_nne );
assert( elemvec.shape()[2] == m_ndim );
// zero-initialize output: matrix of vectors
elemvec.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias (e.g. nodal force)
auto f = xt::adapt(&elemvec(e,0,0), xt::xshape<m_nne,m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape<m_nne ,m_ndim>());
auto sig = xt::adapt(&qtensor(e,q,0,0), xt::xshape<m_ndim,m_ndim>());
auto& vol = m_vol(e,q);
// - evaluate dot product, and assemble
for ( size_t m = 0 ; m < m_nne ; ++m )
{
f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) ) * vol;
f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) ) * vol;
}
}
}
}
// -------------------------------------------------------------------------------------------------
inline void Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor<double,6> &qtensor,
xt::xtensor<double,3> &elemmat) const
{
assert( qtensor.shape()[0] == m_nelem );
assert( qtensor.shape()[1] == m_nip );
assert( qtensor.shape()[2] == m_ndim );
assert( qtensor.shape()[3] == m_ndim );
assert( qtensor.shape()[4] == m_ndim );
assert( qtensor.shape()[5] == m_ndim );
assert( elemmat.shape()[0] == m_nelem );
assert( elemmat.shape()[1] == m_nne*m_ndim );
assert( elemmat.shape()[2] == m_nne*m_ndim );
// zero-initialize output: matrix of vector
elemmat.fill(0.0);
// loop over all elements (in parallel)
#pragma omp parallel for
for ( size_t e = 0 ; e < m_nelem ; ++e )
{
// alias (e.g. nodal force)
auto K = xt::adapt(&elemmat(e,0,0), xt::xshape<m_nne*m_ndim,m_nne*m_ndim>());
// loop over all integration points in element "e"
for ( size_t q = 0 ; q < m_nip ; ++q )
{
// - alias
auto dNx = xt::adapt(&m_dNx(e,q,0,0), xt::xshape<m_nne,m_ndim>());
auto C = xt::adapt(&qtensor(e,q,0,0,0,0), xt::xshape<m_ndim,m_ndim,m_ndim,m_ndim>());
auto& vol = m_vol(e,q);
// - evaluate dot product, and assemble
for ( size_t m = 0 ; m < m_nne ; ++m )
for ( size_t n = 0 ; n < m_nne ; ++n )
for ( size_t i = 0 ; i < m_ndim ; ++i )
for ( size_t j = 0 ; j < m_ndim ; ++j )
for ( size_t k = 0 ; k < m_ndim ; ++k )
for ( size_t l = 0 ; l < m_ndim ; ++l )
K(m*m_ndim+j, n*m_ndim+k) += dNx(m,i) * C(i,j,k,l) * dNx(n,l) * vol;
}
}
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,2> Quadrature::dV() const
{
xt::xtensor<double,2> out = xt::empty<double>({m_nelem, m_nip});
this->dV(out);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,4> Quadrature::dVtensor() const
{
xt::xtensor<double,4> out = xt::empty<double>({m_nelem, m_nip, m_ndim, m_ndim});
this->dV(out);
return out;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,4> Quadrature::gradN_vector(const xt::xtensor<double,3> &elemvec) const
{
xt::xtensor<double,4> qtensor = xt::empty<double>({m_nelem, m_nip, m_ndim, m_ndim});
this->gradN_vector(elemvec, qtensor);
return qtensor;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,4> Quadrature::gradN_vector_T(const xt::xtensor<double,3> &elemvec) const
{
xt::xtensor<double,4> qtensor = xt::empty<double>({m_nelem, m_nip, m_ndim, m_ndim});
this->gradN_vector_T(elemvec, qtensor);
return qtensor;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,4> Quadrature::symGradN_vector(const xt::xtensor<double,3> &elemvec) const
{
xt::xtensor<double,4> qtensor = xt::empty<double>({m_nelem, m_nip, m_ndim, m_ndim});
this->symGradN_vector(elemvec, qtensor);
return qtensor;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,3> Quadrature::int_N_scalar_NT_dV(
const xt::xtensor<double,2> &qscalar) const
{
xt::xtensor<double,3> elemmat = xt::empty<double>({m_nelem, m_nne*m_ndim, m_nne*m_ndim});
this->int_N_scalar_NT_dV(qscalar, elemmat);
return elemmat;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,3> Quadrature::int_gradN_dot_tensor2_dV(
const xt::xtensor<double,4> &qtensor) const
{
xt::xtensor<double,3> elemvec = xt::empty<double>({m_nelem, m_nne, m_ndim});
this->int_gradN_dot_tensor2_dV(qtensor, elemvec);
return elemvec;
}
// -------------------------------------------------------------------------------------------------
inline xt::xtensor<double,3> Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV(
const xt::xtensor<double,6> &qtensor) const
{
xt::xtensor<double,3> elemmat = xt::empty<double>({m_nelem, m_ndim*m_nne, m_ndim*m_nne});
this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat);
return elemmat;
}
// -------------------------------------------------------------------------------------------------
}}} // namespace ...
// =================================================================================================
#endif
diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h
index ee0145e..1f94a6c 100644
--- a/include/GooseFEM/MatrixPartitioned.h
+++ b/include/GooseFEM/MatrixPartitioned.h
@@ -1,134 +1,134 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MATRIXPARTITIONED_H
#define GOOSEFEM_MATRIXPARTITIONED_H
// -------------------------------------------------------------------------------------------------
#include "config.h"
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
// =================================================================================================
namespace GooseFEM {
// -------------------------------------------------------------------------------------------------
class MatrixPartitioned
{
public:
// constructors
MatrixPartitioned() = default;
MatrixPartitioned(const xt::xtensor<size_t,2> &conn, const xt::xtensor<size_t,2> &dofs,
const xt::xtensor<size_t,1> &iip);
// dimensions
size_t nelem() const; // number of elements
size_t nne() const; // number of nodes per element
size_t nnode() const; // number of nodes
size_t ndim() const; // number of dimensions
size_t ndof() const; // number of DOFs
size_t nnu() const; // number of unknown DOFs
size_t nnp() const; // number of prescribed DOFs
// DOF lists
xt::xtensor<size_t,2> dofs() const; // DOFs
xt::xtensor<size_t,1> iiu() const; // unknown DOFs
xt::xtensor<size_t,1> iip() const; // prescribed DOFs
// assemble from matrices stored per element [nelem, nne*ndim, nne*ndim]
void assemble(const xt::xtensor<double,3> &elemmat);
// solve: x = A \ b
// x_u = A_uu \ ( b_u - A_up * x_p )
// b_p = A_pu * x_u + A_pp * x_p
void solve(xt::xtensor<double,2> &b,
xt::xtensor<double,2> &x);
void solve(xt::xtensor<double,1> &b,
xt::xtensor<double,1> &x);
void solve_u(const xt::xtensor<double,1> &b_u, const xt::xtensor<double,1> &x_p,
xt::xtensor<double,1> &x_u);
// auto allocation of the functions above
xt::xtensor<double,1> solve_u(const xt::xtensor<double,1> &b_u, const xt::xtensor<double,1> &x_p);
private:
// the matrix
Eigen::SparseMatrix<double> m_data_uu;
Eigen::SparseMatrix<double> m_data_up;
Eigen::SparseMatrix<double> m_data_pu;
Eigen::SparseMatrix<double> m_data_pp;
// the matrix to assemble
std::vector<Eigen::Triplet<double>> m_trip_uu;
std::vector<Eigen::Triplet<double>> m_trip_up;
std::vector<Eigen::Triplet<double>> m_trip_pu;
std::vector<Eigen::Triplet<double>> m_trip_pp;
// solver (re-used to solve different RHS)
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> m_solver;
// signal changes to data compare to the last inverse
bool m_change=false;
// bookkeeping
xt::xtensor<size_t,2> m_conn; // connectivity [nelem, nne ]
xt::xtensor<size_t,2> m_dofs; // DOF-numbers per node [nnode, ndim]
xt::xtensor<size_t,2> m_part; // DOF-numbers per node, renumbered [nnode, ndim]
xt::xtensor<size_t,1> m_iiu; // DOF-numbers that are unknown [nnu]
xt::xtensor<size_t,1> m_iip; // DOF-numbers that are prescribed [nnp]
// dimensions
size_t m_nelem; // number of elements
size_t m_nne; // number of nodes per element
size_t m_nnode; // number of nodes
size_t m_ndim; // number of dimensions
size_t m_ndof; // number of DOFs
size_t m_nnu; // number of unknown DOFs
size_t m_nnp; // number of prescribed DOFs
// compute inverse (automatically evaluated by "solve")
void factorize();
- // convert arrays (see VectorPartitioned, which contains public functions)
+ // convert arrays (Eigen version of VectorPartitioned, which contains public functions)
Eigen::VectorXd asDofs_u(const xt::xtensor<double,1> &dofval) const;
Eigen::VectorXd asDofs_u(const xt::xtensor<double,2> &nodevec) const;
Eigen::VectorXd asDofs_p(const xt::xtensor<double,1> &dofval) const;
Eigen::VectorXd asDofs_p(const xt::xtensor<double,2> &nodevec) const;
};
// -------------------------------------------------------------------------------------------------
} // namespace ...
// =================================================================================================
#include "MatrixPartitioned.hpp"
// =================================================================================================
#endif
diff --git a/include/GooseFEM/python.cpp b/include/GooseFEM/python.cpp
index a4deb7e..257287a 100644
--- a/include/GooseFEM/python.cpp
+++ b/include/GooseFEM/python.cpp
@@ -1,698 +1,727 @@
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#include <Eigen/Eigen>
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <cppmat/pybind11.h>
#include <pyxtensor/pyxtensor.hpp>
#include "GooseFEM.h"
// =================================================================================================
namespace py = pybind11;
namespace M = GooseFEM;
// =================================================================================================
PYBIND11_MODULE(GooseFEM, m) {
m.doc() = "Some simple finite element meshes and operations";
// =================================================================================================
+py::class_<GooseFEM::Vector>(m, "Vector")
+
+ .def(py::init<const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,2> &>(), "Switch between dofval/nodevec/elemvec", py::arg("conn"), py::arg("dofs"))
+
+ .def("nelem", &M::Vector::nelem, "Return number of element")
+ .def("nne" , &M::Vector::nne , "Return number of nodes per element")
+ .def("nnode", &M::Vector::nnode, "Return number of nodes")
+ .def("ndim" , &M::Vector::ndim , "Return number of dimensions")
+ .def("ndof" , &M::Vector::ndof , "Return number of degrees-of-freedom")
+ .def("dofs" , &M::Vector::dofs , "Return degrees-of-freedom")
+
+ .def("asDofs" , py::overload_cast<const xt::xtensor<double,2>&>(&M::Vector::asDofs , py::const_), "Set 'dofval" , py::arg("nodevec"))
+ .def("asDofs" , py::overload_cast<const xt::xtensor<double,3>&>(&M::Vector::asDofs , py::const_), "Set 'dofval" , py::arg("elemvec"))
+
+ .def("asNode" , py::overload_cast<const xt::xtensor<double,1>&>(&M::Vector::asNode , py::const_), "Set 'nodevec", py::arg("dofval"))
+ .def("asNode" , py::overload_cast<const xt::xtensor<double,3>&>(&M::Vector::asNode , py::const_), "Set 'nodevec", py::arg("elemvec"))
+
+ .def("asElement", py::overload_cast<const xt::xtensor<double,1>&>(&M::Vector::asElement, py::const_), "Set 'elemvec", py::arg("dofval"))
+ .def("asElement", py::overload_cast<const xt::xtensor<double,2>&>(&M::Vector::asElement, py::const_), "Set 'elemvec", py::arg("nodevec"))
+
+ .def("assembleDofs", py::overload_cast<const xt::xtensor<double,2>&>(&M::Vector::assembleDofs, py::const_), "Assemble 'dofval'" , py::arg("nodevec"))
+ .def("assembleDofs", py::overload_cast<const xt::xtensor<double,3>&>(&M::Vector::assembleDofs, py::const_), "Assemble 'dofval'" , py::arg("elemvec"))
+
+ .def("assembleNode", py::overload_cast<const xt::xtensor<double,3>&>(&M::Vector::assembleNode, py::const_), "Assemble 'nodevec'", py::arg("elemvec"))
+
+ .def("__repr__", [](const GooseFEM::Vector &){ return "<GooseFEM.Vector>"; });
+
+// =================================================================================================
+
py::class_<GooseFEM::VectorPartitioned>(m, "VectorPartitioned")
.def(py::init<const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,1> &>(), "Switch between dofval/nodevec/elemvec", py::arg("conn"), py::arg("dofs"), py::arg("iip"))
.def("nelem", &M::VectorPartitioned::nelem, "Return number of element")
.def("nne" , &M::VectorPartitioned::nne , "Return number of nodes per element")
.def("nnode", &M::VectorPartitioned::nnode, "Return number of nodes")
.def("ndim" , &M::VectorPartitioned::ndim , "Return number of dimensions")
.def("ndof" , &M::VectorPartitioned::ndof , "Return number of degrees-of-freedom")
.def("nnu" , &M::VectorPartitioned::nnu , "Return number of unknown degrees-of-freedom")
.def("nnp" , &M::VectorPartitioned::nnp , "Return number of prescribed degrees-of-freedom")
.def("dofs" , &M::VectorPartitioned::dofs , "Return degrees-of-freedom")
.def("iiu" , &M::VectorPartitioned::iiu , "Return unknown degrees-of-freedom")
.def("iip" , &M::VectorPartitioned::iip , "Return prescribed degrees-of-freedom")
.def("asDofs" , py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(&M::VectorPartitioned::asDofs , py::const_), "Set 'dofval" , py::arg("dofval_u"), py::arg("dofval_p"))
.def("asDofs" , py::overload_cast<const xt::xtensor<double,2>& >(&M::VectorPartitioned::asDofs , py::const_), "Set 'dofval" , py::arg("nodevec"))
.def("asDofs" , py::overload_cast<const xt::xtensor<double,3>& >(&M::VectorPartitioned::asDofs , py::const_), "Set 'dofval" , py::arg("elemvec"))
.def("asDofs_u" , py::overload_cast<const xt::xtensor<double,2>& >(&M::VectorPartitioned::asDofs_u , py::const_), "Set 'dofval" , py::arg("nodevec"))
.def("asDofs_u" , py::overload_cast<const xt::xtensor<double,3>& >(&M::VectorPartitioned::asDofs_u , py::const_), "Set 'dofval" , py::arg("elemvec"))
.def("asDofs_p" , py::overload_cast<const xt::xtensor<double,2>& >(&M::VectorPartitioned::asDofs_p , py::const_), "Set 'dofval" , py::arg("nodevec"))
.def("asDofs_p" , py::overload_cast<const xt::xtensor<double,3>& >(&M::VectorPartitioned::asDofs_p , py::const_), "Set 'dofval" , py::arg("elemvec"))
.def("asNode" , py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(&M::VectorPartitioned::asNode , py::const_), "Set 'nodevec", py::arg("dofval_u"), py::arg("dofval_p"))
.def("asNode" , py::overload_cast<const xt::xtensor<double,1>& >(&M::VectorPartitioned::asNode , py::const_), "Set 'nodevec", py::arg("dofval"))
.def("asNode" , py::overload_cast<const xt::xtensor<double,3>& >(&M::VectorPartitioned::asNode , py::const_), "Set 'nodevec", py::arg("elemvec"))
.def("asElement", py::overload_cast<const xt::xtensor<double,1>&,const xt::xtensor<double,1>&>(&M::VectorPartitioned::asElement, py::const_), "Set 'elemvec", py::arg("dofval_u"), py::arg("dofval_p"))
.def("asElement", py::overload_cast<const xt::xtensor<double,1>& >(&M::VectorPartitioned::asElement, py::const_), "Set 'elemvec", py::arg("dofval"))
.def("asElement", py::overload_cast<const xt::xtensor<double,2>& >(&M::VectorPartitioned::asElement, py::const_), "Set 'elemvec", py::arg("nodevec"))
.def("assembleDofs" , py::overload_cast<const xt::xtensor<double,2>&>(&M::VectorPartitioned::assembleDofs , py::const_), "Assemble 'dofval'" , py::arg("nodevec"))
.def("assembleDofs" , py::overload_cast<const xt::xtensor<double,3>&>(&M::VectorPartitioned::assembleDofs , py::const_), "Assemble 'dofval'" , py::arg("elemvec"))
.def("assembleDofs_u", py::overload_cast<const xt::xtensor<double,2>&>(&M::VectorPartitioned::assembleDofs_u, py::const_), "Assemble 'dofval'" , py::arg("nodevec"))
.def("assembleDofs_u", py::overload_cast<const xt::xtensor<double,3>&>(&M::VectorPartitioned::assembleDofs_u, py::const_), "Assemble 'dofval'" , py::arg("elemvec"))
.def("assembleDofs_p", py::overload_cast<const xt::xtensor<double,2>&>(&M::VectorPartitioned::assembleDofs_p, py::const_), "Assemble 'dofval'" , py::arg("nodevec"))
.def("assembleDofs_p", py::overload_cast<const xt::xtensor<double,3>&>(&M::VectorPartitioned::assembleDofs_p, py::const_), "Assemble 'dofval'" , py::arg("elemvec"))
.def("assembleNode" , py::overload_cast<const xt::xtensor<double,3>&>(&M::VectorPartitioned::assembleNode , py::const_), "Assemble 'nodevec'", py::arg("elemvec"))
- .def("__repr__", [](const GooseFEM::VectorPartitioned &){ return "<GooseFEM.Vector>"; });
+ .def("__repr__", [](const GooseFEM::VectorPartitioned &){ return "<GooseFEM.VectorPartitioned>"; });
// =================================================================================================
py::class_<GooseFEM::MatrixDiagonalPartitioned>(m, "MatrixDiagonalPartitioned")
.def(py::init<const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,2> &, const xt::xtensor<size_t,1> &>(), "Diagonal matrix", py::arg("conn"), py::arg("dofs"), py::arg("iip"))
.def("nelem", &M::MatrixDiagonalPartitioned::nelem, "Return number of element")
.def("nne" , &M::MatrixDiagonalPartitioned::nne , "Return number of nodes per element")
.def("nnode", &M::MatrixDiagonalPartitioned::nnode, "Return number of nodes")
.def("ndim" , &M::MatrixDiagonalPartitioned::ndim , "Return number of dimensions")
.def("ndof" , &M::MatrixDiagonalPartitioned::ndof , "Return number of degrees-of-freedom")
.def("nnu" , &M::MatrixDiagonalPartitioned::nnu , "Return number of unknown degrees-of-freedom")
.def("nnp" , &M::MatrixDiagonalPartitioned::nnp , "Return number of prescribed degrees-of-freedom")
.def("assemble", &M::MatrixDiagonalPartitioned::assemble, "Assemble matrix from 'elemmat", py::arg("elemmat"))
.def("dofs" , &M::MatrixDiagonalPartitioned::dofs , "Return degrees-of-freedom")
.def("iiu" , &M::MatrixDiagonalPartitioned::iiu , "Return unknown degrees-of-freedom")
.def("iip" , &M::MatrixDiagonalPartitioned::iip , "Return prescribed degrees-of-freedom")
.def("dot" , py::overload_cast<const xt::xtensor<double,1>& >(&M::MatrixDiagonalPartitioned::dot , py::const_), "Dot product 'b_i = A_ij * x_j" , py::arg("x"))
.def("dot_u", py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(&M::MatrixDiagonalPartitioned::dot_u, py::const_), "Dot product 'b_i = A_ij * x_j (b_u = A_uu * x_u + A_up * x_p == A_uu * x_u)", py::arg("x_u"), py::arg("x_p"))
.def("dot_p", py::overload_cast<const xt::xtensor<double,1>&, const xt::xtensor<double,1>&>(&M::MatrixDiagonalPartitioned::dot_p, py::const_), "Dot product 'b_i = A_ij * x_j (b_p = A_pu * x_u + A_pp * x_p == A_pp * x_p)", py::arg("x_u"), py::arg("x_p"))
// .def("solve_u", &M::MatrixDiagonalPartitioned::solve_u, "Solve 'x_u = A_uu \\ ( b_u - A_up * x_p ) == A_uu \\ b_u'", py::arg("b_u"), py::arg("x_p"))
.def("asDiagonal", &M::MatrixDiagonalPartitioned::asDiagonal, "Return as diagonal matrix (column)")
.def("__repr__", [](const GooseFEM::MatrixDiagonalPartitioned &){ return "<GooseFEM.MatrixDiagonalPartitioned>"; });
// =================================================================================================
py::module mElement = m.def_submodule("Element", "Generic element routines");
// -------------------------------------------------------------------------------------------------
mElement.def("asElementVector" , &GooseFEM::Element::asElementVector , "Covert to 'elemvec'", py::arg("conn"), py::arg("nodevec"));
mElement.def("assembleElementVector", &GooseFEM::Element::assembleNodeVector, "Assemble 'nodevec'" , py::arg("conn"), py::arg("elemvec"));
// =================================================================================================
{
py::module sm = mElement.def_submodule("Quad4", "Linear quadrilateral elements (2D)");
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Element::Quad4::Quadrature>(sm, "Quadrature")
.def(py::init<const xt::xtensor<double,3> &>(), "Quadrature", py::arg("x"))
.def(py::init<const xt::xtensor<double,3> &, const xt::xtensor<double,2> &, const xt::xtensor<double,1> &>(), "Quadrature", py::arg("x"), py::arg("xi"), py::arg("w"))
.def("nelem", &GooseFEM::Element::Quad4::Quadrature::nelem, "Return number of elements" )
.def("nne" , &GooseFEM::Element::Quad4::Quadrature::nne , "Return number of nodes per element" )
.def("ndim" , &GooseFEM::Element::Quad4::Quadrature::ndim , "Return number of dimensions" )
.def("nip" , &GooseFEM::Element::Quad4::Quadrature::nip , "Return number of integration points")
.def("dV" , py::overload_cast<>(&GooseFEM::Element::Quad4::Quadrature::dV , py::const_), "Integration point volume (qscalar)")
.def("dVtensor", py::overload_cast<>(&GooseFEM::Element::Quad4::Quadrature::dVtensor, py::const_), "Integration point volume (qtensor)")
.def("gradN_vector" , py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Quad4::Quadrature::gradN_vector , py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("gradN_vector_T" , py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Quad4::Quadrature::gradN_vector_T , py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("symGradN_vector", py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Quad4::Quadrature::symGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("int_N_scalar_NT_dV" , py::overload_cast<const xt::xtensor<double,2> &>(&GooseFEM::Element::Quad4::Quadrature::int_N_scalar_NT_dV , py::const_), "Integration, returns 'elemmat'", py::arg("qscalar"))
.def("int_gradN_dot_tensor2_dV", py::overload_cast<const xt::xtensor<double,4> &>(&GooseFEM::Element::Quad4::Quadrature::int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor"))
.def("__repr__", [](const GooseFEM::Element::Quad4::Quadrature &){ return "<GooseFEM.Element.Quad4.Quadrature>"; });
// -------------------------------------------------------------------------------------------------
{
py::module ssm = sm.def_submodule("Gauss", "Gauss quadrature");
ssm.def("nip", &GooseFEM::Element::Quad4::Gauss::nip, "Return number of integration point" );
ssm.def("xi" , &GooseFEM::Element::Quad4::Gauss::xi , "Return integration point coordinates");
ssm.def("w" , &GooseFEM::Element::Quad4::Gauss::w , "Return integration point weights" );
}
// -------------------------------------------------------------------------------------------------
{
py::module ssm = sm.def_submodule("Nodal", "Nodal quadrature");
ssm.def("nip", &GooseFEM::Element::Quad4::Nodal::nip, "Return number of integration point" );
ssm.def("xi" , &GooseFEM::Element::Quad4::Nodal::xi , "Return integration point coordinates");
ssm.def("w" , &GooseFEM::Element::Quad4::Nodal::w , "Return integration point weights" );
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
{
py::module sm = mElement.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)");
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Element::Hex8::Quadrature>(sm, "Quadrature")
.def(py::init<const xt::xtensor<double,3> &>(), "Quadrature", py::arg("x"))
.def(py::init<const xt::xtensor<double,3> &, const xt::xtensor<double,2> &, const xt::xtensor<double,1> &>(), "Quadrature", py::arg("x"), py::arg("xi"), py::arg("w"))
.def("nelem", &GooseFEM::Element::Hex8::Quadrature::nelem, "Return number of elements" )
.def("nne" , &GooseFEM::Element::Hex8::Quadrature::nne , "Return number of nodes per element" )
.def("ndim" , &GooseFEM::Element::Hex8::Quadrature::ndim , "Return number of dimensions" )
.def("nip" , &GooseFEM::Element::Hex8::Quadrature::nip , "Return number of integration points")
.def("dV" , py::overload_cast<>(&GooseFEM::Element::Hex8::Quadrature::dV , py::const_), "Integration point volume (qscalar)")
.def("dVtensor", py::overload_cast<>(&GooseFEM::Element::Hex8::Quadrature::dVtensor, py::const_), "Integration point volume (qtensor)")
.def("gradN_vector" , py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Hex8::Quadrature::gradN_vector , py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("gradN_vector_T" , py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Hex8::Quadrature::gradN_vector_T , py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("symGradN_vector", py::overload_cast<const xt::xtensor<double,3> &>(&GooseFEM::Element::Hex8::Quadrature::symGradN_vector, py::const_), "Dyadic product, returns 'qtensor'", py::arg("elemvec"))
.def("int_N_scalar_NT_dV" , py::overload_cast<const xt::xtensor<double,2> &>(&GooseFEM::Element::Hex8::Quadrature::int_N_scalar_NT_dV , py::const_), "Integration, returns 'elemmat'", py::arg("qscalar"))
.def("int_gradN_dot_tensor2_dV", py::overload_cast<const xt::xtensor<double,4> &>(&GooseFEM::Element::Hex8::Quadrature::int_gradN_dot_tensor2_dV, py::const_), "Integration, returns 'elemvec'", py::arg("qtensor"))
.def("__repr__", [](const GooseFEM::Element::Hex8::Quadrature &){ return "<GooseFEM.Element.Hex8.Quadrature>"; });
// -------------------------------------------------------------------------------------------------
{
py::module ssm = sm.def_submodule("Gauss", "Gauss quadrature");
ssm.def("nip", &GooseFEM::Element::Hex8::Gauss::nip, "Return number of integration point" );
ssm.def("xi" , &GooseFEM::Element::Hex8::Gauss::xi , "Return integration point coordinates");
ssm.def("w" , &GooseFEM::Element::Hex8::Gauss::w , "Return integration point weights" );
}
// -------------------------------------------------------------------------------------------------
{
py::module ssm = sm.def_submodule("Nodal", "Nodal quadrature");
ssm.def("nip", &GooseFEM::Element::Hex8::Nodal::nip, "Return number of integration point" );
ssm.def("xi" , &GooseFEM::Element::Hex8::Nodal::xi , "Return integration point coordinates");
ssm.def("w" , &GooseFEM::Element::Hex8::Nodal::w , "Return integration point weights" );
}
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
py::module mMesh = m.def_submodule("Mesh", "Generic mesh routines");
// -------------------------------------------------------------------------------------------------
mMesh.def("dofs", &GooseFEM::Mesh::dofs, "List with DOF-numbers (in sequential order)", py::arg("nnode"), py::arg("ndim"));
mMesh.def("renumber", &GooseFEM::Mesh::renumber, "Renumber DOF-list to use the lowest possible index", py::arg("dofs"));
mMesh.def("renumber_index", &GooseFEM::Mesh::renumber_index, "Index-list to renumber", py::arg("dofs"));
mMesh.def("reorder", &GooseFEM::Mesh::reorder, "Renumber DOF-list to begin or end with 'idx'", py::arg("dofs"), py::arg("idx"), py::arg("location")="end");
mMesh.def("reorder_index", &GooseFEM::Mesh::reorder_index, "Index-list to reorder", py::arg("dofs"), py::arg("idx"), py::arg("location")="end");
mMesh.def("coordination", &GooseFEM::Mesh::coordination, "Coordination number of each node", py::arg("conn"));
mMesh.def("elem2node", &GooseFEM::Mesh::elem2node, "Elements connect to each node", py::arg("conn"));
// =================================================================================================
{
py::module sm = mMesh.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)");
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Mesh::Hex8::Regular>(sm, "Regular")
.def(py::init<size_t,size_t,size_t,double>(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1.)
.def("nelem" , &GooseFEM::Mesh::Hex8::Regular::nelem )
.def("nnode" , &GooseFEM::Mesh::Hex8::Regular::nnode )
.def("nne" , &GooseFEM::Mesh::Hex8::Regular::nne )
.def("ndim" , &GooseFEM::Mesh::Hex8::Regular::ndim )
.def("coor" , &GooseFEM::Mesh::Hex8::Regular::coor )
.def("conn" , &GooseFEM::Mesh::Hex8::Regular::conn )
.def("nodesFront" , &GooseFEM::Mesh::Hex8::Regular::nodesFront )
.def("nodesBack" , &GooseFEM::Mesh::Hex8::Regular::nodesBack )
.def("nodesLeft" , &GooseFEM::Mesh::Hex8::Regular::nodesLeft )
.def("nodesRight" , &GooseFEM::Mesh::Hex8::Regular::nodesRight )
.def("nodesBottom" , &GooseFEM::Mesh::Hex8::Regular::nodesBottom )
.def("nodesTop" , &GooseFEM::Mesh::Hex8::Regular::nodesTop )
.def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace )
.def("nodesBackFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBackFace )
.def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace )
.def("nodesRightFace" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFace )
.def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace )
.def("nodesTopFace" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFace )
.def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge )
.def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge )
.def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge )
.def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge )
.def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge )
.def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge )
.def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge )
.def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge )
.def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge )
.def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge )
.def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge )
.def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge )
.def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge )
.def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge )
.def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge )
.def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge )
.def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge )
.def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge )
.def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge )
.def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge )
.def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge )
.def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge )
.def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge )
.def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge )
.def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge )
.def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge )
.def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge )
.def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge )
.def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge )
.def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge )
.def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge )
.def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge )
.def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge )
.def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge )
.def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge )
.def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge )
.def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge )
.def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge )
.def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge )
.def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge )
.def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge )
.def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge )
.def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge )
.def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge )
.def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge )
.def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge )
.def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge )
.def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge )
.def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner )
.def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner)
.def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner )
.def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner )
.def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner )
.def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner )
.def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner )
.def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner )
.def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner )
.def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner )
.def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner )
.def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner )
.def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner )
.def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner)
.def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner)
.def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner)
.def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner)
.def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner)
.def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner )
.def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner )
.def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner )
.def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner )
.def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner )
.def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner )
.def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner )
.def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner )
.def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner )
.def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner )
.def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner )
.def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner )
.def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner )
.def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner )
.def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner )
.def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner )
.def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner )
.def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner )
.def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner )
.def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner )
.def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner )
.def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner )
.def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner )
.def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner )
.def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner )
.def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner )
.def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner )
.def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner )
.def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner )
.def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Hex8::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Hex8::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic )
.def("__repr__", [](const GooseFEM::Mesh::Hex8::Regular &){ return "<GooseFEM.Mesh.Hex8.Regular>"; });
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Mesh::Hex8::FineLayer>(sm, "FineLayer")
.def(py::init<size_t,size_t,size_t,double,size_t>(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1., py::arg("nfine")=1)
.def("nelem" , &GooseFEM::Mesh::Hex8::FineLayer::nelem )
.def("nnode" , &GooseFEM::Mesh::Hex8::FineLayer::nnode )
.def("nne" , &GooseFEM::Mesh::Hex8::FineLayer::nne )
.def("ndim" , &GooseFEM::Mesh::Hex8::FineLayer::ndim )
.def("shape" , &GooseFEM::Mesh::Hex8::FineLayer::shape )
.def("coor" , &GooseFEM::Mesh::Hex8::FineLayer::coor )
.def("conn" , &GooseFEM::Mesh::Hex8::FineLayer::conn )
.def("elementsMiddleLayer" , &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer )
.def("nodesFront" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFront )
.def("nodesBack" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBack )
.def("nodesLeft" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft )
.def("nodesRight" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRight )
.def("nodesBottom" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom )
.def("nodesTop" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTop )
.def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace )
.def("nodesBackFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace )
.def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace )
.def("nodesRightFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace )
.def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace )
.def("nodesTopFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace )
.def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge )
.def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge )
.def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge )
.def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge )
.def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge )
.def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge )
.def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge )
.def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge )
.def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge )
.def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge )
.def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge )
.def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge )
.def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge )
.def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge )
.def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge )
.def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge )
.def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge )
.def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge )
.def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge )
.def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge )
.def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge )
.def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge )
.def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge )
.def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge )
.def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge )
.def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge )
.def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge )
.def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge )
.def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge )
.def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge )
.def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge )
.def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge )
.def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge )
.def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge )
.def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge )
.def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge )
.def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge )
.def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge )
.def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge )
.def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge )
.def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge )
.def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge )
.def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge )
.def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge )
.def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge )
.def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge )
.def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge )
.def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge )
.def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner )
.def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner)
.def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner )
.def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner )
.def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner )
.def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner )
.def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner )
.def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner )
.def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner )
.def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner )
.def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner )
.def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner )
.def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner )
.def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner)
.def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner)
.def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner)
.def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner)
.def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner)
.def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner )
.def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner )
.def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner )
.def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner )
.def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner )
.def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner )
.def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner )
.def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner )
.def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner )
.def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner )
.def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner )
.def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner )
.def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner )
.def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner )
.def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner )
.def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner )
.def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner )
.def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner )
.def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner )
.def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner )
.def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner )
.def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner )
.def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner )
.def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner )
.def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner )
.def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner )
.def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner )
.def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner )
.def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner )
.def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Hex8::FineLayer::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic )
.def("__repr__", [](const GooseFEM::Mesh::Hex8::FineLayer &){ return "<GooseFEM.Mesh.Hex8.FineLayer>"; });
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
{
py::module sm = mMesh.def_submodule("Quad4", "Linear quadrilateral elements (2D)");
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Mesh::Quad4::Regular>(sm, "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("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("nodesPeriodic" , &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Quad4::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Quad4::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic )
.def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular &){ return "<GooseFEM.Mesh.Quad4.Regular>"; });
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Mesh::Quad4::FineLayer>(sm, "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("shape" , &GooseFEM::Mesh::Quad4::FineLayer::shape )
.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("elementsMiddleLayer" , &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer )
.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("nodesPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Quad4::FineLayer::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic )
.def("__repr__",
[](const GooseFEM::Mesh::Quad4::FineLayer &){ return "<GooseFEM.Mesh.Quad4.FineLayer>"; }
);
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
{
py::module sm = mMesh.def_submodule("Tri3" , "Linear triangular elements (2D)");
// -------------------------------------------------------------------------------------------------
py::class_<GooseFEM::Mesh::Tri3::Regular>(sm, "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::Tri3::Regular::coor )
.def("conn" , &GooseFEM::Mesh::Tri3::Regular::conn )
.def("nelem" , &GooseFEM::Mesh::Tri3::Regular::nelem )
.def("nnode" , &GooseFEM::Mesh::Tri3::Regular::nnode )
.def("nne" , &GooseFEM::Mesh::Tri3::Regular::nne )
.def("ndim" , &GooseFEM::Mesh::Tri3::Regular::ndim )
.def("nodesBottomEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge )
.def("nodesTopEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge )
.def("nodesLeftEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge )
.def("nodesRightEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge )
.def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge )
.def("nodesTopOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge )
.def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge )
.def("nodesRightOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge )
.def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner )
.def("nodesBottomRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner)
.def("nodesTopLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner )
.def("nodesTopRightCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner )
.def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner )
.def("nodesLeftTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner )
.def("nodesRightBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner)
.def("nodesRightTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner )
.def("nodesPeriodic" , &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic )
.def("nodesOrigin" , &GooseFEM::Mesh::Tri3::Regular::nodesOrigin )
.def("dofs" , &GooseFEM::Mesh::Tri3::Regular::dofs )
.def("dofsPeriodic" , &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic )
.def("__repr__", [](const GooseFEM::Mesh::Tri3::Regular &){ return "<GooseFEM.Mesh.Tri3.Regular>"; });
// -------------------------------------------------------------------------------------------------
sm.def("getOrientation", &GooseFEM::Mesh::Tri3::getOrientation, "Get the orientation of each element", py::arg("coor"), py::arg("conn"));
sm.def("retriangulate", &GooseFEM::Mesh::Tri3::retriangulate, "Re-triangulate existing mesh", py::arg("coor"), py::arg("conn"), py::arg("orientation")=-1);
// -------------------------------------------------------------------------------------------------
}
// =================================================================================================
}

Event Timeline

c4science · Help