diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 964e3a7..685c94c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,7 +13,7 @@ jobs: - uses: actions/checkout@v3 - uses: actions/setup-python@v4 with: - python-version: "3.10" + python-version: "3.11" cache: pip - name: Run pre-commit uses: pre-commit/action@v3.0.0 @@ -26,7 +26,7 @@ jobs: - ubuntu-latest # - windows-latest Reason: get it running on ubuntu-latest first # - macOS-latest Reason: get it running on ubuntu-latest first - python-version: ["3.10",] + python-version: ["3.11",] name: Run tests runs-on: ${{ matrix.os }} steps: diff --git a/.gitignore b/.gitignore index b0d6253..a3faa7b 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ __pycache__/ *.py[cod] *$py.class +/bin # C extensions *.so diff --git a/README.md b/README.md index 1e2691c..604a2e6 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ The corresponding implementation is available in our [AdaptiveThermoMechROM](htt ### Online phase: Usage of the thermo-mechanical NTFA in simulations on the macroscale -1. Load the tabular data for the NTFA matrices $`\underline{\underline{A}}(\theta_j)`$, $`\underline{\underline{D}}(\theta_j)`$, $`\bar{\underline{\underline{C}}}(\theta_j)`$, $`\underline{\tau}\_{\mathrm{\theta}}(\theta_j)`$, and $`\underline{\tau}_{\mathsf{p}}(\theta_j)`$ that are generated in the offline phase based on direct numerical simulations on the microscale. +1. Load the tabular data for the NTFA matrices $`\underline{\underline{A}}(\theta_j)`$, $`\underline{\underline{D}}(\theta_j)`$, $`\bar{\underline{\underline{C}}}(\theta_j)`$, $`\underline{\tau}_{\mathrm{\theta}}(\theta_j)`$, and $`\underline{\tau}_{\mathsf{p}}(\theta_j)`$ that are generated in the offline phase based on direct numerical simulations on the microscale. Optionally truncate the NTFA modes $`N_{\mathrm{modes}}`$ to be used. 2. Perform a linear interpolation to determine the NTFA matrices at the current model temperature based on the tabular data. @@ -86,7 +86,11 @@ pip install -e . ### Requirements -TODO: update +TODO: update dependencies + +TODO: upload dataset to Darus + +TODO: provide functionality for download from Darus - Python 3.9 or later - `numpy` and `h5py` (installed as part of the `thermontfa` PIP package) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..f9d2932 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,19 @@ +# Building the ThermoNTFA Python documentation + +To build the documentation: + +1. Install ThermoNTFA (`thermontfa` pip package). It must be possible to import + the module ``thermontfa``. +2. Run in this directory: + + make html + +## Processing demo/example programs + +Python demo programs are written in Python, with MyST-flavoured +Markdown syntax for comments. + +1. `jupytext` reads the Python demo code, then converts to and writes a + Markdown file. +2. `myst_parser` allows Sphinx to process the Markdown file. + diff --git a/test.h5 b/test.h5 new file mode 100644 index 0000000..2e623a8 Binary files /dev/null and b/test.h5 differ diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 8708b8b..20b4679 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -3,36 +3,35 @@ Test tabular interpolation """ -# %% import h5py import numpy as np - -# %% from thermontfa import TabularInterpolation -# %% def test_tab_inter_init(): + theta = np.linspace(0, 1, 10) + val = np.linspace(2, 5, 10).reshape((-1, 1)) + + tab_inter = TabularInterpolation(temps=theta, data=val) + assert tab_inter.temp_min == theta.min() + assert tab_inter.temp_max == theta.max() + assert np.allclose(tab_inter.temps.ravel(), theta.ravel()) + assert np.allclose(tab_inter.data.ravel(), val.ravel()) + + +def test_tab_inter_from_h5(): theta = np.arange(5) val = theta**2 idx = np.random.permutation(np.arange(5)) - theta = theta[idx] - val = val[idx] - val.reshape((5, 1)) - F = h5py.File("test.h5", "w") - F.create_dataset("/theta", data=theta) - F.create_dataset("/data", data=val) - F.close() - - tab_inter = TabularInterpolation() - tab_inter.init_h5("test.h5", "/theta", "/data") - print(tab_inter.t_min) - print(tab_inter.t_max) - print(tab_inter.t) - print(tab_inter.data) + theta_perm = theta[idx] + val_perm = val[idx] + val_perm = val_perm.reshape((-1, 1)) + with h5py.File("test.h5", "w") as file: + file.create_dataset("/theta", data=theta_perm) + file.create_dataset("/data", data=val_perm) - -# %% -a = np.linspace(2, 5, 10).reshape((-1, 1)) -t = np.linspace(0, 1, 10) -inter = TabularInterpolation(a, t) + tab_inter = TabularInterpolation.from_h5(file_name="test.h5", dset_temps="/theta", dset_data="/data") + assert tab_inter.temp_min == theta.min() + assert tab_inter.temp_max == theta.max() + assert np.allclose(tab_inter.temps.ravel(), theta.ravel()) + assert np.allclose(tab_inter.data.ravel(), val.ravel()) diff --git a/thermontfa/tabular_interpolation.py b/thermontfa/tabular_interpolation.py index 1c24bb8..02c768b 100644 --- a/thermontfa/tabular_interpolation.py +++ b/thermontfa/tabular_interpolation.py @@ -19,7 +19,7 @@ 390740016 (EXC-2075); 406068690 (FR2702/8-1); 517847245 (FR2702/10-1). """ -from typing import Optional, Tuple +from typing import Optional, Tuple, Self import h5py import numpy as np @@ -27,117 +27,130 @@ class TabularInterpolation: """ - Tabular interpolation for the NTFA + Tabular interpolation for the thermo-mechanical NTFA + + Performs a linear interpolation of the NTFA system matrices for a given temperature + given tabular data for sufficiently many temperature points. + It can be initialized for given data or based on a HDF5 file. """ - t_min: float = 0.0 - t_max: float = 0.0 + temp_min: float = 0.0 + temp_max: float = 0.0 dim: Tuple[int, ...] = () - t: np.ndarray = np.array([]) + temps: np.ndarray = np.array([]) data: np.ndarray = np.array([]) const_extrapolate: bool = False def __init__( - self, data: Optional[np.ndarray] = None, param: Optional[np.ndarray] = None - ): + self, temps: np.ndarray = None, data: np.ndarray = None, const_extrapolate: bool = False + ) -> None: """ - Initialize the tabular interpolator - - :param data: - :param param: + Initialize the tabular interpolator for given `data` at prescribed temperatures `temps`. + + :param temps: temperature points on which tabular data is available. + The shape of the `numpy` array is expected to be `(N_t,)`. + :type temps: np.ndarray + :param data: tabular data, e.g., a batch of NTFA system matrices with shape `(N_t, ...)`. + :type data: np.ndarray + :param const_extrapolate: If true, a constant extrapolation instead of a linear extrapolation is performed. + The default value is false. + :type const: bool """ - if data is not None and param is not None: - idx = np.argsort(param) + if temps is None or data is None: + raise ValueError("The arguments `temps` and `data` are required.") + + idx = np.argsort(temps) + + self.data = data[idx, ...] + self.temps = temps[idx] + self.temp_min = self.temps[0] + self.temp_max = self.temps[-1] + self.dim = data[0].shape + self.const_extrapolate = const_extrapolate - self.data = data[idx, :] - self.t = param[idx] - self.t_min = self.t[0] - self.t_max = self.t[-1] - self.dim = data[0].shape + n = self.temps.shape[0] - def init_h5( - self, + if data.ndim == 1: + self.data = self.data.reshape((-1, 1)) + assert ( + self.temps.size == n + ), "Parameter must be an array of scalars (i.e. shape=[n] or [n, 1] or [n, 1, 1] or ...)" + assert ( + self.data.shape[0] == n + ), f"Number of scalar parameters not matching dimension of available data ({n} vs. {self.data.shape[0]}." + + @classmethod + def from_h5( + cls, file_name: str, - dset_param: str, + dset_temps: str, dset_data: str, - extrapolate: str = "linear", - transpose: Optional[Tuple[int, ...]] = None, - ): + transpose_dims: Optional[Tuple[int, ...]] = None, + const_extrapolate: bool = False, + ) -> Self: """ - Initialize the tabular interpolator using data in a H5 file - - :param file_name: path of the H5 file + Initialize the tabular interpolator based on tabular data stored in a HDF5 file (*.h5). + + This is a factory method and returns a new instance of the `TabularInterpolation` class. + It is expected that the HDF5 file contains a data set with path `dset_temps` that contains + a list of the temperature points on which tabular data is available. The shape of this + dataset is expected to be `(N_t,)`. + Additionaly, the HDF5 file must contain a data set with path `dset_data` that contains the + tabular data, e.g., a batch of NTFA system matrices with shape `(N_t, ...)`. + The order of axes/dimensions of the data set with path `dset_data` can be changed by transposing + to the axis order given in `transpose_dims`. + + :param file_name: path of the HDF5 file :type file_name: str - :param dset_param: path to the desired dataset in the H5 file + :param dset_param: path to the desired dataset in the HDF5 file :type dset_param: str :param dset_data: :type dset_param: str - :param extrapolate: "linear" or "constant" - :type extrapolate: str - :param transpose: axis order for transposition - :type transpose: Tuple[int, ...], optional + :param const_extrapolate: "linear" or "constant" + :type extrapolate: bool + :param transpose_dims: axis order for transposition + :type transpose_dims: Tuple[int, ...], optional + :return: new instance of the `TabularInterpolation` class + :rtype: TabularInterpolation """ - F = h5py.File(file_name, "r") - self.t = np.array(F[dset_param]) - if transpose is None: - self.data = np.array(F[dset_data]) - else: - self.data = np.array(F[dset_data]).transpose(transpose) - F.close() - self.dim = self.data.shape[1:] - - if extrapolate == "linear": - self.const_extrapolate = False - else: - if extrapolate == "constant": - self.const_extrapolate = True + with h5py.File(file_name, "r") as file: + temps = np.array(file[dset_temps]) + + if transpose_dims is None: + data = np.array(file[dset_data]) else: - raise ValueError( - f'extrapolate can only take values "linear" or "constant" (received: "{extrapolate}"' - ) + data = np.array(file[dset_data]).transpose(transpose_dims) + + return cls(temps=temps, data=data, const_extrapolate=const_extrapolate) - if self.data.ndim == 1: - self.data = self.data.reshape((-1, 1)) - n = self.t.shape[0] - assert ( - self.t.size == n - ), "ERROR: parameter must be an array of scalars (i.e. shape=[n] or [n, 1] or [n, 1, 1] or ...)" - assert ( - self.data.shape[0] == n - ), f"ERROR: number of scalar parameters not matching dimension of available data ({n} vs. {self.data.shape[0]}." - idx = np.argsort(self.t) - self.t = self.t[idx] - self.data = self.data[idx, :] - self.t_min = self.t[0] - self.t_max = self.t[-1] - - def interpolate(self, t: float) -> np.ndarray: + def interpolate(self, temp: float) -> np.ndarray: """ - Perform a linear interpolation at a given temperature t + Perform a linear interpolation based on the available tabular data at a given temperature `temp` - :param t: temperature point for interpolation - :type t: float + :param temp: temperature point for interpolation + :type temp: float :return: interpolated quantity :rtype: np.ndarray """ - if t < self.t_min: + if temp < self.temp_min: if self.const_extrapolate: - return self.data[0, :] + return self.data[0, ...] else: # linear extrapolation - alpha = (self.t_min - t) / (self.t[1] - self.t[0]) - return self.data[0, :] - alpha * (self.data[1, :] - self.data[0, :]) + alpha = (self.temp_min - temp) / (self.temp[1] - self.temp[0]) + return self.data[0, ...] - alpha * (self.data[1, ...] - self.data[0, ...]) - if t >= self.t_max: + if temp >= self.temp_max: if self.const_extrapolate: - return self.data[-1, :] + return self.data[-1, ...] else: # linear extrapolation - alpha = (t - self.t_max) / (self.t[-1] - self.t[-2]) - return self.data[-1, :] + alpha * (self.data[-1, :] - self.data[-2, :]) + alpha = (temp - self.temp_max) / (self.temp[-1] - self.temp[-2]) + return self.data[-1, ...] + alpha * (self.data[-1, ...] - self.data[-2, ...]) - idx = np.searchsorted(self.t > t, 1) - 1 + idx = np.searchsorted(self.temp > temp, 1) - 1 t1 = self.t[idx] t2 = self.t[idx + 1] - alpha = (t - t1) / (t2 - t1) - return self.data[idx, :] + alpha * (self.data[idx + 1, :] - self.data[idx, :]) + alpha = (temp - t1) / (t2 - t1) + interp_data = self.data[idx, ...] + alpha * (self.data[idx + 1, ...] - self.data[idx, ...]) + return interp_data diff --git a/thermontfa/thermoNTFA.py b/thermontfa/thermoNTFA.py index 0c4ee26..1d99344 100644 --- a/thermontfa/thermoNTFA.py +++ b/thermontfa/thermoNTFA.py @@ -38,6 +38,8 @@ class ThermoMechNTFA: with temperature-dependent material parameters in both phases. """ + # TODO: Document properties + def __init__( self, file_name: str, @@ -46,12 +48,12 @@ def __init__( N_max: Optional[int] = None, tol: float = 1e-4, verbose: bool = False, - ): + ) -> None: """ - Initialize the thermo-mechanical NTFA from an H5 file + Initialize the thermo-mechanical NTFA from an HDF5 file (*.h5) - Seek the data in H5-file named `file_name` within the group `group_name`. - The following datasets containing tabular data for the NTFA are expected: + Seek the data in HDF5 file named `file_name` within the group `group_name`. + The following data sets containing tabular data for the NTFA are expected in the group: - `temperatures`: list of temperature points of the tabular data, shape: `(N_temp,)` - `A_bar`: shape: `(N_temp, N_modes, 6)` - `C_bar`: shape: `(N_temp, 6, 6)` @@ -59,9 +61,13 @@ def __init__( - `tau_theta`: shape: `(N_temp, 6)` - `tau_xi`: shape: `(N_temp, N_modes)` - :param file_name: path to the H5 file + In addition, the group in the HDF5 file must contain the following data sets: + - `v_frac`: volume fraction of the different phases + - `SIG_phases`: stress data different phases + + :param file_name: path to the HDF5 file :type file_name: str - :param group_name: group in the H5 file that contains the NTFA tabular data + :param group_name: group in the HDF5 file that contains the NTFA tabular data :type group_name: str :param sig_y: function/callable that returns the yield stress `sig_y(theta, q_n, derivative)` given the temperature `theta` and the current isotropic hardening variable `q_n`. @@ -79,57 +85,55 @@ def __init__( self.N_max = N_max self.verbose = verbose - self.A_bar = TabularInterpolation() - self.A_bar.init_h5( - self.file_name, - self.group_name + "/temperatures", - self.group_name + "/A_bar", + self.A_bar = TabularInterpolation.from_h5( + file_name=self.file_name, + dset_temps=self.group_name + "/temperatures", + dset_data=self.group_name + "/A_bar", ) - self.C_bar = TabularInterpolation() - self.C_bar.init_h5( - self.file_name, - self.group_name + "/temperatures", - self.group_name + "/C_bar", + self.C_bar = TabularInterpolation.from_h5( + file_name=self.file_name, + dset_temps=self.group_name + "/temperatures", + dset_data=self.group_name + "/C_bar", ) - self.D_xi = TabularInterpolation() - self.D_xi.init_h5( - self.file_name, self.group_name + "/temperatures", self.group_name + "/D_xi" + self.D_xi = TabularInterpolation.from_h5( + file_name=self.file_name, + dset_temps=self.group_name + "/temperatures", + dset_data=self.group_name + "/D_xi", ) - self.tau_th = TabularInterpolation() - self.tau_th.init_h5( - self.file_name, - self.group_name + "/temperatures", - self.group_name + "/tau_theta", + self.tau_theta = TabularInterpolation.from_h5( + file_name=self.file_name, + dset_temps=self.group_name + "/temperatures", + dset_data=self.group_name + "/tau_theta", ) - self.tau_xi = TabularInterpolation() - self.tau_xi.init_h5( - self.file_name, - self.group_name + "/temperatures", - self.group_name + "/tau_xi", + self.tau_xi = TabularInterpolation.from_h5( + file_name=self.file_name, + dset_temps=self.group_name + "/temperatures", + dset_data=self.group_name + "/tau_xi", ) - with h5py.File(self.file_name, "r") as F: - self.v_frac = np.array(F[self.group_name + "/v_frac"]) - sig_phases_data = np.array(F[self.group_name + "/SIG_phases"]) + with h5py.File(self.file_name, "r") as file: + self.v_frac = np.array(file[self.group_name + "/v_frac"]) + sig_phases_data = np.array(file[self.group_name + "/SIG_phases"]) - N_input = self.A_bar.data.shape[1] + N_input = self.A_bar.dim[0] if self.N_max is not None: # truncate modes if needed assert ( self.N_max <= N_input ), f"N_modes on input ({N_input}) is smaller than the provided N_max ({self.N_max})" + # TODO: truncate in tabular interpolation? self.A_bar.data = self.A_bar.data[:, : self.N_max, :] self.D_xi.data = self.D_xi.data[:, : self.N_max, : self.N_max] self.tau_xi.data = self.tau_xi.data[:, : self.N_max] sig_phases_data = sig_phases_data[:, :, :, : (self.N_max + 7)] self.sig_phases = [] - for d in sig_phases_data: - self.sig_phases.append(TabularInterpolation(self.A_bar.t, d)) + for sig_data in sig_phases_data: + self.sig_phases.append(TabularInterpolation(self.A_bar.temps, sig_data)) self.theta_i = 0.0 # last interpolation temperature self.C = self.C_bar.interpolate(self.theta_i) @@ -139,7 +143,7 @@ def __init__( self.t_xi = self.tau_xi.interpolate(self.theta_i) self.n_modes = self.D_xi.data.shape[-1] - def interpolate(self, theta: float): + def interpolate(self, theta: float) -> None: """ Interpolate NTFA matrices to current temperature `theta` if the given tolerance is exceeded @@ -181,6 +185,100 @@ def stress( zeta = np.hstack((eps, 1, xi)) return self.sig_phases[i_phase].interpolate(theta) @ zeta + def UMAT_mixed( + self, + eps_idx: np.ndarray, + eps_n: np.ndarray, + deps: np.ndarray, + sig_bc: np.ndarray, + theta: float, + q_n: float, + xi_n: np.ndarray, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float, np.ndarray]: + """ + Run the UMAT using partial eps-BC, e.g., uniaxial tension test. + + :param eps_idx: Indices of eps that are prescribed. If `None`, then all components of `sig_bc` are prescribed. + :type eps_idx: nd-array, dtype=int + :param eps_n: Strain at the beginning of the increment. + :type eps_n: nd-array, dtype=float, shape=[6] + :param deps: Strain increment. Only deps[eps_idx] is used. + :type deps: nd-array, dtype=float, shape=[6] + :param sig_bc: Stress at the end of the increment. Only non `eps_idx` components are used. + :type sig_bc: nd-array, dtype=float, shape=[6] + :param theta: Temperature at the end of the time increment. + :type theta: float + :param q_n: Hardening variable at the beginning of the time increment. + :type q_n: float + :param xi_n: Reduced coefficients at the beginning of the time increment. + :type xi_n: nd-array, dtype=float + :return: **E** - + Full strain tensor at the end of the increment, i.e. + the entries not within eps_idx are set + :rtype: np.ndarray, dtype=float, shape=[6] + :return: **S** - + Stress at the end of the increment. Only S[eps_idx] is non-zero. + :rtype: np.ndarray, dtype=float + :return: **C** - + Stiffness tensor at the end of the increment + :rtype: `np.ndarray`, `dtype=float`, `shape=[6, 6]` + :return: **q** - + Hardening variable at the end of the time increment + :rtype: float + :return: **xi** - + Reduced coefficients at the end of the time increment + :rtype: `nd.ndarray`, `dtype=float` + """ + q = q_n + # partial strain update: + sig_idx = np.setdiff1d(np.arange(6), eps_idx) + + eps = eps_n.copy() + eps[eps_idx] += deps[eps_idx] + deps[sig_idx] = 0.0 + + self.tol_sig = 1e-7 * self.sig_y(theta, q_n, derivative=False) + err_sig = 1000.0 * self.tol_sig + it = 0 + it_max = 100 + + # eleminate thermo-elastic effects and approximate deps + sig = self.stress(eps_n, theta, xi_n) + + # initial guess for eps asserting elastic behaviour + C = self.C + deps[sig_idx] += np.linalg.solve( + C[sig_idx][:, sig_idx], + -(sig[sig_idx] - sig_bc[sig_idx]) - C[sig_idx, :] @ deps, + ) + + # compute the actual stress state + eps = eps_n + deps + sig, q, xi, C = self.solve(eps, deps, theta, q_n, xi_n) + err_sig = np.linalg.norm(sig[sig_idx]) + while (err_sig > self.tol_sig) and (it < it_max): + ddeps = -np.linalg.solve( + C[sig_idx][:, sig_idx], sig[sig_idx] - sig_bc[sig_idx] + ) + factor = 1.0 + deps[sig_idx] += factor * ddeps + eps = eps_n + deps + sig, q, xi, C = self.solve(eps, deps, theta, q_n, xi_n) + err_sig = np.linalg.norm(sig[sig_idx]) + it = it + 1 + + if self.verbose: + err_eps = np.linalg.norm(ddeps) + print(f'it {it:3d} .. err_eps {err_eps:10.3e} (self.tol: {self.tol_eps:8.3e}) " + + ".. err_sig {err_sig:10.3e} (self.tol: {self.tol_sig:8.3e})') + + eps = eps_n + deps + + if self.verbose: + print(f'needed {it} iterations...') + + return eps, sig, C, q, xi + def solve( self, eps: np.ndarray, @@ -303,89 +401,3 @@ def solve( # self.dxi_deps = np.zeros((self.n_modes, 6)) S = self.stress(eps, theta, xi_n) return S, q, xi, C - - def UMAT_mixed( - self, - eps_idx: np.ndarray, - eps_n: np.ndarray, - deps: np.ndarray, - sig_bc: np.ndarray, - theta: float, - q_n: float, - xi_n: np.ndarray, - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float, np.ndarray]: - """ - Run the UMAT using partial eps-BC, e.g., uniaxial tension test. - - :param eps_idx: Indices of eps that are prescribed. If None, then all components of sigma are prescribed. - :type eps_idx: nd-array, dtype=int - :param eps_n: Strain at the beginning of the increment. - :type eps_n: nd-array, dtype=float, shape=[6] - :param deps: Strain increment. Only deps[eps_idx] is used. - :type deps: nd-array, dtype=float, shape=[6] - :param sig_bc: Stress at the end of the increment. Only non eps_idx entries are used. - :type sig_bc: nd-array, dtype=float, shape=[6] - :param theta: Temperature at the end of the time increment. - :type theta: float - :param q_n: Hardening variable at the beginning of the time increment. - :type q_n: float - :param xi_n: Reduced coefficients at the beginning of the time increment. - :type xi_n: nd-array, dtype=float - :return: **E** - - Full strain tensor at the end of the increment, i.e. - the entries not within eps_idx are set - :rtype: np.ndarray, dtype=float, shape=[6] - :return: **S** - - Stress at the end of the increment. Only S[eps_idx] is non-zero. - :rtype: np.ndarray, dtype=float - :return: **C** - - Stiffness tensor at the end of the increment - :rtype: `np.ndarray`, `dtype=float`, `shape=[6, 6]` - :return: **q** - - Hardening variable at the end of the time increment - :rtype: float - :return: **xi** - - Reduced coefficients at the end of the time increment - :rtype: `nd.ndarray`, `dtype=float` - """ - q = q_n - # partial strain update: - sig_idx = np.setdiff1d(np.arange(6), eps_idx) - - eps = eps_n.copy() - eps[eps_idx] += deps[eps_idx] - deps[sig_idx] = 0.0 - - self.tol_sig = 1e-7 * self.sig_y(theta, q_n, derivative=False) - err_sig = 1000.0 * self.tol_sig - it = 0 - it_max = 100 - - # eleminate thermo-elastic effects and approximate deps - sig = self.stress(eps_n, theta, xi_n) - # initial guess for eps asserting elastic behaviour - C = self.C - deps[sig_idx] += np.linalg.solve( - C[sig_idx][:, sig_idx], - -(sig[sig_idx] - sig_bc[sig_idx]) - C[sig_idx, :] @ deps, - ) - # compute the actual stress state - eps = eps_n + deps - sig, q, xi, C = self.solve(eps, deps, theta, q_n, xi_n) - err_sig = np.linalg.norm(sig[sig_idx]) - while (err_sig > self.tol_sig) and (it < it_max): - ddeps = -np.linalg.solve( - C[sig_idx][:, sig_idx], sig[sig_idx] - sig_bc[sig_idx] - ) - factor = 1.0 - deps[sig_idx] += factor * ddeps - eps = eps_n + deps - # err_eps = np.linalg.norm(ddeps) - sig, q, xi, C = self.solve(eps, deps, theta, q_n, xi_n) - err_sig = np.linalg.norm(sig[sig_idx]) - it = it + 1 - # print(f'it {it:3d} .. err_eps {err_eps:10.3e} (self.tol: {self.tol_eps:8.3e}) " - # + ".. err_sig {err_sig:10.3e} (self.tol: {self.tol_sig:8.3e}) ') - eps = eps_n + deps - # print(f'needed {it} iterations...') - return eps, sig, C, q, xi