OSCR

Sparse polynomial surrogates for F-actin networks with compliant crosslinkers.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches
  1. [1] § Methods › Homogenization into an affine network › Kinematics and strain-energy functions ↔ generate.py, lines 993–1034 · score 0.62 · Neo Hookean, Mooney Rivlin, oriented, models, network
  2. [2] § Methods › Homogenization into an affine network › Kinematics and strain-energy functions ↔ src/mod_hyperelastic.f90, lines 1–58 · score 0.55 · Mooney Rivlin, Hookean, hyperelastic, Neo, isochoric, Kinematics
  3. [3] § Methods › F-actin with compliant crosslinkers ↔ src/mod_network.f90, lines 1–62 · score 0.50 · Filament force, bending, persistence, Brent, stiffness, contour

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,154 lines · 40 KB · no license · 1 match

  1. #!/usr/bin/env python3
  2. """Generate a self-contained material law directory from a JSON configuration.
  4. Usage:
  5. python generate.py config.json # Generate from config
  6. python generate.py --example neo_hooke # Generate example config + material
  7. python generate.py --list # List available model types
  8. """
  9. import argparse
  10. import json
  11. import os
  12. import sys
  13. import stat
  14. from pathlib import Path
  16. SCRIPT_DIR = Path(__file__).resolve().parent
  17. SRC_DIR = SCRIPT_DIR / "src"
  19. # Source files in concatenation order
  20. SOURCE_FILES = [
  21. "mod_constants.f90",
  22. "mod_tensor.f90",
  23. "mod_kinematics.f90",
  24. "mod_continuum.f90",
  25. "mod_hyperelastic.f90",
  26. "mod_icosahedron.f90",
  27. "mod_anisotropic.f90",
  28. "mod_network.f90",
  29. "mod_damage.f90",
  30. "mod_viscosity.f90",
  31. "umat_builder.f90",
  32. "uexternaldb.f90",
  33. ]
  35. # Element-layer sources appended after SOURCE_FILES when emitting uel.f90
  36. UEL_SOURCE_FILES = [
  37. "element/mod_uel_config.f90",
  38. "element/mod_uel_shape.f90",
  39. "element/mod_uel_element.f90",
  40. "element/uel_entry.f90",
  41. ]
  43. # Defaults for the optional "element" config block (UEL emission)
  44. DEFAULT_ELEMENT = {
  45. "type": "u3d8", # 8-node brick (only type currently)
  46. "nint": 8, # volume integration points (1 or 8)
  47. "fbar": True, # F-bar locking treatment (active on the 8-pt brick)
  48. "num_elem": 1, # UEL elements in the real mesh (sizes globalSdv)
  49. "elem_offset": 1000, # dummy-mesh element-number offset (must match deck)
  50. }
  52. # --- Example configurations ---------------------------------------------------
  54. EXAMPLES = {
  55. "neo_hooke": {
  56. "name": "neo_hooke",
  57. "kbulk": 1000.0,
  58. "iso_type": 1, "iso_params": [10.0],
  59. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  60. "network_type": 0, "network_params": [],
  61. "damage_type": 0, "damage_params": [],
  62. "n_visco": 0, "visco_params": [],
  63. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  64. },
  65. "mooney_rivlin": {
  66. "name": "mooney_rivlin",
  67. "kbulk": 1000.0,
  68. "iso_type": 2, "iso_params": [6.3, 0.012],
  69. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  70. "network_type": 0, "network_params": [],
  71. "damage_type": 0, "damage_params": [],
  72. "n_visco": 0, "visco_params": [],
  73. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  74. },
  75. "humphrey_hgo": {
  76. "name": "humphrey_hgo",
  77. "kbulk": 500.0,
  78. "iso_type": 4, "iso_params": [2.0, 1.5],
  79. "aniso_type": 1, "n_fiber_fam": 1,
  80. "aniso_params": [100.0, 10.0, 0.226, 1.0, 0.0, 0.0],
  81. "network_type": 0, "network_params": [],
  82. "damage_type": 0, "damage_params": [],
  83. "n_visco": 0, "visco_params": [],
  84. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  85. },
  86. "ogden_3term": {
  87. "name": "ogden_3term",
  88. "kbulk": 1000.0,
  89. "iso_type": 3,
  90. "iso_params": [3, 1.3, 5.0, 0.5, -2.0, 0.012, 2.0],
  91. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  92. "network_type": 0, "network_params": [],
  93. "damage_type": 0, "damage_params": [],
  94. "n_visco": 0, "visco_params": [],
  95. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  96. },
  97. "neo_hooke_damage": {
  98. "name": "neo_hooke_damage",
  99. "kbulk": 1000.0,
  100. "iso_type": 1, "iso_params": [10.0],
  101. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  102. "network_type": 0, "network_params": [],
  103. "damage_type": 1, "damage_params": [5.0, 50.0],
  104. "n_visco": 0, "visco_params": [],
  105. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  106. },
  107. "neo_hooke_visco": {
  108. "name": "neo_hooke_visco",
  109. "kbulk": 1000.0,
  110. "iso_type": 1, "iso_params": [10.0],
  111. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  112. "network_type": 0, "network_params": [],
  113. "damage_type": 0, "damage_params": [],
  114. "n_visco": 1, "visco_params": [0.5, 0.25],
  115. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  116. },
  117. "affine_network": {
  118. "name": "affine_network",
  119. "kbulk": 500.0,
  120. "iso_type": 0, "iso_params": [],
  121. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  122. "network_type": 5, "network_params": [
  123. 0.5, 1.0e6, 2.0, 1.0, 6, 1.0, 0.0, 0.0,
  124. 1.0, 0.1, 0.01, 2.0, 0.1, 1.0, 0.0, 0.0,
  125. ],
  126. "damage_type": 0, "damage_params": [],
  127. "n_visco": 0, "visco_params": [],
  128. "test": {"stretch_max": 1.3, "gamma_max": 0.3, "nsteps": 200, "dtime": 0.01},
  129. },
  130. "humphrey_fiber": {
  131. "name": "humphrey_fiber",
  132. "kbulk": 500.0,
  133. "iso_type": 4, "iso_params": [2.0, 1.5],
  134. "aniso_type": 2, "n_fiber_fam": 1,
  135. "aniso_params": [100.0, 10.0, 1.0, 0.0, 0.0],
  136. "network_type": 0, "network_params": [],
  137. "damage_type": 0, "damage_params": [],
  138. "n_visco": 0, "visco_params": [],
  139. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  140. },
  141. "humphrey_hgo_damage": {
  142. "name": "humphrey_hgo_damage",
  143. "kbulk": 500.0,
  144. "iso_type": 4, "iso_params": [2.0, 1.5],
  145. "aniso_type": 1, "n_fiber_fam": 1,
  146. "aniso_params": [100.0, 10.0, 0.226, 1.0, 0.0, 0.0],
  147. "network_type": 0, "network_params": [],
  148. "damage_type": 1, "damage_params": [5.0, 50.0],
  149. "n_visco": 0, "visco_params": [],
  150. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  151. },
  152. "humphrey_fiber_damage": {
  153. "name": "humphrey_fiber_damage",
  154. "kbulk": 500.0,
  155. "iso_type": 4, "iso_params": [2.0, 1.5],
  156. "aniso_type": 2, "n_fiber_fam": 1,
  157. "aniso_params": [100.0, 10.0, 1.0, 0.0, 0.0],
  158. "network_type": 0, "network_params": [],
  159. "damage_type": 1, "damage_params": [5.0, 50.0],
  160. "n_visco": 0, "visco_params": [],
  161. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  162. },
  163. "mooney_rivlin_visco": {
  164. "name": "mooney_rivlin_visco",
  165. "kbulk": 1000.0,
  166. "iso_type": 2, "iso_params": [6.3, 0.012],
  167. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  168. "network_type": 0, "network_params": [],
  169. "damage_type": 0, "damage_params": [],
  170. "n_visco": 1, "visco_params": [0.5, 0.25],
  171. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  172. },
  173. "ogden_visco": {
  174. "name": "ogden_visco",
  175. "kbulk": 1000.0,
  176. "iso_type": 3,
  177. "iso_params": [3, 1.3, 5.0, 0.5, -2.0, 0.012, 2.0],
  178. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  179. "network_type": 0, "network_params": [],
  180. "damage_type": 0, "damage_params": [],
  181. "n_visco": 1, "visco_params": [0.5, 0.25],
  182. "test": {"stretch_max": 1.5, "gamma_max": 0.6, "nsteps": 400, "dtime": 0.01},
  183. },
  184. "humphrey_hgo_visco": {
  185. "name": "humphrey_hgo_visco",
  186. "kbulk": 500.0,
  187. "iso_type": 4, "iso_params": [2.0, 1.5],
  188. "aniso_type": 1, "n_fiber_fam": 1,
  189. "aniso_params": [100.0, 10.0, 0.226, 1.0, 0.0, 0.0],
  190. "network_type": 0, "network_params": [],
  191. "damage_type": 0, "damage_params": [],
  192. "n_visco": 1, "visco_params": [0.5, 0.25],
  193. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  194. },
  195. "humphrey_fiber_visco": {
  196. "name": "humphrey_fiber_visco",
  197. "kbulk": 500.0,
  198. "iso_type": 4, "iso_params": [2.0, 1.5],
  199. "aniso_type": 2, "n_fiber_fam": 1,
  200. "aniso_params": [100.0, 10.0, 1.0, 0.0, 0.0],
  201. "network_type": 0, "network_params": [],
  202. "damage_type": 0, "damage_params": [],
  203. "n_visco": 1, "visco_params": [0.5, 0.25],
  204. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  205. },
  206. "contractile_network": {
  207. "name": "contractile_network",
  208. "kbulk": 500.0,
  209. "iso_type": 0, "iso_params": [],
  210. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  211. "network_type": 4, "network_params": [
  212. 0.5, # PHI
  213. 0.2, 2.0, 1.0, # N, B_orient, EFI
  214. 11.0, 11.0, # FRIC, FFMAX
  215. 6, 1.0, 0.0, 0.0, # factor, prefdir
  216. 0.988, 0.804, 38600.0, 0.438, # L, R0F, mu0, beta
  217. 0.065, 1.007, # B0, lambda0
  218. 0.014, 0.667, # R0C, ETAC
  219. 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, # KCH(7) rate constants
  220. ],
  221. "damage_type": 0, "damage_params": [],
  222. "n_visco": 0, "visco_params": [],
  223. "test": {"stretch_max": 1.2, "gamma_max": 0.2, "nsteps": 100, "dtime": 0.01},
  224. },
  225. "mixed_network": {
  226. "name": "mixed_network",
  227. "kbulk": 500.0,
  228. "iso_type": 0, "iso_params": [],
  229. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  230. "network_type": 3, "network_params": [
  231. 0.5, # PHI
  232. 5.0e5, 2.0, # N_naff, PP
  233. 5.0e5, 2.0, 1.0, # N_aff, B_orient, EFI
  234. 6, 1.0, 0.0, 0.0, # factor, prefdir
  235. 1.0, 0.1, 0.01, 2.0, 0.1, 1.0, 0.0, 0.0, # filprops(8)
  236. ],
  237. "damage_type": 0, "damage_params": [],
  238. "n_visco": 0, "visco_params": [],
  239. "test": {"stretch_max": 1.3, "gamma_max": 0.3, "nsteps": 200, "dtime": 0.01},
  240. },
  241. "humphrey_hgo_ai": {
  242. "name": "humphrey_hgo_ai",
  243. "kbulk": 500.0,
  244. "iso_type": 4, "iso_params": [2.0, 1.5],
  245. "aniso_type": 3, "n_fiber_fam": 1,
  246. "aniso_params": [100.0, 10.0, 5.0, 6, 1.0, 0.0, 0.0],
  247. "network_type": 0, "network_params": [],
  248. "damage_type": 0, "damage_params": [],
  249. "n_visco": 0, "visco_params": [],
  250. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 200, "dtime": 0.01},
  251. },
  252. "affine_network_linkers": {
  253. "name": "affine_network_linkers",
  254. "kbulk": 500.0,
  255. "iso_type": 0, "iso_params": [],
  256. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  257. "network_type": 5, "network_params": [
  258. 0.5, 1.0e6, 2.0, 1.0, 6, 1.0, 0.0, 0.0,
  259. 1.0, 0.1, 0.01, 2.0, 0.1, 1.0, 0.014, 0.6667,
  260. ],
  261. "damage_type": 0, "damage_params": [],
  262. "n_visco": 0, "visco_params": [],
  263. "test": {"stretch_max": 1.3, "gamma_max": 0.3, "nsteps": 200, "dtime": 0.01},
  264. },
  265. "humphrey_muscle": {
  266. "name": "humphrey_muscle",
  267. "kbulk": 500.0,
  268. "iso_type": 4, "iso_params": [2.0, 1.5],
  269. "aniso_type": 5, "n_fiber_fam": 1,
  270. "aniso_params": [100.0, 10.0, 1.0, 0.0, 0.0, 50.0],
  271. "network_type": 0, "network_params": [],
  272. "damage_type": 0, "damage_params": [],
  273. "n_visco": 0, "visco_params": [],
  274. "test": {"stretch_max": 1.3, "gamma_max": 0.4, "nsteps": 400, "dtime": 0.01},
  275. },
  276. "nonaffine_network": {
  277. "name": "nonaffine_network",
  278. "kbulk": 500.0,
  279. "iso_type": 0, "iso_params": [],
  280. "aniso_type": 0, "n_fiber_fam": 0, "aniso_params": [],
  281. "network_type": 6, "network_params": [
  282. 0.5, 1.0e6, 2.0, 6,
  283. 1.0, 0.1, 0.01, 2.0, 0.1, 1.0, 0.0, 0.0,
  284. ],
  285. "damage_type": 0, "damage_params": [],
  286. "n_visco": 0, "visco_params": [],
  287. "test": {"stretch_max": 1.3, "gamma_max": 0.3, "nsteps": 200, "dtime": 0.01},
  288. },
  289. }
  291. # --- Helpers ------------------------------------------------------------------
  293. def compute_nstatev(cfg):
  294. n = 1
  295. if cfg["damage_type"] > 0:
  296. n += 2
  297. n += 9 * cfg["n_visco"]
  298. if cfg["network_type"] == 4: # contractile
  299. # Extract factor from network_params (position 6, 0-indexed)
  300. nparams = cfg["network_params"]
  301. factor = int(nparams[6])
  302. # Icosahedron: 20 faces, each subdivided into factor^2 subtriangles
  303. nwp = 20 * factor * factor
  304. n += 4 + nwp # FRAC(4) + RU0(nwp)
  305. return n
  308. def build_props(cfg):
  309. """Return flat list of all PROPS values."""
  310. p = [
  311. cfg["kbulk"],
  312. float(cfg["iso_type"]),
  313. float(cfg["aniso_type"]),
  314. float(cfg["n_fiber_fam"]),
  315. float(cfg["network_type"]),
  316. float(cfg["damage_type"]),
  317. float(cfg["n_visco"]),
  318. ]
  319. p.extend(cfg["iso_params"])
  320. for _ in range(cfg["n_fiber_fam"]):
  321. p.extend(cfg["aniso_params"])
  322. p.extend(cfg["network_params"])
  323. p.extend(cfg["damage_params"])
  324. p.extend(cfg["visco_params"])
  325. return p
  328. def fmt_props_fortran(props):
  329. """Generate Fortran PROPS assignment lines."""
  330. lines = []
  331. for i, v in enumerate(props, 1):
  332. if v == int(v) and abs(v) < 1e10:
  333. lines.append(f" props({i}) = {v:.1f}d0")
  334. else:
  335. lines.append(f" props({i}) = {v:g}d0")
  336. return "\n".join(lines)
  339. def fmt_props_abaqus(props):
  340. """Format PROPS values for ABAQUS *User Material card (8 per line)."""
  341. lines = []
  342. for i in range(0, len(props), 8):
  343. chunk = props[i : i + 8]
  344. line = ", ".join(f"{v:g}" for v in chunk)
  345. if i + 8 < len(props):
  346. line += ","
  347. lines.append(line)
  348. return "\n".join(lines)
  351. # --- File generators ----------------------------------------------------------
  353. def generate_umat_f90(outdir):
  354. """Concatenate all source modules into a single umat.f90."""
  355. out = outdir / "umat.f90"
  356. with open(out, "w") as f:
  357. f.write("! Auto-generated UMAT — do not edit manually.\n")
  358. f.write("! Regenerate with: python generate.py <config.json>\n\n")
  359. for src in SOURCE_FILES:
  360. path = SRC_DIR / src
  361. f.write(f"! {'='*72}\n")
  362. f.write(f"! SOURCE: {src}\n")
  363. f.write(f"! {'='*72}\n")
  364. f.write(path.read_text())
  365. f.write("\n\n")
  368. def generate_aba_param(outdir):
  369. src = SRC_DIR / "aba_param.inc"
  370. (outdir / "aba_param.inc").write_text(src.read_text())
  373. # --- UEL emission ---------------------------------------------------------
  375. def element_cfg(cfg):
  376. """Merged element block (or None if UEL emission is not requested)."""
  377. if "element" not in cfg:
  378. return None
  379. elem = dict(DEFAULT_ELEMENT)
  380. elem.update(cfg["element"] or {})
  381. return elem
  384. def subst_uel_config(text, elem):
  385. """Rewrite the parameter values in mod_uel_config.f90 from the element block."""
  386. import re
  387. text = re.sub(r"(numElem\s*=\s*)\d+", rf"\g<1>{elem['num_elem']}", text)
  388. text = re.sub(r"(ElemOffset\s*=\s*)\d+", rf"\g<1>{elem['elem_offset']}", text)
  389. fbar = ".true." if elem["fbar"] else ".false."
  390. text = re.sub(r"(use_fbar\s*=\s*)\.\w+\.", rf"\g<1>{fbar}", text)
  391. text = re.sub(r"(nIntPt\s*=\s*)\d+", rf"\g<1>{elem['nint']}", text)
  392. return text
  395. def uel_source(elem):
  396. """Concatenated uel.f90 text: material modules + element layer."""
  397. parts = []
  398. for src in SOURCE_FILES + UEL_SOURCE_FILES:
  399. text = (SRC_DIR / src).read_text()
  400. if src.endswith("mod_uel_config.f90"):
  401. text = subst_uel_config(text, elem)
  402. parts.append(f"! {'='*72}\n! SOURCE: {src}\n! {'='*72}\n" + text)
  403. return ("! Auto-generated UEL (user element + material library) — do not edit.\n"
  404. "! Regenerate with: python generate.py <config.json>\n\n"
  405. + "\n\n".join(parts) + "\n")
  408. def generate_uel_f90(outdir, cfg, elem):
  409. (outdir / "uel.f90").write_text(uel_source(elem))
  412. # Unit-cube nodes in the element's local ordering (sketch in mod_uel_element):
  413. # bottom face z=0: 1(0,0,0) 2(1,0,0) 3(1,1,0) 4(0,1,0); top face z=1: 5..8
  414. UEL_CUBE_COORDS = [
  415. (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 1.0, 0.0),
  416. (0.0, 0.0, 1.0), (1.0, 0.0, 1.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0),
  417. ]
  420. def generate_uel_test_driver(outdir, cfg, elem):
  421. """Standalone single-element driver: ramps the unit cube through affine
  422. deformation histories (same load cases as test_umat) and records the
  423. internal force on the x+ face (nodes 2,3,6,7)."""
  424. props = build_props(cfg)
  425. nprops = len(props)
  426. nstatev = compute_nstatev(cfg)
  427. nint = elem["nint"]
  428. test = cfg.get("test", {})
  429. stretch_max = test.get("stretch_max", 1.5)
  430. gamma_max = test.get("gamma_max", 0.6)
  431. nsteps = test.get("nsteps", 400)
  432. dtime = test.get("dtime", 0.01)
  434. coords_lines = "\n".join(
  435. f" coords(1,{i+1}) = {x:.1f}d0; coords(2,{i+1}) = {y:.1f}d0; coords(3,{i+1}) = {z:.1f}d0"
  436. for i, (x, y, z) in enumerate(UEL_CUBE_COORDS))
  438. code = f"""\
  439. ! Auto-generated UEL test driver for material: {cfg["name"]}
  440. ! Drives one U3D8 element through affine deformation ramps (no ABAQUS).
  441. ! Output columns: step time, control parameter, x+ face force (fx, fy, fz).
  442. program test_uel
  443. implicit none
  445. integer, parameter :: nnode = 8, ndofel = 24, mlvarx = 24, nrhs = 1
  446. integer, parameter :: nprops = {nprops}, nsdv = {nstatev}, nintp = {nint}
  447. integer, parameter :: nsvars = nintp*nsdv, njprop = 2
  448. integer, parameter :: mcrd = 3, jtype = 3, jelem = 1
  449. integer, parameter :: ndload = 0, mdload = 1, npredf = 1
  451. double precision :: rhs(mlvarx,1), amatrx(ndofel,ndofel), svars(nsvars)
  452. double precision :: energy(8), props(nprops), coords(mcrd,nnode)
  453. double precision :: u(ndofel), du(mlvarx,1), v(ndofel), a(ndofel)
  454. double precision :: time(2), dtime, params(1)
  455. double precision :: adlmag(mdload,1), ddlmag(mdload,1)
  456. double precision :: predef(2,npredf,nnode), pnewdt, period
  457. integer :: jdltyp(mdload,1), lflags(4), jprops(njprop)
  459. double precision :: Ftar(3,3)
  460. double precision :: stretch_max, gamma_max
  461. integer :: nsteps
  462. double precision, parameter :: zero = 0.0d0, one = 1.0d0
  464. ! --- Element/material setup ---
  465. {fmt_props_fortran(props)}
  466. jprops(1) = nsdv ! local SDVs per integration point
  467. jprops(2) = nsdv ! global SDVs per integration point (UVARM)
  468. lflags = 0
  469. lflags(1) = 1 ! static general step
  470. lflags(2) = 1 ! nlgeom=yes
  471. dtime = {dtime}d0
  472. nsteps = {nsteps}
  473. stretch_max = {stretch_max}d0
  474. gamma_max = {gamma_max}d0
  475. v = zero; a = zero; params = zero; energy = zero
  476. adlmag = zero; ddlmag = zero; predef = zero; jdltyp = 0
  477. period = zero
  479. {coords_lines}
  481. time = zero
  482. call uexternaldb(0, 0, time, zero, 0, 0)
  483. call execute_command_line('mkdir -p results')
  485. ! Uniaxial: F = diag(s, 1/sqrt(s), 1/sqrt(s))
  486. Ftar = zero
  487. Ftar(1,1) = stretch_max
  488. Ftar(2,2) = one/sqrt(stretch_max); Ftar(3,3) = one/sqrt(stretch_max)
  489. call run_case(Ftar, 'results/uel_uniaxial.dat', 'Uniaxial ')
  491. ! Biaxial: F = diag(s, s, 1/s^2)
  492. Ftar = zero
  493. Ftar(1,1) = stretch_max; Ftar(2,2) = stretch_max
  494. Ftar(3,3) = one/(stretch_max*stretch_max)
  495. call run_case(Ftar, 'results/uel_biaxial.dat', 'Biaxial ')
  497. ! Pure shear: F12 = F21 = gamma
  498. Ftar = zero; Ftar(1,1) = one; Ftar(2,2) = one; Ftar(3,3) = one
  499. Ftar(1,2) = gamma_max; Ftar(2,1) = gamma_max
  500. call run_case(Ftar, 'results/uel_shear.dat', 'Shear ')
  502. ! Simple shear: F12 = gamma
  503. Ftar = zero; Ftar(1,1) = one; Ftar(2,2) = one; Ftar(3,3) = one
  504. Ftar(1,2) = gamma_max
  505. call run_case(Ftar, 'results/uel_simple_shear.dat', 'Simple sh')
  507. contains
  509. !> Ramp the element from I to Ftar in nsteps affine increments,
  510. !> carrying svars (state) across increments.
  511. subroutine run_case(Ft, fname, label)
  512. double precision, intent(in) :: Ft(3,3)
  513. character(*), intent(in) :: fname, label
  514. double precision :: F(3,3), uold(ndofel), fface(3), s
  515. integer :: i, n, k, kk, face_nodes(4)
  517. face_nodes = (/2, 3, 6, 7/) ! x+ face (X=1)
  518. svars = zero; u = zero; uold = zero; time = zero
  519. pnewdt = one
  521. open(unit=21, file=fname, status='replace')
  522. do i = 1, nsteps
  523. s = dble(i)/dble(nsteps)
  524. F = identity3() + s*(Ft - identity3())
  525. ! Affine nodal displacements u_a = (F - I) X_a
  526. do n = 1, nnode
  527. do k = 1, 3
  528. u(3*(n-1)+k) = sum((F(k,:) - identity_row(k))*coords(:,n))
  529. end do
  530. end do
  531. du(:,1) = u - uold
  532. rhs = zero; amatrx = zero
  533. call uel(rhs, amatrx, svars, energy, ndofel, nrhs, nsvars, &
  534. props, nprops, coords, mcrd, nnode, u, du, v, a, jtype, &
  535. time, dtime, 1, i, jelem, params, ndload, jdltyp, adlmag, &
  536. predef, npredf, lflags, mlvarx, ddlmag, mdload, pnewdt, &
  537. jprops, njprop, period)
  538. ! Internal force on the x+ face: f = -sum(RHS) over face nodes
  539. fface = zero
  540. do kk = 1, 4
  541. n = face_nodes(kk)
  542. do k = 1, 3
  543. fface(k) = fface(k) - rhs(3*(n-1)+k, 1)
  544. end do
  545. end do
  546. write(21, '(5ES20.10)') time(1), s, fface(1), fface(2), fface(3)
  547. time(1) = time(1) + dtime
  548. uold = u
  549. end do
  550. close(21)
  551. write(*,'(A,A,A)') label, ' -> ', fname
  552. end subroutine run_case
  554. pure function identity3() result(iden)
  555. double precision :: iden(3,3)
  556. integer :: ii
  557. iden = 0.0d0
  558. do ii = 1, 3
  559. iden(ii,ii) = 1.0d0
  560. end do
  561. end function identity3
  563. pure function identity_row(k) result(row)
  564. integer, intent(in) :: k
  565. double precision :: row(3)
  566. row = 0.0d0
  567. row(k) = 1.0d0
  568. end function identity_row
  570. end program test_uel
  572. ! Stub for ABAQUS-provided routine (standalone builds only)
  573. subroutine getoutdir(outdir, lenoutdir)
  574. implicit none
  575. character(len=256), intent(out) :: outdir
  576. integer, intent(out) :: lenoutdir
  577. outdir = '.'
  578. lenoutdir = 1
  579. end subroutine getoutdir
  580. """
  581. (outdir / "test_uel.f90").write_text(code)
  584. def generate_uel_abaqus(outdir, cfg, elem):
  585. """ABAQUS single-element UEL deck: U3 real mesh + dummy mesh for UVARM
  586. visualization. Reuses the bcs_*.inp files written by generate_abaqus_dir
  587. (same node numbering and node sets)."""
  588. abq = outdir / "abaqus"
  589. abq.mkdir(exist_ok=True)
  591. props = build_props(cfg)
  592. nprops = len(props)
  593. nstatev = compute_nstatev(cfg)
  594. nvars = elem["nint"] * nstatev
  595. offset = elem["elem_offset"]
  596. dummy_type = "C3D8" if elem["nint"] == 8 else "C3D8R"
  598. # *UEL PROPERTY data: reals first, the two integer properties last
  599. uel_props = fmt_props_abaqus(list(props) + [nstatev, nstatev])
  601. deck = f"""\
  602. *Heading
  603. UEL single-element test — {cfg["name"]}
  604. ** Real mesh: user element U3 (8-node brick, F-bar={'on' if elem['fbar'] else 'off'}, {elem['nint']}-pt)
  605. ** Dummy mesh: {dummy_type} at element offset {offset}, carries UVARM output
  606. *Node, nset=all_nodes
  607. 1, 1., 1., 1.
  608. 2, 1., 0., 1.
  609. 3, 1., 1., 0.
  610. 4, 1., 0., 0.
  611. 5, 0., 1., 1.
  612. 6, 0., 0., 1.
  613. 7, 0., 1., 0.
  614. 8, 0., 0., 0.
  615. *User Element, type=U3, nodes=8, coordinates=3, properties={nprops}, iproperties=2, variables={nvars}, unsymm
  616. 1,2,3
  617. *Element, type=U3, elset=main_element
  618. 1, 5, 6, 8, 7, 1, 2, 4, 3
  619. *Element, type={dummy_type}, elset=dummy_mesh
  620. {1 + offset}, 5, 6, 8, 7, 1, 2, 4, 3
  621. *Nset, nset=Set-1, generate
  622. 2, 8, 2
  623. *Nset, nset=Set-2, generate
  624. 1, 7, 2
  625. *Nset, nset=Set-3
  626. 1, 2, 5, 6
  627. *Nset, nset=Set-4, generate
  628. 5, 8, 1
  629. *Nset, nset=Set-5
  630. 2, 4, 6, 8
  631. *Nset, nset=Set-6
  632. 3, 4, 7, 8
  633. *Nset, nset=Set-7, generate
  634. 1, 4, 1
  635. *Uel Property, elset=main_element
  636. {uel_props}
  637. *Solid Section, elset=dummy_mesh, material=dummy_material
  638. *Material, name=dummy_material
  639. *User output variables
  640. {nstatev},
  641. *Elastic
  642. 1.e-20
  643. *Step, name=static, nlgeom=YES, unsymm=YES, inc=200
  644. *Static
  645. 0.01, 1., 1e-05, 0.1
  646. *INCLUDE, file=bcs_uni.inp
  647. *OUTPUT,FIELD,VARIABLE=PRESELECT,FREQ=1
  648. *ELEMENT OUTPUT, elset=dummy_mesh
  649. UVARM
  650. *OUTPUT,HISTORY,VARIABLE=PRESELECT,FREQ=1
  651. *End Step
  652. """
  653. (abq / "uel_cube.inp").write_text(deck)
  655. run_sh = """\
  656. #!/bin/bash
  657. # Run ABAQUS single-element UEL test
  658. # Usage: ./run_uel.sh [bcs_file]
  659. BCS=${1:-bcs_uni.inp}
  660. sed -i "s/INCLUDE, file=bcs_.*/INCLUDE, file=${BCS}/" uel_cube.inp
  661. abaqus job=uel_cube user=../uel.f90 interactive
  662. """
  663. run_path = abq / "run_uel.sh"
  664. run_path.write_text(run_sh)
  665. run_path.chmod(run_path.stat().st_mode | stat.S_IEXEC)
  668. def generate_test_driver(outdir, cfg):
  669. props = build_props(cfg)
  670. nprops = len(props)
  671. nstatev = compute_nstatev(cfg)
  672. test = cfg.get("test", {})
  673. stretch_max = test.get("stretch_max", 1.5)
  674. gamma_max = test.get("gamma_max", 0.6)
  675. nsteps = test.get("nsteps", 400)
  676. dtime = test.get("dtime", 0.01)
  678. code = f"""\
  679. ! Auto-generated test driver for material: {cfg["name"]}
  680. ! Runs uniaxial, biaxial, pure shear, and simple shear tests.
  681. program test_umat
  682. implicit none
  684. integer, parameter :: ntens = 6, ndi = 3, nshr = 3
  685. integer, parameter :: nprops = {nprops}, nstatev = {nstatev}
  686. integer, parameter :: noel = 1, npt = 1
  688. double precision :: stress(ntens), statev(nstatev), ddsdde(ntens, ntens)
  689. double precision :: ddsddt(ntens), drplde(ntens)
  690. double precision :: stran(ntens), dstran(ntens)
  691. double precision :: time(2), predef(1), dpred(1)
  692. double precision :: props(nprops), coords(3), drot(3,3)
  693. double precision :: dfgrd0(3,3), dfgrd1(3,3)
  694. double precision :: sse, spd, scd, rpl, drpldt, dtime
  695. double precision :: temp, dtemp, pnewdt, celent
  696. character(len=8) :: cmname
  697. integer :: layer, kspt, kstep, kinc
  699. integer :: nsteps, i
  700. double precision :: stretch, dstretch, gamma_val, dgamma
  701. double precision, parameter :: zero = 0.0d0, one = 1.0d0
  703. ! --- Initialize all arrays ---
  704. stress = zero; statev = zero; ddsdde = zero
  705. stran = zero; dstran = zero; time = zero
  706. drot = zero; dfgrd0 = zero; dfgrd1 = zero
  707. coords = zero; predef = zero; dpred = zero
  708. temp = zero; dtemp = zero; pnewdt = one; celent = one
  709. rpl = zero; drpldt = zero; ddsddt = zero; drplde = zero
  710. layer = 1; kspt = 1; kinc = 1; cmname = 'UMAT'
  711. dfgrd0(1,1) = one; dfgrd0(2,2) = one; dfgrd0(3,3) = one
  713. ! --- Initialize external database (for RW network models) ---
  714. call uexternaldb(0, 0, time, zero, 0, 0)
  716. ! --- Material properties ---
  717. {fmt_props_fortran(props)}
  719. ! --- Test parameters ---
  720. nsteps = {nsteps}
  721. dtime = {dtime}d0
  722. dstretch = ({stretch_max}d0 - one) / nsteps
  723. dgamma = (2.0d0 * {gamma_max}d0) / nsteps
  725. call execute_command_line('mkdir -p results')
  727. ! ========================== UNIAXIAL ==========================
  728. stress = zero; statev = zero; time = zero
  729. dfgrd1 = zero; dfgrd1(1,1) = one; dfgrd1(2,2) = one; dfgrd1(3,3) = one
  730. stretch = one
  731. open(unit=10, file='results/uniaxial.dat', status='replace')
  732. do i = 1, nsteps
  733. dfgrd1(1,1) = stretch
  734. dfgrd1(2,2) = one / sqrt(stretch)
  735. dfgrd1(3,3) = one / sqrt(stretch)
  736. call umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, &
  737. drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, &
  738. predef, dpred, cmname, ndi, nshr, ntens, nstatev, props, &
  739. nprops, coords, drot, pnewdt, celent, dfgrd0, dfgrd1, &
  740. noel, npt, layer, kspt, kstep, kinc)
  741. write(10, '(3ES20.10)') time(1), stretch, stress(1)
  742. time(1) = time(1) + dtime
  743. stretch = stretch + dstretch
  744. end do
  745. close(10)
  746. write(*, '(A)') 'Uniaxial -> results/uniaxial.dat'
  748. ! ========================== BIAXIAL ==========================
  749. stress = zero; statev = zero; time = zero
  750. dfgrd1 = zero; dfgrd1(1,1) = one; dfgrd1(2,2) = one; dfgrd1(3,3) = one
  751. stretch = one
  752. open(unit=11, file='results/biaxial.dat', status='replace')
  753. do i = 1, nsteps
  754. dfgrd1(1,1) = stretch
  755. dfgrd1(2,2) = stretch
  756. dfgrd1(3,3) = one / (stretch * stretch)
  757. call umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, &
  758. drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, &
  759. predef, dpred, cmname, ndi, nshr, ntens, nstatev, props, &
  760. nprops, coords, drot, pnewdt, celent, dfgrd0, dfgrd1, &
  761. noel, npt, layer, kspt, kstep, kinc)
  762. write(11, '(3ES20.10)') time(1), stretch, stress(1)
  763. time(1) = time(1) + dtime
  764. stretch = stretch + dstretch
  765. end do
  766. close(11)
  767. write(*, '(A)') 'Biaxial -> results/biaxial.dat'
  769. ! ========================== PURE SHEAR ==========================
  770. stress = zero; statev = zero; time = zero
  771. dfgrd1 = zero; dfgrd1(1,1) = one; dfgrd1(2,2) = one; dfgrd1(3,3) = one
  772. gamma_val = -{gamma_max}d0
  773. open(unit=12, file='results/shear.dat', status='replace')
  774. do i = 1, nsteps
  775. dfgrd1(1,2) = gamma_val
  776. dfgrd1(2,1) = gamma_val
  777. call umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, &
  778. drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, &
  779. predef, dpred, cmname, ndi, nshr, ntens, nstatev, props, &
  780. nprops, coords, drot, pnewdt, celent, dfgrd0, dfgrd1, &
  781. noel, npt, layer, kspt, kstep, kinc)
  782. write(12, '(3ES20.10)') time(1), gamma_val, stress(4)
  783. time(1) = time(1) + dtime
  784. gamma_val = gamma_val + dgamma
  785. end do
  786. close(12)
  787. write(*, '(A)') 'Shear -> results/shear.dat'
  789. ! ======================== SIMPLE SHEAR ========================
  790. stress = zero; statev = zero; time = zero
  791. dfgrd1 = zero; dfgrd1(1,1) = one; dfgrd1(2,2) = one; dfgrd1(3,3) = one
  792. gamma_val = -{gamma_max}d0
  793. open(unit=13, file='results/simple_shear.dat', status='replace')
  794. do i = 1, nsteps
  795. dfgrd1(1,2) = gamma_val
  796. call umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, &
  797. drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, &
  798. predef, dpred, cmname, ndi, nshr, ntens, nstatev, props, &
  799. nprops, coords, drot, pnewdt, celent, dfgrd0, dfgrd1, &
  800. noel, npt, layer, kspt, kstep, kinc)
  801. write(13, '(3ES20.10)') time(1), gamma_val, stress(4)
  802. time(1) = time(1) + dtime
  803. gamma_val = gamma_val + dgamma
  804. end do
  805. close(13)
  806. write(*, '(A)') 'Simple sh -> results/simple_shear.dat'
  808. end program test_umat
  810. ! Stub for ABAQUS-provided routine (standalone builds only)
  811. subroutine getoutdir(outdir, lenoutdir)
  812. implicit none
  813. character(len=256), intent(out) :: outdir
  814. integer, intent(out) :: lenoutdir
  815. outdir = '.'
  816. lenoutdir = 1
  817. end subroutine getoutdir
  818. """
  819. (outdir / "test_umat.f90").write_text(code)
  822. def generate_makefile(outdir, elem=None):
  823. uel_all = " test_uel" if elem else ""
  824. mk = f"""\
  825. FC = gfortran
  826. FFLAGS = -O2 -ffree-form
  828. all: test_umat{uel_all}
  830. test_umat: umat.f90 test_umat.f90
  831. \t$(FC) $(FFLAGS) -o $@ $^
  833. run: test_umat
  834. \t./test_umat
  835. """
  836. if elem:
  837. mk += """
  838. test_uel: uel.f90 test_uel.f90
  839. \t$(FC) $(FFLAGS) -o $@ $^
  841. run_uel: test_uel
  842. \t./test_uel
  843. """
  844. mk += """
  845. clean:
  846. \trm -f test_umat test_uel *.mod
  847. \trm -rf results
  848. """
  849. (outdir / "Makefile").write_text(mk)
  852. def generate_abaqus_dir(outdir, cfg):
  853. abq = outdir / "abaqus"
  854. abq.mkdir(exist_ok=True)
  856. props = build_props(cfg)
  857. nprops = len(props)
  858. nstatev = compute_nstatev(cfg)
  860. # --- cube.inp ---
  861. cube = """\
  862. ** Unit cube, single C3D8 element
  863. *Node, nset=all_nodes
  864. 1, 1., 1., 1.
  865. 2, 1., 0., 1.
  866. 3, 1., 1., 0.
  867. 4, 1., 0., 0.
  868. 5, 0., 1., 1.
  869. 6, 0., 0., 1.
  870. 7, 0., 1., 0.
  871. 8, 0., 0., 0.
  872. *Element, type=C3D8, elset=main_element
  873. 1, 5, 6, 8, 7, 1, 2, 4, 3
  874. *Nset, nset=Set-1, generate
  875. 2, 8, 2
  876. *Nset, nset=Set-2, generate
  877. 1, 7, 2
  878. *Nset, nset=Set-3
  879. 1, 2, 5, 6
  880. *Nset, nset=Set-4, generate
  881. 5, 8, 1
  882. *Nset, nset=Set-5
  883. 2, 4, 6, 8
  884. *Nset, nset=Set-6
  885. 3, 4, 7, 8
  886. *Nset, nset=Set-7, generate
  887. 1, 4, 1
  888. *Nset, nset=Set-8
  889. 1,3
  890. *Nset, nset=Set-9
  891. 6,8
  892. *Elset, elset=Surf
  893. 1,
  894. *Surface, type=ELEMENT, name=Surf-1
  895. Surf, S1
  896. *Surface, type=ELEMENT, name=Surf-2
  897. Surf, S2
  898. *Surface, type=ELEMENT, name=Surf-3
  899. Surf, S3
  900. *Surface, type=ELEMENT, name=Surf-4
  901. Surf, S4
  902. *Surface, type=ELEMENT, name=Surf-5
  903. Surf, S5
  904. *Surface, type=ELEMENT, name=Surf-6
  905. Surf, S6
  906. *INCLUDE, file=sec.inp
  907. *Step, name=static, nlgeom=Yes, inc=200
  908. *Static
  909. 0.01, 1., 1e-05, 0.1
  910. *INCLUDE, file=bcs_uni.inp
  911. *OUTPUT,FIELD,VARIABLE=PRESELECT,FREQ=1
  912. *ELEMENT OUTPUT, elset=main_element
  913. SDV
  914. *OUTPUT,HISTORY,VARIABLE=PRESELECT,FREQ=1
  915. *End Step
  916. """
  917. (abq / "cube.inp").write_text(cube)
  919. # --- sec.inp ---
  920. sdv_lines = f"{nstatev},\n1, DET, \"DET\""
  921. if cfg["damage_type"] > 0:
  922. sdv_lines += "\n2, DMG, \"damage\"\n3, MAXSEF, \"max SEF\""
  923. sec = f"""\
  924. *Solid Section, elset=main_element, material=UD
  925. *Material, name=UD
  926. *User Material, constants={nprops}
  927. {fmt_props_abaqus(props)}
  928. *DEPVAR
  929. {sdv_lines}
  930. """
  931. (abq / "sec.inp").write_text(sec)
  933. # --- bcs_uni.inp ---
  934. bcs_uni = """\
  935. *Boundary
  936. Set-3, ZSYMM
  937. *Boundary
  938. Set-4, XSYMM
  939. *Boundary
  940. Set-1, YSYMM
  941. *Boundary, type=displacement
  942. Set-2, 2,2, 0.6
  943. """
  944. (abq / "bcs_uni.inp").write_text(bcs_uni)
  946. # --- bcs_bi.inp ---
  947. bcs_bi = """\
  948. *Boundary
  949. Set-1, YSYMM
  950. *Boundary
  951. Set-3, ZSYMM
  952. *Boundary
  953. Set-4, XSYMM
  954. *Boundary
  955. Set-2, 2,2, 0.6
  956. *Boundary
  957. Set-7, 1,1, 0.6
  958. """
  959. (abq / "bcs_bi.inp").write_text(bcs_bi)
  961. # --- bcs_sh.inp ---
  962. bcs_sh = """\
  963. *Boundary
  964. Set-3, ZSYMM
  965. *Boundary
  966. Set-1, 1,2, 0.0
  967. *Boundary
  968. Set-2, 1,1, 0.6
  969. Set-2, 2,2, 0.0
  970. """
  971. (abq / "bcs_sh.inp").write_text(bcs_sh)
  973. # --- run.sh ---
  974. run_sh = """\
  975. #!/bin/bash
  976. # Run ABAQUS single-element test
  977. # Usage: ./run.sh [bcs_file]
  978. # ./run.sh # uniaxial (default)
  979. # ./run.sh bcs_bi.inp # biaxial
  980. # ./run.sh bcs_sh.inp # shear
  981. BCS=${1:-bcs_uni.inp}
  982. # Update the included BCS file
  983. sed -i "s/INCLUDE, file=bcs_.*/INCLUDE, file=${BCS}/" cube.inp
  984. abaqus job=cube user=../umat.f90 interactive
  985. """
  986. run_path = abq / "run.sh"
  987. run_path.write_text(run_sh)
  988. run_path.chmod(run_path.stat().st_mode | stat.S_IEXEC)
  991. # --- Main ---------------------------------------------------------------------
  993. def list_models():
  994. print("""
  995. Available model types:
  997. ISO_TYPE (isotropic):
  998. 0 None
  999. 1 Neo-Hookean params: C10
  1000. 2 Mooney-Rivlin params: C10, C01
  1001. 3 Ogden (N-term) params: N, mu1, alpha1, ..., muN, alphaN
  1002. 4 Humphrey exponential params: C10, C01
  1004. ANISO_TYPE (anisotropic, per fiber family):
  1005. 0 None
  1006. 1 HGO (dispersed) params: K1, K2, kappa, fiber_x, fiber_y, fiber_z
  1007. 2 Humphrey fiber params: K1, K2, fiber_x, fiber_y, fiber_z
  1008. 3 HGO (AI discrete) params: K1, K2, bdisp, factor, fiber_x, fiber_y, fiber_z
  1009. 4 Humphrey (AI discrete) params: K1, K2, bdisp, factor, fiber_x, fiber_y, fiber_z
  1010. 5 Humphrey + activation params: K1, K2, fiber_x, fiber_y, fiber_z, T0M
  1012. NETWORK_TYPE:
  1013. 0 None
  1014. 1 Affine (RW) params: PHI, N, B_orient, EFI, pdir(3),
  1015. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC
  1016. 2 Non-affine (RW) params: PHI, N, B_orient, EFI, PP,
  1017. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC
  1018. 3 Mixed (AI) params: PHI, N_naff, PP, N_aff, B_orient, EFI, factor, pdir(3),
  1019. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC
  1020. 4 Contractile (AI) params: PHI, N, B_orient, EFI, FRIC, FFMAX, factor, pdir(3),
  1021. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC, KCH(7)
  1022. 5 Affine (AI) params: PHI, N, B_orient, EFI, factor, pdir(3),
  1023. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC
  1024. 6 Non-affine (AI) params: PHI, N, PP, factor,
  1025. L, R0F, mu0, beta, B0, lambda0, R0C, ETAC
  1027. DAMAGE_TYPE:
  1028. 0 None
  1029. 1 Sigmoid params: beta_d, psi_half
  1031. VISCO (per Maxwell branch, max 3):
  1032. params: tau, theta (per branch)
  1034. Example configs: """ + ", ".join(EXAMPLES.keys()))
  1037. def sphere_quadrature(n=60):
  1038. """Quasi-uniform unit-sphere quadrature in the format UEXTERNALDB reads
  1039. ('x y z weight' per line, weights = 1/n). Fibonacci spiral; a generic
  1040. orientation quadrature for the RW network types (1, 2). Replace with a
  1041. higher-order spherical design for production accuracy if needed."""
  1042. import math
  1043. ga = math.pi * (3.0 - math.sqrt(5.0))
  1044. w = 1.0 / n
  1045. lines = []
  1046. for i in range(n):
  1047. z = 1.0 - 2.0 * (i + 0.5) / n
  1048. r = math.sqrt(max(0.0, 1.0 - z * z))
  1049. th = ga * i
  1050. lines.append(f"{r*math.cos(th):.13f} {r*math.sin(th):.13f} {z:.13f} {w:.13f}")
  1051. return "\n".join(lines) + "\n"
  1054. def generate_quadrature(outdir, cfg):
  1055. """Ship a sphere quadrature for the RW network types (1, 2), which read it
  1056. via UEXTERNALDB. AI types (3-6) integrate over an icosahedron and need none."""
  1057. if cfg["network_type"] in (1, 2):
  1058. content = sphere_quadrature(60)
  1059. (outdir / "sphere_int60c.inp").write_text(content)
  1060. (outdir / "abaqus" / "sphere_int60c.inp").write_text(content)
  1063. def generate(cfg):
  1064. name = cfg["name"]
  1065. outdir = SCRIPT_DIR / name
  1066. outdir.mkdir(exist_ok=True)
  1068. elem = element_cfg(cfg)
  1070. generate_umat_f90(outdir)
  1071. generate_aba_param(outdir)
  1072. generate_test_driver(outdir, cfg)
  1073. generate_makefile(outdir, elem)
  1074. generate_abaqus_dir(outdir, cfg)
  1075. generate_quadrature(outdir, cfg)
  1076. if elem:
  1077. generate_uel_f90(outdir, cfg, elem)
  1078. generate_uel_test_driver(outdir, cfg, elem)
  1079. generate_uel_abaqus(outdir, cfg, elem)
  1081. # Save the config for reproducibility
  1082. (outdir / "config.json").write_text(json.dumps(cfg, indent=2) + "\n")
  1084. nprops = len(build_props(cfg))
  1085. nstatev = compute_nstatev(cfg)
  1087. print(f"\nGenerated material: {name}/")
  1088. print(f" umat.f90 Concatenated UMAT source ({nprops} PROPS, {nstatev} STATEV)")
  1089. print(f" test_umat.f90 Standalone test driver")
  1090. print(f" Makefile Build & run: make run")
  1091. print(f" aba_param.inc ABAQUS include")
  1092. print(f" config.json Configuration (for regeneration)")
  1093. print(f" abaqus/ ABAQUS single-element test")
  1094. print(f" cube.inp C3D8 mesh + step definition")
  1095. print(f" sec.inp *User Material card")
  1096. print(f" bcs_uni.inp Uniaxial boundary conditions")
  1097. print(f" bcs_bi.inp Biaxial boundary conditions")
  1098. print(f" bcs_sh.inp Shear boundary conditions")
  1099. print(f" run.sh ABAQUS submission script")
  1100. if elem:
  1101. nvars = elem["nint"] * nstatev
  1102. print(f" uel.f90 User element + material library "
  1103. f"(U3, {elem['nint']}-pt, F-bar={'on' if elem['fbar'] else 'off'}, Variables={nvars})")
  1104. print(f" test_uel.f90 Standalone single-element driver")
  1105. print(f" abaqus/uel_cube.inp + run_uel.sh ABAQUS UEL test")
  1107. print(f"\nStandalone test: cd {name} && make run")
  1108. if elem:
  1109. print(f"UEL test: cd {name} && make run_uel")
  1110. print(f"ABAQUS test: cd {name}/abaqus && ./run.sh")
  1113. def main():
  1114. parser = argparse.ArgumentParser(
  1115. description="Generate a self-contained UMAT material law directory."
  1116. )
  1117. parser.add_argument("config", nargs="?", help="JSON configuration file")
  1118. parser.add_argument("--example", metavar="NAME",
  1119. help="Generate from built-in example: " + ", ".join(EXAMPLES.keys()))
  1120. parser.add_argument("--list", action="store_true", help="List available model types")
  1121. parser.add_argument("--uel", action="store_true",
  1122. help="Also emit a user element (uel.f90 + test_uel.f90 + ABAQUS "
  1123. "UEL deck) driven by this material; configs may instead "
  1124. "carry an explicit \"element\" block")
  1125. args = parser.parse_args()
  1127. if args.list:
  1128. list_models()
  1129. return
  1131. if args.example:
  1132. if args.example not in EXAMPLES:
  1133. print(f"Unknown example: {args.example}")
  1134. print(f"Available: {', '.join(EXAMPLES.keys())}")
  1135. sys.exit(1)
  1136. cfg = dict(EXAMPLES[args.example])
  1137. if args.uel:
  1138. cfg.setdefault("element", {})
  1139. generate(cfg)
  1140. return
  1142. if args.config:
  1143. with open(args.config) as f:
  1144. cfg = json.load(f)
  1145. if args.uel:
  1146. cfg.setdefault("element", {})
  1147. generate(cfg)
  1148. return
  1150. parser.print_help()
  1153. if __name__ == "__main__":
  1154. main()

