From 533efa1dec2666f6077a0f1ae712b60a29315dd8 Mon Sep 17 00:00:00 2001 From: Sjard Mathis Rosenbusch <69848361+srosenbu@users.noreply.github.com> Date: Wed, 14 Aug 2024 10:08:11 +0200 Subject: [PATCH] Hardening (#20) * Start on JHR model * Parameter for tensile collapse of yield surface * Functional programming * exponential damage law * gradient model in terms of total nonlocal strain * fix some errors * fix another error * WIP: nonlocal balance in init conf * fix key error * initial_config in implicit * Fix error in initial connfig * error in mesh update * fix condition for init config * Fix * local sound speed * stuff * Overnonlocal model * Fix negative damage * e_0 in damage * Include hardening * Fix nan error * logic error in damage * Fix lambda calculation * Add Cargo lock file to version control --------- Co-authored-by: srosenbu --- Cargo.lock | 496 ++++++++++++++++++++++++++++++++++++++++++++ python/comfe/cdm.py | 66 ++++-- src/gradient_jh2.rs | 120 +++++++---- src/interfaces.rs | 3 + src/jh2.rs | 68 ++++-- src/jhr.rs | 278 +++++++++++++++++++++++++ src/lib.rs | 21 +- src/smallstrain.rs | 140 +++++++++---- test/test_jh2_2d.py | 11 +- 9 files changed, 1087 insertions(+), 116 deletions(-) create mode 100644 Cargo.lock create mode 100644 src/jhr.rs diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..7065a30 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,496 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bytemuck" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "17febce684fd15d89027105661fec94afb475cb995fbc59d2865198446ba2eea" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "comfe" +version = "0.1.0" +dependencies = [ + "nalgebra", + "numpy", + "pyo3", + "simba", + "strum", + "strum_macros", +] + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "indoc" +version = "1.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bfa799dd5ed20a7e349f3b4639aa80d74549c81716d9ec4f994c9b5815598306" + +[[package]] +name = "libc" +version = "0.2.147" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4668fb0ea861c1df094127ac5f1da3409a82116a4ba74fca2e58ef927159bb3" + +[[package]] +name = "lock_api" +version = "0.4.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c1cc9717a20b1bb222f333e6a92fd32f7d8a18ddc5a3191a11af45dcbf4dcd16" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "090126dc04f95dc0d1c1c91f61bdd474b3930ca064c1edc8a849da2c6cbe1e77" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "memoffset" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "nalgebra" +version = "0.32.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "307ed9b18cc2423f29e83f84fd23a8e73628727990181f18641a8b5dc2ab1caa" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91761aed67d03ad966ef783ae962ef9bbaca728d2dd7ceb7939ec110fffad998" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "ndarray" +version = "0.15.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb12d4e967ec485a5f71c6311fe28158e9d6f4bc4a447b474184d0f91a8fa32" +dependencies = [ + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "rawpointer", +] + +[[package]] +name = "num-complex" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ba157ca0885411de85d6ca030ba7e2a83a28636056c7c699b07c8b6f7383214" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f30b0abd723be7e2ffca1272140fac1a2f084c77ec3e123c192b66af1ee9e6c2" +dependencies = [ + "autocfg", +] + +[[package]] +name = "numpy" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "437213adf41bbccf4aeae535fbfcdad0f6fed241e1ae182ebe97fa1f3ce19389" +dependencies = [ + "libc", + "nalgebra", + "ndarray", + "num-complex", + "num-integer", + "num-traits", + "pyo3", + "rustc-hash", +] + +[[package]] +name = "once_cell" +version = "1.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" + +[[package]] +name = "parking_lot" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3742b2c103b9f06bc9fff0a37ff4912935851bee6d36f3c02bcc755bcfec228f" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93f00c865fe7cabf650081affecd3871070f26767e7b2070a3ffae14c654b447" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "paste" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c" + +[[package]] +name = "proc-macro2" +version = "1.0.66" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "18fb31db3f9bddb2ea821cde30a9f70117e3f119938b5ee630b7403aa6e2ead9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e681a6cfdc4adcc93b4d3cf993749a4552018ee0a9b65fc0ccfad74352c72a38" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "parking_lot", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "076c73d0bc438f7a4ef6fdd0c3bb4732149136abd952b110ac93e4edb13a6ba5" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e53cee42e77ebe256066ba8aa77eff722b3bb91f3419177cf4cd0f304d3284d9" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dfeb4c99597e136528c6dd7d5e3de5434d1ceaf487436a3f03b2d56b6fc9efd1" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.19.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "947dc12175c254889edc0c02e399476c2f652b4b9ebd123aa655c224de259536" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "quote" +version = "1.0.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "redox_syscall" +version = "0.3.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "567664f262709473930a4bf9e51bf2ebf3348f2e748ccc50dea20646858f8f29" +dependencies = [ + "bitflags", +] + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustversion" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ffc183a10b4478d04cbbbfc96d0873219d962dd5accaff2ffbd4ceb7df837f4" + +[[package]] +name = "safe_arch" +version = "0.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "simba" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "smallvec" +version = "1.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62bb4feee49fdd9f707ef802e22365a35de4b7b299de4763d44bfea899442ff9" + +[[package]] +name = "strum" +version = "0.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "290d54ea6f91c969195bdbcd7442c8c2a2ba87da8bf60a7ee86a235d4bc1e125" + +[[package]] +name = "strum_macros" +version = "0.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ad8d03b598d3d0fff69bf533ee3ef19b8eeb342729596df84bcc7e1f96ec4059" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "rustversion", + "syn 2.0.29", +] + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.29" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c324c494eba9d92503e6f1ef2e6df781e78f6a7705a0202d9801b198807d518a" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9d0e916b1148c8e263850e1ebcbd046f333e0683c724876bb0da63ea4373dc8a" + +[[package]] +name = "typenum" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba" + +[[package]] +name = "unicode-ident" +version = "1.0.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "301abaae475aa91687eb82514b328ab47a211a533026cb25fc3e519b86adfc3c" + +[[package]] +name = "unindent" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1766d682d402817b5ac4490b3c3002d91dfa0d22812f341609f97b08757359c" + +[[package]] +name = "wide" +version = "0.7.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa469ffa65ef7e0ba0f164183697b89b854253fd31aeb92358b7b6155177d62f" +dependencies = [ + "bytemuck", + "safe_arch", +] + +[[package]] +name = "windows-targets" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a2fa6e2155d7247be68c096456083145c183cbbbc2764150dda45a87197940c" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b38e32f0abccf9987a4e3079dfb67dcd799fb61361e53e2882c3cbaf0d905d8" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dc35310971f3b2dbbf3f0690a219f40e2d9afcf64f9ab7cc1be722937c26b4bc" + +[[package]] +name = "windows_i686_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a75915e7def60c94dcef72200b9a8e58e5091744960da64ec734a6c6e9b3743e" + +[[package]] +name = "windows_i686_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f55c233f70c4b27f66c523580f78f1004e8b5a8b659e05a4eb49d4166cca406" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53d40abd2583d23e4718fddf1ebec84dbff8381c07cae67ff7768bbf19c6718e" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b7b52767868a23d5bab768e390dc5f5c55825b6d30b86c844ff2dc7414044cc" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.48.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed94fce61571a4006852b7389a063ab983c02eb1bb37b47f8272ce92d06d9538" diff --git a/python/comfe/cdm.py b/python/comfe/cdm.py index 4ed33ae..e2e55b8 100644 --- a/python/comfe/cdm.py +++ b/python/comfe/cdm.py @@ -361,7 +361,7 @@ class CDMNonlocal(NonlocalInterface): M: df.fem.Function # fields: dict[str, df.fem.Function] # q_fields: dict[str, df.fem.Function | ufl.core.expr.Expr] - rate_evaluator: QuadratureEvaluator + strain_evaluator: QuadratureEvaluator parameters: dict[str, float] form: df.fem.FormMetaClass t: float @@ -379,17 +379,17 @@ def __init__( parameters: dict[str, float], quadrature_rule: QuadratureRule, q_fields_local: dict[str, df.fem.Function], - q_field_nonlocal_rate: df.fem.Function, + q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, + displacements: df.fem.Function | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} if Q_nonlocal_rate[-4:] != "rate": raise ValueError("Q_nonlocal_rate must be some rate for CDM, you provided " + Q_nonlocal_rate) - q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate - Q_nonlocal = Q_nonlocal_rate[:-5] + q_fields[Q_nonlocal] = q_field_nonlocal fields = { Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), @@ -408,6 +408,9 @@ def __init__( else: g = 1.0 + if displacements is not None: + fields["u"] = displacements + f_int_ufl = ( g * parameters["l"] ** 2 * ufl.inner(ufl.grad(fields[Q_nonlocal]), ufl.grad(test_function)) + fields[Q_nonlocal] * test_function @@ -419,7 +422,7 @@ def __init__( f_form = df.fem.form(f_ufl) - rate_evaluator = QuadratureEvaluator(fields[Q_nonlocal_rate], function_space.mesh, quadrature_rule) + strain_evaluator = QuadratureEvaluator(fields[Q_nonlocal], function_space.mesh, quadrature_rule) super().__init__( Q_local=Q_local, Q_nonlocal=Q_nonlocal, @@ -431,13 +434,17 @@ def __init__( q_fields=q_fields, fields=fields, form=f_form, - rate_evaluator=rate_evaluator, + strain_evaluator=strain_evaluator, ) def step(self, h: float) -> None: with self.fields["nonlocal_force"].vector.localForm() as f_local: f_local.set(0.0) + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, -self.fields["u"].x.array, mode="add") + # self.fields["u"].function_space.mesh.geometry.x[:] -= self.fields["u"].x.array + df.fem.petsc.assemble_vector(self.fields["nonlocal_force"].vector, self.form) self.fields["nonlocal_force"].x.scatter_reverse(ScatterMode.add) @@ -454,8 +461,12 @@ def step(self, h: float) -> None: self.fields[self.Q_nonlocal].vector.array[:] += h * self.fields[self.Q_nonlocal_rate].vector.array self.fields[self.Q_nonlocal].x.scatter_forward() - self.rate_evaluator(self.q_fields[self.Q_nonlocal_rate]) - self.q_fields[self.Q_nonlocal_rate].x.scatter_forward() + self.strain_evaluator(self.q_fields[self.Q_nonlocal]) + self.q_fields[self.Q_nonlocal].x.scatter_forward() + + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, self.fields["u"].x.array, mode="add") + # self.fields["u"].function_space.mesh.geometry.x[:] += self.fields["u"].x.array self.t += h @@ -464,7 +475,7 @@ class ImplicitNonlocal(NonlocalInterface): # M: df.fem.Function # fields: dict[str, df.fem.Function] # q_fields: dict[str, df.fem.Function | ufl.core.expr.Expr] - rate_evaluator: QuadratureEvaluator + strain_evaluator: QuadratureEvaluator parameters: dict[str, float] # form: df.fem.FormMetaClass problem: df.fem.petsc.LinearProblem @@ -483,20 +494,22 @@ def __init__( parameters: dict[str, float], quadrature_rule: QuadratureRule, q_fields_local: dict[str, df.fem.Function], - q_field_nonlocal_rate: df.fem.Function, + q_field_nonlocal: df.fem.Function, Q_local_damage: str | None = None, + displacements: df.fem.Function | None = None, ): q_fields = {Q_local: q_fields_local[Q_local]} if Q_nonlocal_rate[-4:] != "rate": raise ValueError("Q_nonlocal_rate must be some rate for CDM, you provided " + Q_nonlocal_rate) - q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate + # q_fields[Q_nonlocal_rate] = q_field_nonlocal_rate Q_nonlocal = Q_nonlocal_rate[:-5] + q_fields[Q_nonlocal] = q_field_nonlocal fields = { - Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), + # Q_nonlocal_rate: df.fem.Function(function_space, name=Q_nonlocal_rate), Q_nonlocal: df.fem.Function(function_space, name=Q_nonlocal), "nonlocal_force": df.fem.Function(function_space, name="nonlocal_force"), } @@ -510,6 +523,9 @@ def __init__( else: g = 1.0 + if displacements is not None: + fields["u"] = displacements + test_function = ufl.TestFunction(function_space) trial_function = ufl.TrialFunction(function_space) @@ -533,7 +549,7 @@ def __init__( # f_form = df.fem.form(f_ufl) - rate_evaluator = QuadratureEvaluator(fields[Q_nonlocal_rate], function_space.mesh, quadrature_rule) + strain_evaluator = QuadratureEvaluator(fields[Q_nonlocal], function_space.mesh, quadrature_rule) super().__init__( Q_local=Q_local, Q_nonlocal=Q_nonlocal, @@ -545,17 +561,22 @@ def __init__( q_fields=q_fields, fields=fields, # form=f_form, - rate_evaluator=rate_evaluator, + strain_evaluator=strain_evaluator, problem=problem, ) def step(self, h: float) -> None: - old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() + # old_nonlocal = self.fields[self.Q_nonlocal].vector.array.copy() + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, -self.fields["u"].x.array, mode="add") + self.problem.solve() - self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( - self.fields[self.Q_nonlocal].vector.array - old_nonlocal - ) / h - self.rate_evaluator(self.q_fields[self.Q_nonlocal_rate]) + # self.q_fields[self.Q_nonlocal_rate].vector.array[:] = ( + # self.fields[self.Q_nonlocal].vector.array - old_nonlocal + # ) / h + self.strain_evaluator(self.q_fields[self.Q_nonlocal]) + if "u" in self.fields: + set_mesh_coordinates(self.fields["u"].function_space.mesh, self.fields["u"].x.array, mode="add") self.t += h @@ -647,6 +668,7 @@ def __init__( mass_mechanics: df.fem.Function | None = None, mass_nonlocal: df.fem.Function | None = None, calculate_bulk_viscosity: bool = False, + nonlocal_initial_config: bool = True, ) -> None: mass_mechanics = ( mass_mechanics @@ -655,7 +677,7 @@ def __init__( ) mass_nonlocal = ( mass_nonlocal - if mass_nonlocal is not None or nonlocal_parameters["solver"] == "implicit" + if mass_nonlocal is not None or nonlocal_solver == ImplicitNonlocal else diagonal_mass(nonlocal_space, nonlocal_parameters["zeta"], invert=True) ) mechanics_solver = mechanics_solver( @@ -670,6 +692,7 @@ def __init__( calculate_bulk_viscosity=calculate_bulk_viscosity, ) + displacements = mechanics_solver.fields["u"] if nonlocal_initial_config else None nonlocal_solver = nonlocal_solver( Q_local, Q_nonlocal_rate, @@ -679,8 +702,9 @@ def __init__( nonlocal_parameters, quadrature_rule, mechanics_solver.q_fields, - mechanics_solver.model.input[Q_nonlocal_rate], + mechanics_solver.model.input[Q_nonlocal], Q_local_damage=Q_local_damage, + displacements=displacements, ) # add all fields from the solver to this class for easier postprocessing diff --git a/src/gradient_jh2.rs b/src/gradient_jh2.rs index d61e8f6..9f3bae6 100644 --- a/src/gradient_jh2.rs +++ b/src/gradient_jh2.rs @@ -7,13 +7,19 @@ use crate::jh2::JH2ConstParameters; use std::collections::HashMap; +#[derive(Debug)] +pub struct NonlocalParameters { + m_overnonlocal: f64, +} + #[derive(Debug)] pub struct GradientJH23D { parameters: JH2ConstParameters, + nonlocal_parameters: NonlocalParameters, } impl ConstitutiveModel for GradientJH23D { - fn new(parameters: &HashMap) -> Option{ + fn new(parameters: &HashMap) -> Option { Some(Self { parameters: JH2ConstParameters { RHO: *parameters.get("RHO").unwrap(), @@ -35,6 +41,15 @@ impl ConstitutiveModel for GradientJH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), + E_F: *parameters.get("E_F").unwrap_or(&0.0), + E_0: *parameters.get("E_0").unwrap_or(&0.0), + LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), + HARDENING: *parameters.get("HARDENING").unwrap_or(&0.0), + //M_OVERNONLOCAL: *parameters.get("M_OVERNONLOCAL").unwrap_or(&0.0), + //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), + }, + nonlocal_parameters: NonlocalParameters { + m_overnonlocal: *parameters.get("M_OVERNONLOCAL").unwrap_or(&1.0), }, }) } @@ -47,15 +62,21 @@ impl ConstitutiveModel for GradientJH23D { let d_eps = mandel_rate_from_velocity_gradient(&velocity_gradient); let (mut d_eps_vol, d_eps_dev) = mandel_decomposition(&d_eps); d_eps_vol *= -1.0; - + let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); - let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); + //let del_lambda_nonlocal = del_t * input.get_scalar(Q::EqNonlocalPlasticStrainRate, ip).max(0.0); + let nonlocal_plastic_strain = input.get_scalar(Q::EqNonlocalPlasticStrain, ip); + let history_maximum_0 = input.get_scalar(Q::HistoryMaximum, ip); + //let del_lambda_nonlocal = input.get_scalar(Q::EqNonlocalPlasticStrain, ip)- input.get_scalar(Q::HistoryMaximum, ip); + let history_maximum_1 = history_maximum_0.max(nonlocal_plastic_strain); + let del_lambda_nonlocal = history_maximum_1 - history_maximum_0; + output.set_scalar(Q::HistoryMaximum, ip, history_maximum_1); + let mut del_lambda = 0.0; - let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); - + let s_tr = s_0 + 2. * self.parameters.SHEAR_MODULUS * d_eps_dev * del_t; let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); @@ -64,17 +85,37 @@ impl ConstitutiveModel for GradientJH23D { let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - - - let e_p_f = (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)) - .max(self.parameters.EFMIN); + + let e_p_f = + (self.parameters.D1 * (p_s + t_s).powf(self.parameters.D2)).max(self.parameters.EFMIN); let damage_0 = input.get_scalar(Q::Damage, ip); - let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); + let mut damage_1 = 0.0; + if self.parameters.E_F > 0.0 { + let m = self.nonlocal_parameters.m_overnonlocal; + let kappa = (m * history_maximum_1 + + (1.0 - m) * input.get_scalar(Q::EqPlasticStrain, ip)) + .max(0.0); + damage_1 = (1. - ((self.parameters.E_0 - kappa) / self.parameters.E_F).exp()).max(0.0); + } else { + damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); + } + //let damage_1 = (damage_0 + del_lambda_nonlocal / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); + let hardening = self.parameters.HARDENING; + let hardening_slope = hardening / ((1.0 - hardening) * self.parameters.E_0); + let hardening_factor = { + if hardening_slope.is_nan() { + 1.0 + } else { + hardening_slope * input.get_scalar(Q::EqPlasticStrain, ip) + 1.0 + } + }; + //let t_s_factor = (1.-damage_1).powf(self.parameters.REDUCE_T); + let yield_factor = 1.0 - hardening; let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + yield_factor * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -83,9 +124,9 @@ impl ConstitutiveModel for GradientJH23D { } let yield_surface = rate_factor * { if damage_0 == 0.0 { - fracture_surface + fracture_surface * hardening_factor } else { - fracture_surface * (1. - damage_1) + damage_1 * residual_surface + fracture_surface * hardening_factor * (1. - damage_1) + damage_1 * residual_surface } }; if s_tr_eq > yield_surface { @@ -94,8 +135,12 @@ impl ConstitutiveModel for GradientJH23D { } else { alpha = 1.0; } - - output.set_scalar(Q::EqPlasticStrain, ip, input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda); + + output.set_scalar( + Q::EqPlasticStrain, + ip, + input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda, + ); // /*********************************************************************** // * UPDATE DENSITY @@ -109,19 +154,22 @@ impl ConstitutiveModel for GradientJH23D { output.set_scalar(Q::Density, ip, density_1); let mu = density_1 / self.parameters.RHO - 1.; - + let mut local_bulk_modulus = self.parameters.K1; let p_1 = { if mu > 0.0 { + local_bulk_modulus += + 2.0 * self.parameters.K2 * mu + 3.0 * self.parameters.K3 * mu.powi(2); self.parameters.K1 * mu + self.parameters.K2 * mu.powi(2) + self.parameters.K3 * mu.powi(3) + input.get_scalar(Q::BulkingPressure, ip) } else { - let p_trial = self.parameters.K1*mu; + let p_trial = self.parameters.K1 * mu; let p_damaged = -self.parameters.T * (1. - damage_1); if p_trial > p_damaged { p_trial } else { + local_bulk_modulus = 0.0; p_damaged } } @@ -145,10 +193,12 @@ impl ConstitutiveModel for GradientJH23D { // Calculate bulk viscosity let l = input.get_scalar(Q::CellDiameter, ip); let tr_deps = d_eps_vol * 3.; - let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let bulk_modulus = local_bulk_modulus * self.parameters.LOCAL_SOUND_SPEED + + (1.0 - self.parameters.LOCAL_SOUND_SPEED) * self.parameters.K1; + let c = (bulk_modulus / self.parameters.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { - l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) + l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) } else { 0.0 } @@ -159,7 +209,7 @@ impl ConstitutiveModel for GradientJH23D { // * Combine deviatoric and volumetric stresses // **********************************************************************/ let s_1 = s_tr * alpha; - let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1+q_1); + let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1 + q_1); output.set_vector(Q::MandelStress, ip, sigma_1); // *********************************************************************** @@ -178,39 +228,39 @@ impl ConstitutiveModel for GradientJH23D { if output.is_some(Q::Pressure) { output.set_scalar(Q::Pressure, ip, p_1); } - + let elastic_rate = - - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; - + -(1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + let density_mid = 0.5 * (density_0 + density_1); if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); let deviatoric_rate = d_eps_dev - elastic_rate; let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); - let e_1 = e_0 + del_t/density_mid * (s_mid.dot(&deviatoric_rate)); + let e_1 = e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); } if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); - let p_mid = - 0.5 * (p_0 + p_1); + let p_mid = -0.5 * (p_0 + p_1); let deviatoric_rate = elastic_rate; let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); - let e_1 = e_0 + del_t/density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + let e_1 = + e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); output.set_scalar(Q::InternalElasticEnergy, ip, e_1); } if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { - let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let q_mid = -0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); - let e_1 = e_0 + del_t/density_mid * 3.* d_eps_vol * q_mid; + let e_1 = e_0 + del_t / density_mid * 3. * d_eps_vol * q_mid; output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); } if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { let e_0 = input.get_scalar(Q::InternalEnergy, ip); let sigma_mid = 0.5 * (sigma_0 + sigma_1); - let e_1 = e_0 + del_t/density_mid * sigma_mid.dot(&d_eps); + let e_1 = e_0 + del_t / density_mid * sigma_mid.dot(&d_eps); output.set_scalar(Q::InternalEnergy, ip, e_1); } - } /// Returns the physical quantities that are required as input for the @@ -218,7 +268,8 @@ impl ConstitutiveModel for GradientJH23D { fn define_input(&self) -> HashMap { HashMap::from([ (Q::VelocityGradient, QDim::SquareTensor(3)), - (Q::EqNonlocalPlasticStrainRate, QDim::Scalar), + //(Q::EqNonlocalPlasticStrainRate, QDim::Scalar), + (Q::EqNonlocalPlasticStrain, QDim::Scalar), (Q::CellDiameter, QDim::Scalar), ]) } @@ -233,17 +284,16 @@ impl ConstitutiveModel for GradientJH23D { (Q::Density, QDim::Scalar), (Q::EqPlasticStrain, QDim::Scalar), (Q::BulkViscosity, QDim::Scalar), + (Q::HistoryMaximum, QDim::Scalar), ]) } /// Returns the physical quantities that are needed as output, but are not /// necessarily needed in oredr to calculate the constitutive model. An example is - /// the consistent tangent which is not needed for the calculation of the stresses + /// the consistent tangent which is not needed for the calculation of the stresses /// and is therefore purely an output quantity. fn define_output(&self) -> HashMap { - HashMap::from([ - (Q::MandelStress, QDim::Vector(6)), - ]) + HashMap::from([(Q::MandelStress, QDim::Vector(6))]) } /// Returns the physical quantities that are optional output of the constitutive diff --git a/src/interfaces.rs b/src/interfaces.rs index c742c1c..7abf451 100644 --- a/src/interfaces.rs +++ b/src/interfaces.rs @@ -63,6 +63,8 @@ pub enum Q { BulkViscosity, #[strum(serialize = "CellDiameter", serialize = "cell_diameter")] CellDiameter, + #[strum(serialize = "HistoryMaximum", serialize = "history_maximum")] + HistoryMaximum, #[strum(serialize = "_LAST", serialize = "_last")] _LAST, } @@ -234,6 +236,7 @@ impl Q { Q::InternalPlasticEnergyRate => QDim::Scalar, Q::BulkViscosity => QDim::Scalar, Q::CellDiameter => QDim::Scalar, + Q::HistoryMaximum => QDim::Scalar, Q::_LAST => QDim::Scalar, } } diff --git a/src/jh2.rs b/src/jh2.rs index f0498e0..6825ce8 100644 --- a/src/jh2.rs +++ b/src/jh2.rs @@ -3,6 +3,7 @@ use crate::stress_strain::{ mandel_decomposition, mandel_rate_from_velocity_gradient, MANDEL_IDENTITY, }; +use core::panic; use std::collections::HashMap; #[derive(Debug)] pub struct JH2ConstParameters { @@ -25,6 +26,11 @@ pub struct JH2ConstParameters { pub BETA: f64, pub EFMIN: f64, pub DMAX: f64, + pub E_F: f64, + pub E_0: f64, + pub LOCAL_SOUND_SPEED: f64, + pub HARDENING: f64, + //pub REDUCE_T: f64, } #[derive(Debug)] pub struct JH23D { @@ -54,6 +60,11 @@ impl ConstitutiveModel for JH23D { BETA: *parameters.get("BETA").unwrap(), EFMIN: *parameters.get("EFMIN").unwrap(), DMAX: *parameters.get("DMAX").unwrap_or(&1.0), + E_F: *parameters.get("E_F").unwrap_or(&0.0), + E_0: *parameters.get("E_0").unwrap_or(&0.0), + LOCAL_SOUND_SPEED: *parameters.get("LOCAL_SOUND_SPEED").unwrap_or(&0.0), + HARDENING: *parameters.get("HARDENING").unwrap_or(&0.0), + //REDUCE_T: *parameters.get("REDUCE_T").unwrap_or(&0.0), }, }) } @@ -74,20 +85,34 @@ impl ConstitutiveModel for JH23D { let damage_0 = input.get_scalar(Q::Damage, ip); let mut damage_1 = damage_0; + let lambda_old = -self.parameters.E_F * (1. - damage_0).ln() + self.parameters.E_0; + let hardening = self.parameters.HARDENING; + + let hardening_slope = hardening / ((1.0 - hardening) * self.parameters.E_0); + let hardening_factor = { + if hardening_slope.is_nan() { + 1.0 + } else { + hardening_slope * lambda_old + 1.0 + } + }; + let (p_0, s_0) = mandel_decomposition(&sigma_0); let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); let s_tr = s_0 + 2. * self.parameters.SHEAR_MODULUS * d_eps_dev * del_t; let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); - let mut alpha; + let mut alpha: f64; let p_s = p_0 / self.parameters.PHEL; let t_s = self.parameters.T / self.parameters.PHEL; let mut rate_factor = 1.; - let fracture_surface = - (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) + //let t_s_factor = (1.-damage_0).powf(self.parameters.REDUCE_T); + let yield_factor = 1.0 - hardening; + let fracture_surface = yield_factor + * (self.parameters.A * (p_s + t_s).powf(self.parameters.N) * self.parameters.SIGMAHEL) .max(0.0); let residual_surface = (self.parameters.B * (p_s).powf(self.parameters.M) * self.parameters.SIGMAHEL).max(0.0); @@ -96,9 +121,9 @@ impl ConstitutiveModel for JH23D { } let yield_surface = rate_factor * { if damage_0 == 0.0 { - fracture_surface + fracture_surface * hardening_factor } else { - fracture_surface * (1. - damage_0) + damage_0 * residual_surface + fracture_surface * hardening_factor * (1. - damage_0) + damage_0 * residual_surface } }; if s_tr_eq > yield_surface { @@ -108,7 +133,15 @@ impl ConstitutiveModel for JH23D { del_lambda = (s_tr_eq - yield_surface) / (3. * self.parameters.SHEAR_MODULUS); alpha = yield_surface / s_tr_eq; - damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + if self.parameters.E_F > 0.0 { + //let lambda_old = -self.parameters.E_F * (1. - damage_0).ln(); + let lambda_new = lambda_old + del_lambda; + damage_1 = + (1. - ((self.parameters.E_0 - lambda_new) / self.parameters.E_F).exp()).max(0.0); + } else { + damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); + } + //damage_1 = (damage_0 + del_lambda / e_p_f).min(self.parameters.DMAX); output.set_scalar(Q::Damage, ip, damage_1); } else { alpha = 1.0; @@ -131,9 +164,11 @@ impl ConstitutiveModel for JH23D { output.set_scalar(Q::Density, ip, density_1); let mu = density_1 / self.parameters.RHO - 1.; - + let mut local_bulk_modulus = self.parameters.K1; let p_1 = { if mu > 0.0 { + local_bulk_modulus += + 2.0 * self.parameters.K2 * mu + 3.0 * self.parameters.K3 * mu.powi(2); self.parameters.K1 * mu + self.parameters.K2 * mu.powi(2) + self.parameters.K3 * mu.powi(3) @@ -144,6 +179,7 @@ impl ConstitutiveModel for JH23D { if p_trial > p_damaged { p_trial } else { + local_bulk_modulus = 0.0; p_damaged } } @@ -167,7 +203,9 @@ impl ConstitutiveModel for JH23D { // Calculate bulk viscosity let l = input.get_scalar(Q::CellDiameter, ip); let tr_deps = d_eps_vol * 3.; - let c = (self.parameters.K1 / self.parameters.RHO).sqrt(); + let bulk_modulus = local_bulk_modulus * self.parameters.LOCAL_SOUND_SPEED + + (1.0 - self.parameters.LOCAL_SOUND_SPEED) * self.parameters.K1; + let c = (bulk_modulus / self.parameters.RHO).sqrt(); let q_1 = { if tr_deps < 0.0 { l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) @@ -201,15 +239,14 @@ impl ConstitutiveModel for JH23D { } let elastic_rate = - - (1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; - + -(1. - alpha) / (2. * self.parameters.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + let density_mid = 0.5 * (density_0 + density_1); if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { let s_mid = 0.5 * (s_0 + s_1); let deviatoric_rate = d_eps_dev - elastic_rate; let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); - let e_1 = e_0 - + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); + let e_1 = e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); } if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { @@ -217,9 +254,8 @@ impl ConstitutiveModel for JH23D { let p_mid = -0.5 * (p_0 + p_1); let deviatoric_rate = elastic_rate; let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); - let e_1 = e_0 - + del_t / density_mid - * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + let e_1 = + e_0 + del_t / density_mid * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); output.set_scalar(Q::InternalElasticEnergy, ip, e_1); } if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { @@ -230,7 +266,7 @@ impl ConstitutiveModel for JH23D { } if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); - let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let q_mid = -0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); let e_1 = e_0 + del_t / density_mid * 3. * q_mid * d_eps_vol; output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); } diff --git a/src/jhr.rs b/src/jhr.rs new file mode 100644 index 0000000..b96c1f6 --- /dev/null +++ b/src/jhr.rs @@ -0,0 +1,278 @@ + +use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q}; +use crate::stress_strain::{ + deviatoric, mandel_decomposition, mandel_rate_from_velocity_gradient, MANDEL_IDENTITY, +}; + +use nalgebra::{SMatrix, DVectorView}; +use std::cmp; +use std::collections::HashMap; + +#[derive(Debug)] +pub struct JHR3D { + RHO: f64, + SHEAR_MODULUS: f64, + A: f64, + B: f64, + C: f64, + M: f64, + N: f64, + EPS0: f64, + T: f64, + SIGMAHEL: f64, + PHEL: f64, + D1: f64, + D2: f64, + K1: f64, + K2: f64, + K3: f64, + BETA: f64, + EFMIN: f64, + DMAX: f64, +} + +impl ConstitutiveModel for JHR3D { + fn new(parameters: &HashMap) -> Option { + Some(Self { + RHO: *parameters.get("RHO")?, + SHEAR_MODULUS: *parameters.get("SHEAR_MODULUS")?, + A: *parameters.get("A")?, + B: *parameters.get("B")?, + C: *parameters.get("C")?, + M: *parameters.get("M")?, + N: *parameters.get("N").unwrap(), + EPS0: *parameters.get("EPS0").unwrap(), + T: *parameters.get("T").unwrap(), + SIGMAHEL: *parameters.get("SIGMAHEL").unwrap(), + PHEL: *parameters.get("PHEL").unwrap(), + D1: *parameters.get("D1").unwrap(), + D2: *parameters.get("D2").unwrap(), + K1: *parameters.get("K1").unwrap(), + K2: *parameters.get("K2").unwrap(), + K3: *parameters.get("K3").unwrap(), + BETA: *parameters.get("BETA").unwrap(), + EFMIN: *parameters.get("EFMIN").unwrap(), + DMAX: *parameters.get("DMAX").unwrap_or(&1.0), + }) + } + fn evaluate_ip(&self, ip: usize, del_t: f64, input: &QValueInput, output: &mut QValueOutput) { + let velocity_gradient = input + .get_tensor::<{ Q::VelocityGradient.dim() }, { Q::VelocityGradient.size() }>( + Q::VelocityGradient, + ip, + ); + let d_eps = mandel_rate_from_velocity_gradient(&velocity_gradient); + let (mut d_eps_vol, d_eps_dev) = mandel_decomposition(&d_eps); + d_eps_vol *= -1.0; + + let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip); + + let mut del_lambda = 0.0; + + let damage_0 = input.get_scalar(Q::Damage, ip); + let mut damage_1 = damage_0; + + let (p_0, s_0) = mandel_decomposition(&sigma_0); + let p_0 = p_0 - input.get_scalar(Q::BulkViscosity, ip); + + let s_tr = s_0 + 2. * self.SHEAR_MODULUS * d_eps_dev * del_t; + let s_tr_eq = (1.5 * s_tr.norm_squared()).sqrt(); + let d_eps_eq = ((2. / 3.) * d_eps.norm_squared()).sqrt(); + let mut alpha; + + let p_s = p_0 / self.PHEL; + let t_s = self.T / self.PHEL; + let mut rate_factor = 1.; + + let fracture_surface = + (self.A * (p_s + t_s*(1.-damage_0)).powf(self.N) * self.SIGMAHEL) + .max(0.0); + let residual_surface = + (self.B * (p_s).powf(self.M) * self.SIGMAHEL).max(0.0); + if d_eps_eq >= self.EPS0 { + rate_factor += self.C * (d_eps_eq / self.EPS0).ln(); + } + let yield_surface = rate_factor * { + if damage_0 == 0.0 { + fracture_surface + } else { + fracture_surface * (1. - damage_0) + damage_0 * residual_surface + } + }; + if s_tr_eq > yield_surface { + let e_p_f = (self.D1 * (p_s + t_s).powf(self.D2)) + .max(self.EFMIN); + + del_lambda = (s_tr_eq - yield_surface) / (3. * self.SHEAR_MODULUS); + alpha = yield_surface / s_tr_eq; + + damage_1 = (damage_0 + del_lambda / e_p_f).min(self.DMAX); + output.set_scalar(Q::Damage, ip, damage_1); + } else { + alpha = 1.0; + } + + // /*********************************************************************** + // * UPDATE DENSITY + // * The density is updated using the explicit midpoint rule for the + // * deformation gradient. + // TODO: Move this since, it will be calculated outside of the constitutive model + // **********************************************************************/ + let f1 = del_t / 2. * 3. * d_eps_vol; + let density_0 = input.get_scalar(Q::Density, ip); + let density_1 = density_0 * (1. - f1) / (1. + f1); + debug_assert!( + density_1 > 0.0, + "Negative density encountered in JH2 model: {}", + density_1 + ); + output.set_scalar(Q::Density, ip, density_1); + + let mu = density_1 / self.RHO - 1.; + + let p_1 = { + if mu > 0.0 { + self.K1 * mu + + self.K2 * mu.powi(2) + + self.K3 * mu.powi(3) + //+ input.get_scalar(Q::BulkingPressure, ip) + } else { + let p_trial = self.K1 * mu; + let p_damaged = -self.T * (1. - damage_1); + if p_trial > p_damaged { + p_trial + } else { + p_damaged + } + } + }; + + // Calculate bulk viscosity + let l = input.get_scalar(Q::CellDiameter, ip); + let tr_deps = d_eps_vol * 3.; + let c = (self.K1 / self.RHO).sqrt(); + let q_1 = { + if tr_deps < 0.0 { + l * density_1 * (1.5 * l * tr_deps.powi(2) + c * 0.06 * tr_deps.abs()) + } else { + 0.0 + } + }; + output.set_scalar(Q::BulkViscosity, ip, q_1); + // /*********************************************************************** + // * Combine deviatoric and volumetric stresses + // **********************************************************************/ + let s_1 = s_tr * alpha; + let sigma_1 = s_1 - MANDEL_IDENTITY * (p_1 + q_1); + output.set_vector(Q::MandelStress, ip, sigma_1); + + // *********************************************************************** + // Update optional output variables if needed + // ********************************************************************** + + if output.is_some(Q::EqStrainRate) { + output.set_scalar(Q::EqStrainRate, ip, d_eps_eq); + } + if output.is_some(Q::MandelStrainRate) { + output.set_vector(Q::MandelStrainRate, ip, d_eps); + } + if output.is_some(Q::MisesStress) { + output.set_scalar(Q::MisesStress, ip, alpha * s_tr_eq); + } + if output.is_some(Q::Pressure) { + output.set_scalar(Q::Pressure, ip, p_1); + } + + let elastic_rate = + - (1. - alpha) / (2. * self.SHEAR_MODULUS * del_t) * s_0 + alpha * d_eps_dev; + + let density_mid = 0.5 * (density_0 + density_1); + if output.is_some(Q::InternalPlasticEnergy) && input.is_some(Q::InternalPlasticEnergy) { + let s_mid = 0.5 * (s_0 + s_1); + let deviatoric_rate = d_eps_dev - elastic_rate; + let e_0 = input.get_scalar(Q::InternalPlasticEnergy, ip); + let e_1 = e_0 + + del_t / density_mid * (s_mid.dot(&deviatoric_rate)); + output.set_scalar(Q::InternalPlasticEnergy, ip, e_1); + } + if output.is_some(Q::InternalElasticEnergy) && input.is_some(Q::InternalElasticEnergy) { + let s_mid = 0.5 * (s_0 + s_1); + let p_mid = -0.5 * (p_0 + p_1); + let deviatoric_rate = elastic_rate; + let e_0 = input.get_scalar(Q::InternalElasticEnergy, ip); + let e_1 = e_0 + + del_t / density_mid + * (s_mid.dot(&deviatoric_rate) + 3. * d_eps_vol * p_mid); + output.set_scalar(Q::InternalElasticEnergy, ip, e_1); + } + if output.is_some(Q::InternalEnergy) && input.is_some(Q::InternalEnergy) { + let e_0 = input.get_scalar(Q::InternalEnergy, ip); + let sigma_mid = 0.5 * (sigma_0 + sigma_1); + let e_1 = e_0 + del_t / density_mid * sigma_mid.dot(&d_eps); + output.set_scalar(Q::InternalEnergy, ip, e_1); + } + if output.is_some(Q::InternalHeatingEnergy) && input.is_some(Q::InternalHeatingEnergy) { + let e_0 = input.get_scalar(Q::InternalHeatingEnergy, ip); + let q_mid = - 0.5 * (q_1 + input.get_scalar(Q::BulkViscosity, ip)); + let e_1 = e_0 + del_t / density_mid * 3. * q_mid * d_eps_vol; + output.set_scalar(Q::InternalHeatingEnergy, ip, e_1); + } + if output.is_some(Q::EqPlasticStrain) && input.is_some(Q::EqPlasticStrain) { + output.set_scalar( + Q::EqPlasticStrain, + ip, + input.get_scalar(Q::EqPlasticStrain, ip) + del_lambda, + ); + } + } + + /// Returns the physical quantities that are required as input for the + /// constitutive model together with their dimensions. + fn define_input(&self) -> HashMap { + HashMap::from([ + (Q::VelocityGradient, QDim::SquareTensor(3)), + (Q::CellDiameter, QDim::Scalar), + ]) + } + + /// Returns the physical quantities that are needed as internal variables + /// for the constitutive model together with their dimensions. These Variables are + /// stored both in in the input and the output. + fn define_history(&self) -> HashMap { + HashMap::from([ + (Q::MandelStress, QDim::Vector(6)), + (Q::Damage, QDim::Scalar), + (Q::Density, QDim::Scalar), + (Q::BulkViscosity, QDim::Scalar), + ]) + } + + /// Returns the physical quantities that are needed as output, but are not + /// necessarily needed in oredr to calculate the constitutive model. An example is + /// the consistent tangent which is not needed for the calculation of the stresses + /// and is therefore purely an output quantity. + fn define_output(&self) -> HashMap { + HashMap::from([(Q::MandelStress, QDim::Vector(6))]) + } + + /// Returns the physical quantities that are optional output of the constitutive + /// model. These quantities are not needed for the calculation of the stresses + /// but can be useful for postprocessing. + fn define_optional_output(&self) -> HashMap { + HashMap::from([ + (Q::EqStrainRate, QDim::Scalar), + (Q::MandelStrainRate, QDim::Vector(6)), + (Q::MisesStress, QDim::Scalar), + (Q::Pressure, QDim::Scalar), + ]) + } + fn define_optional_history(&self) -> HashMap { + HashMap::from([ + (Q::EqPlasticStrain, QDim::Scalar), + (Q::InternalPlasticEnergy, QDim::Scalar), + (Q::InternalElasticEnergy, QDim::Scalar), + (Q::InternalEnergy, QDim::Scalar), + (Q::InternalHeatingEnergy, QDim::Scalar), + ]) + } +} diff --git a/src/lib.rs b/src/lib.rs index 8893a9b..1fb82ca 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,10 +3,12 @@ use std::collections::HashMap; use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q}; use crate::jh2::JH23D; +use crate::jhr::JHR3D; //use crate::jh_concrete::JHConcrete3D; use crate::generic_jh2::GenericJH23D; use crate::gradient_jh2::GradientJH23D; use crate::smallstrain::linear_elastic::LinearElastic3D; +use crate::smallstrain::{evaluate_model, elasticity_3d}; use crate::hypoelasticity::Hypoelasticity3D; //use crate::stress_strain; use nalgebra::{Const, DVectorView, DVectorViewMut, Dyn, SMatrix}; @@ -20,6 +22,7 @@ pub mod jh2; //pub mod jh_concrete; pub mod generic_jh2; pub mod gradient_jh2; +pub mod jhr; pub mod hypoelasticity; pub mod smallstrain; pub mod stress_strain; @@ -393,7 +396,23 @@ fn py_jaumann_rotation_expensive( stress_strain::jaumann_rotation_expensive(del_t, &velocity_gradient, &mut stress); Ok(()) } - +// #[pyfunction(name="evaluate_elasticity_3d")] +// fn py_evaluate_elasticity_3d( +// del_t: f64, +// stress: PyReadwriteArray1, +// del_strain: PyReadonlyArray1, +// parameters: PyReadonlyArray1, +// history: PyReadwriteArray1, +// tangent: PyReadwriteArray1, +// ) -> PyResult<()> { +// let stress = stress.as_slice_mut()?; +// let del_strain = del_strain.as_slice()?; +// let parameters = parameters.as_slice()?; +// let history = history.as_slice_mut()?; +// let mut tangent = tangent.as_slice_mut()?; +// evaluate_model::<6,0,2>(&elasticity_3d, del_t, stress, del_strain, parameters, history, tangent); +// Ok(()) +// } #[pymodule] fn comfe(_py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; diff --git a/src/smallstrain.rs b/src/smallstrain.rs index 8d65904..1561ed5 100644 --- a/src/smallstrain.rs +++ b/src/smallstrain.rs @@ -1,44 +1,108 @@ -use nalgebra::{DVectorView, DVectorViewMut}; +use nalgebra::{DVectorView, DVectorViewMut, SMatrix, SVector}; +use std::convert::TryInto; pub mod linear_elastic; -pub trait SmallStrainModel { - //fn new(parameters: &HashMap) -> Self; - fn evaluate_ip( - &self, - ip: usize, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, +// define a type for a specific function signature +pub type SmallStrainModelFn< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +> = fn( + del_t: f64, + stress: [f64; STRESS_STRAIN_DIM], + del_strain: [f64; STRESS_STRAIN_DIM], + parameters: [f64; PARAMETER_DIM], + history: [f64; HISTORY_DIM], +) -> ( + [f64; STRESS_STRAIN_DIM], + [[f64; STRESS_STRAIN_DIM]; STRESS_STRAIN_DIM], + [f64; HISTORY_DIM], +); + +pub fn elasticity_3d( + del_t: f64, + stress: [f64; 6], + del_strain: [f64; 6], + parameters: [f64; 2], + history: [f64; 0], +) -> ([f64; 6], [[f64; 6]; 6], [f64; 0]) { + let E = parameters[0]; + let nu = parameters[1]; + let mu = E / (2.0 * (1.0 + nu)); + let lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + let c1 = lambda + 2.0 * mu; + let c2 = 2. * mu; + let D = SMatrix::::from([ + [c1, lambda, lambda, 0.0, 0.0, 0.0], + [lambda, c1, lambda, 0.0, 0.0, 0.0], + [lambda, lambda, c1, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, c2, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, c2, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, c2], + ]); + let stress_vector = SVector::::from(stress); + let del_strain_vector = SVector::::from(del_strain); + let new_stress = stress_vector + D * del_strain_vector; + let D_out = D.data.0; + let new_stress = new_stress.data.0[0]; + let history = history; + (new_stress, D_out, history) +} +pub struct SmallStrainModel< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +> { + pub model: SmallStrainModelFn, + pub parameters: [f64; PARAMETER_DIM], +} + +pub fn evaluate_model< + const STRESS_STRAIN_DIM: usize, + const HISTORY_DIM: usize, + const PARAMETER_DIM: usize, +>( + model: &SmallStrainModelFn, + del_t: f64, + stress: &mut [f64], + del_strain: &[f64], + parameters: &[f64], + history: &mut [f64], + tangent: &mut [f64], +) { + let parameters: [f64; PARAMETER_DIM] = parameters + .try_into() + .expect("Slice length does not match array length"); + + let stress_chunks = stress.chunks_exact_mut(STRESS_STRAIN_DIM); + let del_strain_chunks = del_strain.chunks_exact(STRESS_STRAIN_DIM); + let history_chunks = history.chunks_exact_mut(HISTORY_DIM); + let tangent_chunks = tangent.chunks_exact_mut(STRESS_STRAIN_DIM.pow(2)); + + assert!( + stress_chunks.len() == del_strain_chunks.len() + && stress_chunks.len() == history_chunks.len() + && stress_chunks.len() == tangent_chunks.len() ); - fn evaluate( - &self, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, - ) { - assert_eq!(stress.nrows(), del_strain.nrows()); - let n: usize = stress.nrows() / 6; - assert_eq!(stress.nrows(), n * 6); - for ip in 0..n { - self.evaluate_ip(ip, del_t, stress, del_strain, tangent); - } - } - fn evaluate_some( - &self, - del_t: f64, - stress: &mut DVectorViewMut, - del_strain: &DVectorView, - tangent: &mut DVectorViewMut, - ips: &[usize], - ) { - assert_eq!(stress.nrows(), del_strain.nrows()); - let n: usize = stress.nrows() / 6; - assert_eq!(stress.nrows(), n * 6); - for ip in ips { - self.evaluate_ip(*ip, del_t, stress, del_strain, tangent); - } + for (((stress_chunk, del_strain_chunk), history_chunk), tangent_chunk) in stress_chunks + .zip(del_strain_chunks) + .zip(history_chunks) + .zip(tangent_chunks) + { + let stress: [f64; STRESS_STRAIN_DIM] = stress_chunk + .try_into() + .expect("Slice length does not match array length"); + let del_strain: [f64; STRESS_STRAIN_DIM] = del_strain_chunk + .try_into() + .expect("Slice length does not match array length"); + let history: [f64; HISTORY_DIM] = history_chunk + .try_into() + .expect("Slice length does not match array length"); + let (new_stress, D, new_history) = model(del_t, stress, del_strain, parameters, history); + stress_chunk.copy_from_slice(&new_stress); + history_chunk.copy_from_slice(&new_history); + tangent_chunk.copy_from_slice(unsafe{std::slice::from_raw_parts(D.as_ptr() as *const f64, STRESS_STRAIN_DIM.pow(2))}); } + } diff --git a/test/test_jh2_2d.py b/test/test_jh2_2d.py index adc7659..c72ee92 100644 --- a/test/test_jh2_2d.py +++ b/test/test_jh2_2d.py @@ -503,11 +503,6 @@ def test_single_element_2d(test_case: dict, plot: str | None = None) -> None: p_ = np.array(p_).reshape((-1, 1)) s_eq_ = np.array(s_eq_).reshape((-1, 1)) - points = np.hstack((p_, s_eq_)) - tree = KDTree(points) - distances = tree.query(test_case["points"]) - assert np.mean(distances[0] / np.max(np.abs(test_case["points"][:, 1]))) < 0.05 - if plot is not None: import matplotlib import matplotlib.pyplot as plt @@ -524,7 +519,13 @@ def test_single_element_2d(test_case: dict, plot: str | None = None) -> None: plt.savefig(f"{plot}.png") plt.clf() + points = np.hstack((p_, s_eq_)) + tree = KDTree(points) + distances = tree.query(test_case["points"]) + assert np.mean(distances[0] / np.max(np.abs(test_case["points"][:, 1]))) < 0.05 + if __name__ == "__main__": test_single_element_2d(case_3a, plot="jh2_case_3a") test_single_element_2d(case_1, plot="jh2_case_1") + test_single_element_2d(case_2, plot="jh2_case_2")