forked from 1b-agent/fieldcraft-upstream
verification: activate the two shell benchmark cases (#106)
The Scordelis-Lo roof and the pinched cylinder had their published references cited but no recipe, because mesh.block's 2D form is planar and cannot build a curved patch. Both are now active: the mesh is written node by node (mesh.node.add / mesh.element.add, the route docs/verification.md allows), the same cylinder patch the solver's own benchmarks.rs builds in code, with the quarter/eighth symmetry planes, the rigid diaphragm and the published load. - scordelis-lo-roof-shell: 16x16 flat quad4 (MITC4) shells, R = 25, t = 0.25, self-weight 90 per area; u_z at the free-edge mid-span -0.3024 within 2% (ours -0.29901009815961804, 1.12% low). - pinched-cylinder: 32x32 flat quad4 shells, R = 300, t = 3, a quarter of the unit pinch; u under the load -1.8248e-5 within 2% (ours -1.803845342860719e-05, 1.15% low). Each README carries the generator (so the recipe is reproducible), the refinement study the tolerance reason cites, and how the reference was cited before. pinched-hemisphere and mh-twisted-beam now say what they actually wait for: a doubly curved or helical shell mesh written node by node, not missing shell elements.
This commit is contained in:
Erkin Alp Güney authored and Irakli Kourtis committed
1 parent
a0c7d6ce9d
commit
18be60672b
12 files changed
+2973
-18
No files matched your search
@@ -1,3 +1,115 @@
|
||||
# Pinched cylinder with rigid diaphragms
|
||||
|
||||
Radius 300, length 600, thickness 3, E = 3 × 10⁶, ν = 0.3, rigid end diaphragms, two opposite unit loads at mid-length.
|
||||
Our model of the published problem: a cylinder of radius 300, length 600 and thickness 3,
|
||||
E = 3 × 10⁶, ν = 0.3, capped by rigid diaphragms (u_y = u_z = rot_x = 0 at the ends) and pinched by
|
||||
two opposite unit forces at mid-length, one on each of two diametrically opposite generators.
|
||||
Reference: 1.8248 × 10⁻⁵ radial displacement under the load (Belytschko, Stolarski, Liu, Carpenter
|
||||
and Ong, *Stress projection for membrane and shear locking in shell finite elements*, CMAME 51
|
||||
(1985) 221–258, the shell obstacle course) — the value the case already cited.
|
||||
|
||||
The eighth model: x from the loaded section (x = 0, symmetry in x) to the diaphragm (x = 300), φ
|
||||
from the loaded generator (φ = 0, symmetry in y) to the side (φ = 90°, symmetry in z), with
|
||||
R sin φ = y and R cos φ = z on the mid-surface R = 300. A quarter of the unit force, −0.25 in z,
|
||||
acts on the single loaded node (0, 0, 300) — the model the solver's own pinched-cylinder benchmark
|
||||
builds in code.
|
||||
|
||||
## The mesh
|
||||
|
||||
`mesh.block`'s 2D form takes its four corners at z = 0, so it cannot mesh a patch on a cylindrical
|
||||
surface (it can be curved in its own plane, but it is planar): the recipe writes the 1089 nodes and
|
||||
1024 quad4 elements one by one with `mesh.node.add` and `mesh.element.add`, the route
|
||||
[docs/verification.md](../../../docs/verification.md) allows ("nodes and elements written one
|
||||
by one"). Node 1 + i + 33 j lies at
|
||||
(300 i/32, 300 sin(90° j/32), 300 cos(90° j/32)); element 1 + i + 32 j is
|
||||
[node(i, j), node(i + 1, j), node(i + 1, j + 1), node(i, j + 1)], counter-clockwise seen from
|
||||
outside, so its normal points away from the axis — the same mesh the solver's own pinched-cylinder
|
||||
benchmark builds in code
|
||||
([crates/solver/tests/benchmarks.rs](../../../crates/solver/tests/benchmarks.rs)).
|
||||
|
||||
The recipe was generated with the same script as `scordelis-lo-roof-shell`
|
||||
(`python3 cylinder_recipe.py pinched-cylinder 32`):
|
||||
|
||||
```python
|
||||
import json, math
|
||||
|
||||
def pinched_cylinder(n):
|
||||
label = lambda i, j: 1 + i + (n + 1) * j
|
||||
commands = []
|
||||
for j in range(n + 1):
|
||||
for i in range(n + 1):
|
||||
x, phi = 300.0 * i / n, (math.pi / 2) * j / n
|
||||
commands.append({"command": "mesh.node.add",
|
||||
"params": {"label": label(i, j),
|
||||
"xyz": [x, 300.0 * math.sin(phi), 300.0 * math.cos(phi)]}})
|
||||
for j in range(n):
|
||||
for i in range(n):
|
||||
nodes = [label(i, j), label(i + 1, j), label(i + 1, j + 1), label(i, j + 1)]
|
||||
commands.append({"command": "mesh.element.add",
|
||||
"params": {"label": 1 + i + n * j, "kind": "quad4", "nodes": nodes}})
|
||||
return commands + [
|
||||
{"command": "model.set.create", "params": {"name": "shell", "elements": list(range(1, n * n + 1))}},
|
||||
{"command": "model.set.create", "params": {"name": "diaphragm", "nodes": [label(n, j) for j in range(n + 1)]}},
|
||||
{"command": "model.set.create", "params": {"name": "middle", "nodes": [label(0, j) for j in range(n + 1)]}},
|
||||
{"command": "model.set.create", "params": {"name": "top", "nodes": [label(i, 0) for i in range(n + 1)]}},
|
||||
{"command": "model.set.create", "params": {"name": "side", "nodes": [label(i, n) for i in range(n + 1)]}},
|
||||
{"command": "model.set.create", "params": {"name": "pinch", "nodes": [label(0, 0)]}},
|
||||
{"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 3e6, "poisson": 0.3}}},
|
||||
{"command": "model.section.create", "params": {"name": "shell", "set": "shell", "material": "Mat", "kind": {"type": "shell", "thickness": 3.0}}},
|
||||
{"command": "model.boundary_condition.create", "params": {"name": "diaphragm", "type": "displacement", "target": {"set": "diaphragm"}, "values": [None, 0, 0, 0, None, None]}},
|
||||
{"command": "model.boundary_condition.create", "params": {"name": "middle", "type": "symmetry", "target": {"set": "middle"}, "plane": "x"}},
|
||||
{"command": "model.boundary_condition.create", "params": {"name": "top", "type": "symmetry", "target": {"set": "top"}, "plane": "y"}},
|
||||
{"command": "model.boundary_condition.create", "params": {"name": "side", "type": "symmetry", "target": {"set": "side"}, "plane": "z"}},
|
||||
{"command": "model.load.create", "params": {"name": "pinch", "type": "force", "target": {"set": "pinch"}, "vector": [0, 0, -0.25]}},
|
||||
{"command": "model.step.create", "params": {"name": "Load", "procedure": {"type": "static"}, "boundary_conditions": ["diaphragm", "middle", "top", "side"], "loads": ["pinch"]}},
|
||||
{"command": "solve.static", "params": {}, "id": "solve"},
|
||||
]
|
||||
|
||||
print("[\n " + ",\n ".join(json.dumps(c) for c in pinched_cylinder(32)) + "\n]")
|
||||
```
|
||||
|
||||
## A fixed-mesh benchmark, not a convergence study
|
||||
|
||||
The quantity is the deflection *under a point load*, and there a Reissner–Mindlin (MITC4) shell
|
||||
carries an unbounded transverse shear force: the deflection under the load keeps growing with
|
||||
refinement instead of converging. The published 1.8248 × 10⁻⁵ is a thin-shell value, and this case
|
||||
compares against it on the conventional 32 × 32 eighth mesh — the mesh the solver's own
|
||||
pinched-cylinder benchmark uses — not as a limit. The table shows what refinement actually does
|
||||
(each row is the listing above with the `32` in its last line replaced by the row's divisions and
|
||||
run through `cargo xtask oracle` as a scratch case):
|
||||
|
||||
| cells | u | vs reference |
|
||||
|---|---|---|
|
||||
| 8 × 8 | −1.3559 × 10⁻⁵ | −25.7% |
|
||||
| 16 × 16 | −1.6942 × 10⁻⁵ | −7.2% |
|
||||
| 24 × 24 | −1.7735 × 10⁻⁵ | −2.8% |
|
||||
| 32 × 32 (this case) | −1.803845342860719 × 10⁻⁵ | −1.15% |
|
||||
| 48 × 48 | −1.8275 × 10⁻⁵ | +0.15% |
|
||||
| 64 × 64 | −1.8369 × 10⁻⁵ | +0.66% |
|
||||
|
||||
The sequence rises through the reference between 32 × 32 and 48 × 48 and keeps rising at the two
|
||||
finer meshes, which is what a shear-deformable shell under a point load must do.
|
||||
|
||||
The solver's in-code benchmark of the same problem shows the same ascent on its 8 × 8, 16 × 16 and
|
||||
32 × 32 eighth meshes (it gates 32 × 32 within 2%), so the command route and the in-code route
|
||||
agree on the fixed-mesh comparison.
|
||||
|
||||
## Tolerance
|
||||
|
||||
2% covers the agreement with the thin-shell reference at the 32 × 32 eighth mesh (1.15% low), the
|
||||
conventional mesh of the shell obstacle course and the one the solver's own pinched-cylinder
|
||||
benchmark gates at 2%, so the command route and the in-code route are compared on equal terms. The
|
||||
reference itself is quoted by Belytschko et al. to four figures (1.8248 × 10⁻⁵). Because the point
|
||||
load makes the deflection mesh-dependent rather than convergent (above), this is a fixed-mesh
|
||||
tolerance and not a convergence statement: the two finer meshes sit 0.15% and 0.66% above the
|
||||
reference and stay inside it. The mesh is not marked `slow`: the case solves in about 0.15 s.
|
||||
|
||||
`fixture.json` is the accepted value at 32 × 32 (−1.803845342860719 × 10⁻⁵) with the runner's
|
||||
default band: the mesh is written by the recipe, not generated from geometry, so it is bit-identical
|
||||
on every platform and the result is frozen.
|
||||
|
||||
## Running it
|
||||
|
||||
```sh
|
||||
cargo xtask oracle --case pinched-cylinder
|
||||
cargo xtask oracle --bless pinched-cylinder # only after the case passes
|
||||
```
|
||||
@@ -2,13 +2,29 @@
|
||||
"title": "Pinched cylinder with rigid diaphragms",
|
||||
"kind": "benchmark",
|
||||
"source": "Belytschko, Stolarski, Liu, Carpenter and Ong, 'Stress projection for membrane and shear locking in shell finite elements', Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258 (the shell obstacle course); reference 1.8248 × 10⁻⁵ under the load",
|
||||
"status": "pending",
|
||||
"requires": "shell elements (M8)",
|
||||
"quantities": [
|
||||
{
|
||||
"id": "u",
|
||||
"quantity": "radial displacement under the load",
|
||||
"value": -1.8248e-05
|
||||
"quantity": "radial displacement under the load, mid-surface",
|
||||
"measure": {
|
||||
"command": "solve.probe",
|
||||
"params": {
|
||||
"step": "Load",
|
||||
"at": [
|
||||
[
|
||||
0,
|
||||
0,
|
||||
300
|
||||
]
|
||||
]
|
||||
},
|
||||
"pointer": "/nodes/0/displacement/2"
|
||||
},
|
||||
"value": -1.8248e-05,
|
||||
"tolerance": {
|
||||
"relative": 0.02
|
||||
},
|
||||
"reason": "a fixed-mesh comparison at the conventional 32 × 32 eighth mesh, the mesh the solver's own in-code benchmark uses: the reference is the thin-shell value Belytschko et al. quote to four figures (1.8248 × 10⁻⁵), while the deflection under a point load has no mesh-converged limit in a Reissner–Mindlin (MITC4) shell, whose transverse shear force is unbounded at the load, so the value keeps rising with refinement (8 × 8, 16 × 16, 24 × 24, 32 × 32, 48 × 48 and 64 × 64 give 1.3559, 1.6942, 1.7735, 1.8038, 1.8275 and 1.8369 × 10⁻⁵, crossing the reference between 32 × 32 and 48 × 48 and staying above it); 2% is the margin the in-code benchmark gates 32 × 32 with, applied here to the same mesh (crates/solver/tests/benchmarks.rs, README)"
|
||||
}
|
||||
]
|
||||
}
|
||||
@@ -0,0 +1,9 @@
|
||||
{
|
||||
"note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.",
|
||||
"values": {
|
||||
"u": {
|
||||
"band": 1.8247999999999998e-14,
|
||||
"value": -0.00001803845342860719
|
||||
}
|
||||
}
|
||||
}
|
||||
File diff suppressed because it is too large.
Load diff
Reference in new issue
Block a user