generate.py at commit e597e91, no license · at the source

Overview

Authors: Luís Pacheco1,2, Marco Parente1,2, João Ferreira1
  1. Department of Mechanical Engineering, Faculty of Engineering, University of Porto, R. Dr. Roberto Frias, 4200-465 Porto, Portugal
  2. Institute of Mechanical Engineering and Industrial Management, R. Dr. Roberto Frias, 4200-465 Porto, Portugal
Institutions: Universidade do Porto (Portugal)
Journal: Biomechanics and modeling in mechanobiology, volume 25, issue 3, article 52
Dates: received 9 December 2025; accepted 23 March 2026; published online 3 June 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1007/s10237-026-02067-5 · PMID 42234214 · PMCID PMC13234050 · OpenAlex W4417483186
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.21203/rs.3.rs-8315851/v1
Status: code verified
Categories: computational modeling (no new data) (modality), none (in silico) (organism), computational (subfield)
Methods: Statistics
Keywords: Actin networks, Continuum mechanics, Finite elements, Sparse polynomial chaos expansions, Uncertainty quantification
MeSH: Actins*, Cross-Linking Reagents*, Animals, Computer Simulation, Elasticity, Finite Element Analysis, Models, Biological, Monte Carlo Method, Stochastic Processes, Stress, Mechanical (* major topic)
Topic: Cellular Mechanics and Interactions (Cell Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Universidade do Porto
Citations: not cited yet (Europe PMC); 65 references in the paper

Abstract

Filamentous actin (F-actin) constitutes the primary contributor to cell elasticity and structural integrity, forming dynamic, crosslinked networks in the actin cortex. Existing mechanical models for F-actin and crosslinked filament networks successfully describe filament- and network-level behavior, but are often limited in accounting for biological dynamic processes and inherent material uncertainty and variability. We develop a stochastic modeling framework that integrates Polynomial Chaos Expansion (PCE) surrogates using the Finite Element Method (FEM). These surrogates replace filament-scale equations for compliant crosslinked F-actin networks, efficiently enabling uncertainty quantification and sensitivity analysis of key material parameters. The first and second statistical moments from the PCE are incorporated into a micro-sphere network model and implemented via a user-defined material subroutine. Validation was performed against 10 000 Monte Carlo simulations (MCS) for each of four FEM test cases: three simple deformation modes applied to a unit length cubic element, and a thin gel layer under shear mimicking a parallel plate rheology setup. In every test, the surrogate predicts the expected value of relevant stress quantities at maximum deformation with under 1% relative error versus the MCS reference. Moreover, the surrogate captures the network’s variability as measured by second-order moments, demonstrating its ability to deliver rapid, statistically faithful predictions of both mean response and standard deviation in simple element tests and experimentally relevant rheology geometries. The proposed methodology provides a scalable route for incorporating intrinsic material variability into F-actin mechanical modeling, with implications for studying cell motility, division, and pathologies related to cytoskeletal remodeling.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

jpsferreira/UMAT-ABAQUS_library

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e597e91cfd727b6d27f52389c7b99b90029d9c56, 12 June 2026
Languages: Fortran (585), Shell (63), Python (48), NEURON (15)
Size: 2,669 files, 711 scripts
Software Heritage: not archived
Found in: the text, “Methods”
Holds: README, environment (pyproject.toml, uv.lock), tests, documentation
Not found: license file, CITATION.cff, continuous integration
Tools: NumPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
712 files
At the source: github.com/jpsferreira/UMAT-ABAQUS_library

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 711 scripts, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Springer Science+Business Media

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 5 keywords, 10 MeSH terms, 1 funder, 61 references.

Cite

This paper

Pacheco, L., Parente, M., & Ferreira, J. (2026). Sparse polynomial surrogates for F-actin networks with compliant crosslinkers. Biomechanics and modeling in mechanobiology, 25(3), 52. https://doi.org/10.1007/s10237-026-02067-5

BibTeX

@article{pacheco2026sparse,
author = {Pacheco, Luís and Parente, Marco and Ferreira, João},
title = {{Sparse polynomial surrogates for F-actin networks with compliant crosslinkers}},
journal = {Biomechanics and modeling in mechanobiology},
year = {2026},
month = jun,
volume = {25},
number = {3},
pages = {52},
publisher = {Springer Science+Business Media},
issn = {1617-7959},
doi = {10.1007/s10237-026-02067-5},
url = {https://doi.org/10.1007/s10237-026-02067-5},
pmid = {42234214},
pmcid = {PMC13234050}
}

RIS

TY - JOUR
AU - Pacheco, Luís
AU - Parente, Marco
AU - Ferreira, João
TI - Sparse polynomial surrogates for F-actin networks with compliant crosslinkers
T2 - Biomechanics and modeling in mechanobiology
J2 - Biomech Model Mechanobiol
PY - 2026
DA - 2026/06/03
VL - 25
IS - 3
SP - 52
SN - 1617-7959
PB - Springer Science+Business Media
DO - 10.1007/s10237-026-02067-5
UR - https://doi.org/10.1007/s10237-026-02067-5
LA - en
ER -

CSL-JSON

{
"id": "10.1007/s10237-026-02067-5",
"type": "article-journal",
"title": "Sparse polynomial surrogates for F-actin networks with compliant crosslinkers",
"container-title": "Biomechanics and modeling in mechanobiology",
"author": [
{
"family": "Pacheco",
"given": "Luís"
},
{
"family": "Parente",
"given": "Marco"
},
{
"family": "Ferreira",
"given": "João"
}
],
"container-title-short": "Biomech Model Mechanobiol",
"volume": "25",
"issue": "3",
"page": "52",
"DOI": "10.1007/s10237-026-02067-5",
"PMID": "42234214",
"PMCID": "PMC13234050",
"ISSN": "1617-7959",
"publisher": "Springer Science+Business Media",
"URL": "https://doi.org/10.1007/s10237-026-02067-5",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
3
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76045-x [code]
A manufacturability-informed topology framework for AI-guided design of fibrous network materials.
Journal: Nature communications
In common: NumPy, 2 references
[2] doi:10.1002/2211-5463.70302
KDAC6 alters cell morphology and motility through positive and negative modulation of F-actin distribution.
Journal: FEBS open bio
In common: 2 references
[3] doi: [code]
Going deeper with morphologically detailed neural networks by simulation-based gradient propagation
Journal: Frontiers in computational neuroscience
In common: NumPy, none (in silico), computational, computational modeling (no new data)
[4] doi:10.1038/s41526-026-00644-7 [code]
A computational model of altered neuronal activity in altered gravity.
Journal: NPJ microgravity
In common: NumPy, none (in silico), computational, computational modeling (no new data)
[5] doi:10.1371/journal.pcbi.1014458 [code]
Neuronal excitability and parameter variability in the Hodgkin-Huxley model.
Journal: PLoS computational biology
In common: NumPy, none (in silico), computational, computational modeling (no new data)
[6] doi:10.1038/s41598-026-51212-8 [code]
Uncertainty aware machine learning for bridging simulation and experiment in high throughput materials characterization.
Journal: Scientific reports
In common: NumPy, none (in silico), computational, computational modeling (no new data)
[7] doi:10.1186/s12860-026-00584-w
A role for PaxB in regulating blebbing: experimental insights and theoretical perspectives from Dictyostelium discoideum.
Journal: BMC molecular and cell biology
In common: none (in silico), 1 reference
[8] doi:10.1371/journal.pcbi.1014617 [code]
An in silico framework for dissecting the mechanistic origins of in vivo recorded neuronal activity.
Journal: PLoS computational biology
In common: none (in silico), computational, computational modeling (no new data)
[9] doi:10.1007/s11571-026-10522-3 [code]
Acetylcholine enhances deviance detection in Hodgkin-Huxley neuronal networks.
Journal: Cognitive neurodynamics
In common: none (in silico), computational, computational modeling (no new data)
[10] doi:10.3389/fncom.2026.1799705
Structural synaptogenesis superior to functional modulation in a pruning-based recurrent network model of OCD.
Journal: Frontiers in computational neuroscience
In common: none (in silico), computational, computational modeling (no new data)

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.